|| || || || || || || ||

---

MatPInv.h

Go to the documentation of this file.

#ifndef Impala_Core_Matrix_MatPInv_h
#define Impala_Core_Matrix_MatPInv_h

#include "Core/Matrix/MatFunc.h"
#include "Core/Matrix/Vecs.h"
#include "Core/Matrix/MatMul.h"
#include "Core/Matrix/VecFunc.h"
#include "Core/Matrix/MatSVD.h"
#include "Core/Matrix/MatTranspose.h"
#include "Core/Array/PixMax.h"
#include "Core/Array/SetVal.h"
#include "Core/Array/SetPart.h"

namespace Impala
{
namespace Core
{
namespace Matrix
{

template<class ArrayT>
inline ArrayT*
MatPInv(ArrayT* a, double tolerance=-1)
{
    int m = MatNrRow(a);
    int n = MatNrCol(a);
    
    if(n > m)
        return MatTranspose(MatPInv(MatTranspose(a), tolerance));

    Mat* u=0;
    Mat* s=0;
    Mat* v=0;

    MatSVD(a, u, s, v);
    Mat* res=  MatPInv(a,u,s,v,tolerance);

    delete u;
    delete s;
    delete v;

    return res;
}

template<class ArrayT>
inline ArrayT*
MatPInv(ArrayT* a,ArrayT* u,ArrayT* s,ArrayT* v, double tolerance=-1){
    
    Mat* res;
    VecScalarReal64* Singulars;

    int m = MatNrRow(a);
    int n = MatNrCol(a);

    int ssize = std::min(s->H(),s->W());

    //Store the singular values in a vector
    Singulars = VecCreate<VecScalarReal64>(ssize);
    SetVal(Singulars,0);
    for(int i=0; i<ssize;i++)
        Singulars->SetValue(s->Value(i,i),i,0);
    
    //If No tolerance is given, the default tolerance is:
    //eps*max(m,n)*max(Singulars)
    if(tolerance <= 0)
        tolerance = std::max(m,n) * PixMax(Singulars) * std::numeric_limits<double>::epsilon();

    //Find the number of singular values above tolerance
    int countAboveTolerance = 0;
    for(int i=0 ; i<ssize ; ++i)
        if(Singulars->Value(i,0) > tolerance)
            ++countAboveTolerance;

    if (countAboveTolerance == 0)
    {
        res = MatCreate<Mat>(n,m);
        SetVal(res, 0);
    }
    else
    {
        //SPlus is the square matrix of reciprocals of singular values that are
        //above tolerance - it is assumed that the Singular value matrix is
        //sorted from largest to smallest singular values.
        Mat* SPlus = MatCreate<Mat>(countAboveTolerance,countAboveTolerance);
        SetVal(SPlus, 0);
        for(int i=0 ; i<countAboveTolerance ; ++i)
            SPlus->SetValue(1/Singulars->Value(i,0),i,i);

        //We only need the first countAboveTolerance columns from V
        Mat* V =MatCreate<Mat>(n,countAboveTolerance);
        for(int i=0;i<countAboveTolerance;i++)
            for(int j=0;j<n;j++)
                V->SetValue(v->Value(i,j),i,j);

        //We only need the first countAboveTolerance rows from U'
        //Note that U is transposed:
        Mat* UT =MatCreate<Mat>(countAboveTolerance,m);
        for(int i=0;i<countAboveTolerance;i++)
            for(int j=0;j<m;j++)
                UT->SetValue(u->Value(i,j),j,i);

        Mat* SPlus_X_UT = MatMul(SPlus,UT);
        res = MatMul(V,SPlus_X_UT);
        //res = MatMul(V,MatMul(SPlus,UT));
        delete V;
        delete SPlus;
        delete SPlus_X_UT;
        delete UT;
    }

    delete Singulars;
    return res;
}

} // namespace Matrix
} // namespace Core
} // namespace Impala

#endif

---