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

---

MatKLM.h

Go to the documentation of this file.

#ifndef Impala_Core_Matrix_MatKLM_h
#define Impala_Core_Matrix_MatKLM_h

#include "Core/Matrix/MatCovariance.h"

namespace Impala
{
namespace Core
{
namespace Matrix
{


// KLM helping functions

/*
 Householder reduction of a real, symmetric matrix a[0..n-1][0..n-1]
 Output: a is replaced by the orthogonal matrix Q effecting the
 transformation. d[0..n-1] returns the diagonal elements of the
 tridiagonal matrix, and e[0..n-1] the off-diagonal elements, with e[0]=0;
 */
template<class MatrixT, class VectorT>
void
Tred2(MatrixT* a, VectorT* d, VectorT* e)
{
    int i,j,k,l;
    double scale,hh,h,g,f;

    int n = VecNrElem(d);

    for (i=n-1; i>0; i--)
    {
        l = i-1;
        h = scale = 0.0;
        if (l>0)
        {
            for (k=0; k<=l; k++)
                scale += fabs(*MatE(a, i, k));
            if (scale == 0.0)
                *VecE(e, i) = *MatE(a, i, l);
            else
            {
                for (k=0; k<=l; k++)
                {
                    *MatE(a, i, k) /= scale; //use scale a's for transformation
                    h += *MatE(a, i, k) * *MatE(a, i, k);
                }
                f = *MatE(a, i, l);
                g = (f >=0.0 ? -sqrt(h) : sqrt(h));
                *VecE(e, i) = scale*g;
                h -= f*g;
                *MatE(a, i, l) = f-g;
                f = 0.0;
                for (j=0; j<=l; j++)
                {
                    *MatE(a, j, i) = *MatE(a, i, j)/h;
                    g = 0.0;
                    for (k=0; k<=j; k++)
                        g += *MatE(a, j, k) * *MatE(a, i, k);
                    for (k=j+1; k<=l; k++)
                        g += *MatE(a, k, j) * *MatE(a, i, k);
                    *VecE(e, j) = g/h; //form element of p in temporarily unused element of e
                    f += *VecE(e, j) * *MatE(a, i, j);
                }
                hh = f/(h+h);
                for (j=0; j<=l; j++)
                {
                    f = *MatE(a, i, j);
                    *VecE(e, j) = g = *VecE(e, j)-hh*f;
                    for (k=0; k<=j; k++)
                        *MatE(a, j, k) -= (f * *VecE(e, k) 
                                             + g * *MatE(a, i, k));
                }
            }
        }
        else
            *VecE(e, i) = *MatE(a, i, l);
        *VecE(d, i) = h;
    }

    *VecE(d, 0) = 0.0;
    *VecE(e, 0) = 0.0;

    for (i=0; i<n; i++)
    {
        l = i-1;
        if (*VecE(d, i))
        {
            for (j=0; j<=l; j++)
            {
                g = 0.0;
                for (k=0; k<=l; k++)
                    g += *MatE(a, i, k) * *MatE(a, k, j);
                for (k=0; k<=l; k++)
                    *MatE(a, k, j) -= g * *MatE(a, k, i);
            }
        }
        *VecE(d, i) = *MatE(a, i, i);
        *MatE(a, i, i) = 1.0;
        for (j=0; j<=l; j++)
            *MatE(a, j, i) = *MatE(a, i, j) = 0.0;
    }
}

/*
 QL algorithm with implicit shifts, to determin the eigenvalues and eigenvectors
 of a real, symetric, tridiagonal matrix, of of a real, symmetric matrix previously
 reduced by tred2
 Input: d[0..n-1] contains the diagonal elements of the tridiagonal matrix
 e[0..n-1] inputs the subdiagonal elements of the tridiagonal matrix with e[0]
 arbitrary
 Output: d returns the eigenvalues, e is destroyed.
 If the eigenvectors of a tridiagonal matrix are desired, the matrix
 z[0..n-1][0..n-1] is input as the identity matrix. If the eigenvectors of
 a matrix that has been reduced by tred2 are required, then z is input as
 the matrix output by tred2.
 In either case, the kth column of z returns the normalized eigenvector
 corresponding to d[k]
 */
template<class MatrixT, class VectorT>
void
Tqli(VectorT* d, VectorT* e, MatrixT* z)
{

    int m,l,iter,i,k;
    double s,r,p,g,f,dd,c,b;

    int n = VecNrElem(e);

    for (i=1; i<n; i++)
        *VecE(e, i-1) = *VecE(e, i); //renumber
    *VecE(e, n-1) = 0.0;


    for (l=0; l<n; l++)
    {
        iter=0;
        do {
            for (m=l; m<n-1; m++)
            {
                dd = fabs(*VecE(d, m)) + fabs(*VecE(d, m+1));
                if ( (fabs(*VecE(e, m)) + dd) == dd)
                    break;
            }
            if (m != l)
            {
                if (iter++ == 30)
                    std::cout << "Tqli: Too many iterations" << std::endl;
                g = (*VecE(d, l+1) - *VecE(d, l)) / (2.0 * *VecE(e, l));
                r = sqrt(g*g + 1.0);//=pythag(g,1.0);
                g = *VecE(d, m) - *VecE(d, l) + *VecE(e, l)/(g+SIGN(r,g));
                s = c = 1.0;
                p = 0.0;
                for (i=m-1; i>=l; i--)
                {
                    f = s * *VecE(e, i);
                    b = c * *VecE(e, i);
                    *VecE(e, i+1) = (r=sqrt(f*f+g*g));//=(r=pythag(f,g));
                    if (r == 0.0)
                    {
                        *VecE(d, i+1) -= p;
                        *VecE(e, m) = 0.0;
                        break;
                    }
                    s = f/r;
                    c = g/r;
                    g = *VecE(d, i+1)-p;
                    r = (*VecE(d, i)-g)*s + 2.0*c*b;
                    *VecE(d, i+1) = g+(p=s*r);
                    g = c*r-b;

                    for (k=0; k<n; k++)
                    {
                        f = *MatE(z, k, i+1);
                        *MatE(z, k, i+1) = s * *MatE(z, k, i) + c*f;
                        *MatE(z, k, i) = c * *MatE(z, k, i) - s*f;
                    }
                }
                if (r==0.0 && i>=l)
                    continue;
                *VecE(d, l) -= p;
                *VecE(e, l) = g;
                *VecE(e, m) = 0.0;
            }
        } while (m != l);
    }

    //eigenvalues d[0..n-1], eigenvectors z[0..n-1][0..n-1], sort them now
    int j;
    for (i=0; i<n-1; i++)
    {
        p = *VecE(d, k=i);
        for (j=i+1; j<n; j++)
            if (*VecE(d, j) >= p)
                p = *VecE(d, k=j);
        if (k!=i)
        { //insert
            *VecE(d, k) = *VecE(d, i);
            *VecE(d, i) = p;
            //exchange corresponding eigenvectors
            for (j=0; j<n; j++)
            {
                p = *MatE(z, j, i);
                *MatE(z, j, i) = *MatE(z, j, k);
                *MatE(z, j, k) = p;
            }
        }
    }
}

/*
 *  Karhunen Loeve mapping (PCA)
 *  Principal Componant Analysis to reduce matrix to newDim dimensions
 */
template<class ArrayT>
inline ArrayT*
MatKLM(ArrayT* m, int newDim)
{
    typedef ArrayT Mat;
    typedef ArrayT Vec;

    int nr = MatNrRow(m);
    int nc = MatNrCol(m);

    if (newDim<1 || newDim>nc)
    {
        std::cout << "MatKLM: new dimensionality invalid, set to 1..." << std::endl;
        newDim = 1;
    }

    Mat* newfvec = MatCreate<Mat>(nr, newDim);
    Mat* a = MatCovariance(m, (Mat*) 0);
    Vec* d = VecCreate<Vec>(nc);
    Vec* e = VecCreate<Vec>(nc);

    Tred2(a, d, e);
    Tqli(d, e, a);

    // d[k] is an eigenvalue and a[][k] is a corresponding eigenvector

    for (int k=0; k<nr; k++)
    {
        for (int i=0; i<newDim; i++)
        {
            double vl = 0.0;
            for (int j=0; j<nc; j++)
                vl += *MatE(m, k, j) * *MatE(a, j, i);
            *MatE(newfvec, k, i) = vl;
        }
    }
    delete a;
    delete d;
    delete e;
    return newfvec;
}

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

#endif

---