SVD.h

Go to the documentation of this file.

#ifndef CEDRIC_ANALYSIS_SVD_H
#define CEDRIC_ANALYSIS_SVD_H

#include <algorithm>  // for min, max
#include <vector>
#include <math.h>

#include <assert.h>

#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/vector.hpp>

template<class ValueType>
class SVD
{
public:
        typedef boost::numeric::ublas::matrix<ValueType> Matrix;
        typedef boost::numeric::ublas::vector<ValueType> Vector;

        SVD( Matrix& matrix );

        const Vector& getSingularValues() const { return m_S; }
        
        const Matrix& getV() const { return m_V; }
        
protected:
        void svdcmp( Matrix& a, Vector& w, Matrix& v );

private:
        int m_rows, m_cols;
        Vector  m_S;
        Matrix  m_V;
};


//=====================================================================================================================

//              HELPER FUNCTIONS

//=====================================================================================================================

namespace {  // use anonymous namespace to avoid namespace pollution
        
// TODO: currently there is a bug in my nrutil replacement, so they can't be used for now!
#ifdef NRUTIL_REPLACEMENT
        template<class T>
        inline T SQR( T arg ) { return arg == 0.0 ? 0.0 : arg*arg; }

        template<class T>
        inline T SIGN( T a, T b ) { return b >= 0.0 ? fabs(a) : fabs(b); }
#else
        // #include "nrutil.h" --> extraced needed macros below
        static double maxarg1,maxarg2;
        #define FMAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1) > (maxarg2) ?\
        (maxarg1) : (maxarg2))
        static int iminarg1,iminarg2;
        #define IMIN(a,b) (iminarg1=(a),iminarg2=(b),(iminarg1) < (iminarg2) ?\
        (iminarg1) : (iminarg2))
        #define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
        static double sqrarg;
        #define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)
#endif


template<class T>
T pythag( T a, T b )
{
        T absa,absb;
        absa=fabs(a);
        absb=fabs(b);
        if (absa > absb) 
                return absa*sqrt(1.0+SQR(absb/absa));
        else 
                return (absb == 0.0 ? 0.0 : absb*sqrt(1.0+SQR(absa/absb)));
}

template<class Matrix,class Vector> void swap_columns( Matrix& M, int i, int j )
{
        Vector v = boost::numeric::ublas::column(M,i);
        boost::numeric::ublas::column(M,i) = boost::numeric::ublas::column(M,j);
        boost::numeric::ublas::column(M,j) = v;
}

} // -- end of anonymous namespace

//=====================================================================================================================

//              TEMPLATE IMPLEMENTATION

//=====================================================================================================================

#include <iostream>
template<class T>
SVD<T>::SVD( Matrix& matrix )
        : m_rows( matrix.size1() ), 
          m_cols( matrix.size2() ),
          m_S( m_cols ),
          m_V( m_cols, m_cols )
{
        svdcmp( matrix, m_S, m_V );
        
        // BOGO BUG: sorting of eigenvalues is not correct (-> hint that i f***ed up the code?)
        // use simple bubble sort because i only observed one fehlstand
        for( unsigned int run=0; run < m_S.size(); ++run )
        {
                bool ordered = true;
                
                for( unsigned int i=0; i < m_S.size()-1; ++i )
                {
                        if( m_S[i] < m_S[i+1] )
                        {
                                T tmp = m_S[i];
                                m_S[i] = m_S[i+1];
                                m_S[i+1] = tmp;
                                
                                swap_columns<Matrix,Vector>( matrix, i, i+1 );
                                swap_columns<Matrix,Vector>( m_V,    i, i+1 );
                                
                                ordered = false;
                        }
                }
                
                if( ordered ) break;            
        }
}

template<class T>
void SVD<T>::svdcmp( Matrix& a, Vector& w, Matrix& v )
{
        int m = a.size1();
        int n = a.size2();
        
        int flag,i,its,j,jj,k,l,nm;
        
        double anorm,c,f,g,h,s,scale,x,y,z;

        std::vector<double> rv1( n );

        g=scale=anorm=0.0;  // Householder reduction to bidiagonal form.

        for (i=1;i<=n;i++) 
        {
                l=i+1;
                rv1[i-1]=scale*g;
                g=s=scale=0.0;
                
                if (i <= m) 
                {
                        for (k=i;k<=m;k++) scale += fabs(a(k-1,i-1));
                        if (scale) 
                        {
                                for (k=i;k<=m;k++) 
                                {
                                        a(k-1,i-1) /= scale;
                                        s += a(k-1,i-1)*a(k-1,i-1);
                                }
                                
                                f=a(i-1,i-1);
                                g = -SIGN(sqrt(s),f);
                                h=f*g-s;
                                a(i-1,i-1)=f-g;
                                
                                for (j=l;j<=n;j++) 
                                {
                                        for (s=0.0,k=i;k<=m;k++) s += a(k-1,i-1)*a(k-1,j-1);
                                        f=s/h;
                                        for (k=i;k<=m;k++) a(k-1,j-1) += f*a(k-1,i-1);
                                }
                                
                                for (k=i;k<=m;k++) a(k-1,i-1) *= scale;
                        }
                }

                w[i-1]=scale *g;                
                g=s=scale=0.0;
                
                if (i <= m && i != n) 
                {
                        for (k=l;k<=n;k++) scale += fabs(a(i-1,k-1));
                        if (scale) 
                        {
                                for (k=l;k<=n;k++) 
                                {
                                        a(i-1,k-1) /= scale;
                                        s += a(i-1,k-1)*a(i-1,k-1);
                                }
                                
                                f=a(i-1,l-1);
                                g = -SIGN(sqrt(s),f);
                                h=f*g-s;
                                a(i-1,l-1)=f-g;
                                
                                for (k=l;k<=n;k++) rv1[k-1]=a(i-1,k-1)/h;
                                
                                for (j=l;j<=m;j++) 
                                {
                                        for (s=0.0,k=l;k<=n;k++) s += a(j-1,k-1)*a(i-1,k-1);
                                        for (k=l;k<=n;k++) a(j-1,k-1) += s*rv1[k-1];
                                }
                                
                                for (k=l;k<=n;k++) a(i-1,k-1) *= scale;
                        }
                }

          #ifdef NRUTIL_REPLACEMENT
                anorm=std::max(anorm,(fabs(w[i-1])+fabs(rv1[i-1])));
          #else
                anorm=FMAX(anorm,(fabs(w[i-1])+fabs(rv1[i])));          
          #endif
        }
        
        
        for (i=n;i>=1;i--)   // Accumulation of right-hand transformations.
        {
                if (i < n) 
                {
                        if (g) 
                        {
                                for (j=l;j<=n;j++)  // Double division to avoid possible under .o .
                                v(j-1,i-1)=(a(i-1,j-1)/a(i-1,l-1))/g;
                                
                                for (j=l;j<=n;j++) 
                                {
                                        for (s=0.0,k=l;k<=n;k++) s += a(i-1,k-1)*v(k-1,j-1);
                                        for (k=l;k<=n;k++) v(k-1,j-1) += s*v(k-1,i-1);
                                }
                        }
                        
                        for (j=l;j<=n;j++) v(i-1,j-1)=v(j-1,i-1)=0.0;
                }
                
                v(i-1,i-1)=1.0;
                g=rv1[i-1];
                l=i;
        }

  #ifdef NRUTIL_REPLACEMENT
        for (i=std::min(m,n);i>=1;i--)   // Accumulation of left-hand transformations.
  #else
        for (i=IMIN(m,n);i>=1;i--)   // Accumulation of left-hand transformations.
  #endif
        {
                l=i+1;
                g=w[i-1];
                
                for (j=l;j<=n;j++) a(i-1,j-1)=0.0;

                if (g) 
                {
                        g=1.0/g;
                        
                        for (j=l;j<=n;j++) 
                        {
                                for (s=0.0,k=l;k<=m;k++) s += a(k-1,i-1)*a(k-1,j-1);
                                f=(s/a(i-1,i-1))*g;
                                for (k=i;k<=m;k++) a(k-1,j-1) += f*a(k-1,i-1);
                        }
                        
                        for (j=i;j<=m;j++) a(j-1,i-1) *= g;
                                
                } 
                else for (j=i;j<=m;j++) a(j-1,i-1)=0.0;

                //++a(i-1,i-1);  //?????????
                a(i-1,i-1)+=1;
        }
        
        for (k=n;k>=1;k--)   // Diagonalization of the bidiagonal form:Loop over
        {
                
                for (its=1;its<=30;its++)   // singular values,and over allo ed iterations.
                {
                        flag=1;
                        for (l=k;l>=1;l--)   // Test for splitting.
                        {
                                nm=l-1;  // Note that rv1[1] is always zero.
                                if ((double)(fabs(rv1[l-1])+anorm) == anorm) 
                                {
                                        flag=0;
                                        break;
                                }

                                if ((double)(fabs(w[nm-1])+anorm) == anorm) break;
                        }
                        
                        if (flag) 
                        {
                                c=0.0;  // Cancellation of rv1[l],if l > 1
                                s=1.0;
                                
                                for (i=l;i<=k;i++) 
                                {
                                        f=s*rv1[i-1];
                                        rv1[i-1]=c*rv1[i-1];
                                        if ((double)(fabs(f)+anorm) == anorm) break;
                                        g=w[i-1];
                                        h=pythag(f,g);
                                        w[i-1]=h;
                                        h=1.0/h;
                                        c=g*h;
                                        s = -f*h;

                                        for (j=1;j<=m;j++) 
                                        {
                                                y=a(j-1,nm-1);
                                                z=a(j-1,i-1);
                                                a(j-1,nm-1)=y*c+z*s;
                                                a(j-1,i-1)=z*c-y*s;
                                        }
                                }
                        }
                        
                        z=w[k-1];
                        
                        if (l == k)   // Convergence.
                        {
                                if (z < 0.0)   // Singular value is made nonnegative.
                                {
                                        w[k-1] = -z;
                                        for (j=1;j<=n;j++) v(j-1,k-1) = -v(j-1,k-1);
                                }
                                break;
                        }

                        if (its == 30) assert(false);   // TODO: throw some kind of exception
                                //nrerror("no convergence in 30 svdcmp iterations");

                        x=w[l-1];  // Shift from bottom 2-by-2 minor.
                        nm=k-1;
                        y=w[nm-1];
                        g=rv1[nm-1];
                        h=rv1[k-1];
                        f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
                        g=pythag(f,1.0);
                        f=((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x;
                        c=s=1.0;  // Next QR transformation:
                        for (j=l;j<=nm;j++) 
                        {
                                i=j+1;
                                g=rv1[i-1];
                                y=w[i-1];
                                h=s*g;
                                g=c*g;
                                z=pythag(f,h);
                                rv1[j-1]=z;
                                c=f/z;
                                s=h/z;
                                f=x*c+g*s;
                                g = g*c-x*s;
                                h=y*s;
                                y *=c;
                                
                                for (jj=1;jj<=n;jj++) 
                                {
                                        x=v(jj-1,j-1);
                                        z=v(jj-1,i-1);
                                        v(jj-1,j-1)=x*c+z*s;
                                        v(jj-1,i-1)=z*c-x*s;
                                }
                                
                                z=pythag(f,h);
                                w[j-1]=z;  // Rotation can be arbitrary if z =0

                                if (z) 
                                {
                                        z=1.0/z;
                                        c=f*z;
                                        s=h*z;
                                }
                                
                                f=c*g+s*y;
                                x=c*y-s*g;
                                
                                for (jj=1;jj<=m;jj++) 
                                {
                                        y=a(jj-1,j-1);
                                        z=a(j-1,i-1);
                                        a(jj-1,j-1)=y*c+z*s;
                                        a(jj-1,i-1)=z*c-y*s;
                                }
                        }
                        
                        rv1[l-1]=0.0;
                        rv1[k-1]=f;
                        w[k-1]=x;
                }
        }       
}

#endif  // CEDRIC_ANALYSIS_SVD_H