SVD.h

Go to the documentation of this file.

// -*- c++ -*-

// Copyright (C) 2005,2009 Tom Drummond (twd20@cam.ac.uk)
//
// This file is part of the TooN Library.  This library is free
// software; you can redistribute it and/or modify it under the
// terms of the GNU General Public License as published by the
// Free Software Foundation; either version 2, or (at your option)
// any later version.

// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.

// You should have received a copy of the GNU General Public License along
// with this library; see the file COPYING.  If not, write to the Free
// Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307,
// USA.

// As a special exception, you may use this file as part of a free software
// library without restriction.  Specifically, if other files instantiate
// templates or use macros or inline functions from this file, or you compile
// this file and link it with other files to produce an executable, this
// file does not by itself cause the resulting executable to be covered by
// the GNU General Public License.  This exception does not however
// invalidate any other reasons why the executable file might be covered by
// the GNU General Public License.

#ifndef __SVD_H
#define __SVD_H

#include "TooN.h"
#include "lapack.h"


namespace TooN {

    // TODO - should this depend on precision?
static const double condition_no=1e9; // GK HACK TO GLOBAL







template<int Rows=Dynamic, int Cols=Rows, typename Precision=DefaultPrecision>
class SVD {
    // this is the size of the diagonal
    // NB works for semi-dynamic sizes because -1 < +ve ints
    static const int Min_Dim = Rows<Cols?Rows:Cols;
    
public:

    SVD() {}

    SVD(int rows, int cols)
        : my_copy(rows,cols),
          my_diagonal(std::min(rows,cols)),
          my_square(std::min(rows,cols), std::min(rows,cols))
    {}

    template <int R2, int C2, typename P2, typename B2>
    SVD(const Matrix<R2,C2,P2,B2>& m)
        : my_copy(m),
          my_diagonal(std::min(m.num_rows(),m.num_cols())),
          my_square(std::min(m.num_rows(),m.num_cols()),std::min(m.num_rows(),m.num_cols()))
    {
        do_compute();
    }

    template <int R2, int C2, typename P2, typename B2>
    void compute(const Matrix<R2,C2,P2,B2>& m){
        my_copy=m;
        do_compute();
    }
    
    private:
    void do_compute(){
        Precision* const a = my_copy.my_data;
        int lda = my_copy.num_cols();
        int m = my_copy.num_cols();
        int n = my_copy.num_rows();
        Precision* const uorvt = my_square.my_data;
        Precision* const s = my_diagonal.my_data;
        int ldu;
        int ldvt = lda;
        int LWORK;
        int INFO;
        char JOBU;
        char JOBVT;

        if(is_vertical()){ // u is a
            JOBU='O';
            JOBVT='S';
            ldu = lda;
        } else { // vt is a
            JOBU='S';
            JOBVT='O';
            ldu = my_square.num_cols();
        }

        Precision* wk;

        Precision size;
        LWORK = -1;

        // arguments are scrambled because we use rowmajor and lapack uses colmajor
        // thus u and vt play each other's roles.
        dgesvd_( &JOBVT, &JOBU, &m, &n, a, &lda, s, uorvt,
                 &ldvt, uorvt, &ldu, &size, &LWORK, &INFO);
    
        LWORK = (long int)(size);
        wk = new Precision[LWORK];

        dgesvd_( &JOBVT, &JOBU, &m, &n, a, &lda, s, uorvt,
                 &ldvt, uorvt, &ldu, wk, &LWORK, &INFO);
    
        delete[] wk;
    }
    
    bool is_vertical(){ 
        return (my_copy.num_rows() >= my_copy.num_cols()); 
    }

    int min_dim(){ return std::min(my_copy.num_rows(), my_copy.num_cols()); }
    
    public:

    template <int Rows2, int Cols2, typename P2, typename B2>
    Matrix<Cols,Cols2, typename Internal::MultiplyType<Precision,P2>::type >
    backsub(const Matrix<Rows2,Cols2,P2,B2>& rhs, const Precision condition=condition_no)
    {
        Vector<Min_Dim> inv_diag(min_dim());
        get_inv_diag(inv_diag,condition);
        return (get_VT().T() * diagmult(inv_diag, (get_U().T() * rhs)));
    }

    template <int Size, typename P2, typename B2>
    Vector<Cols, typename Internal::MultiplyType<Precision,P2>::type >
    backsub(const Vector<Size,P2,B2>& rhs, const Precision condition=condition_no)
    {
        Vector<Min_Dim> inv_diag(min_dim());
        get_inv_diag(inv_diag,condition);
        return (get_VT().T() * diagmult(inv_diag, (get_U().T() * rhs)));
    }

    Matrix<Cols,Rows> get_pinv(const Precision condition = condition_no){
        Vector<Min_Dim> inv_diag(min_dim());
        get_inv_diag(inv_diag,condition);
        return diagmult(get_VT().T(),inv_diag) * get_U().T();
    }

    Precision determinant() {
        Precision result = my_diagonal[0];
        for(int i=1; i<my_diagonal.size(); i++){
            result *= my_diagonal[i];
        }
        return result;
    }
    
    int rank(const Precision condition = condition_no) {
        if (my_diagonal[0] == 0) return 0;
        int result=1;
        for(int i=0; i<min_dim(); i++){
            if(my_diagonal[i] * condition <= my_diagonal[0]){
                result++;
            }
        }
        return result;
    }

    Matrix<Rows,Min_Dim,Precision,Reference::RowMajor> get_U(){
        if(is_vertical()){
            return Matrix<Rows,Min_Dim,Precision,Reference::RowMajor>
                (my_copy.my_data,my_copy.num_rows(),my_copy.num_cols());
        } else {
            return Matrix<Rows,Min_Dim,Precision,Reference::RowMajor>
                (my_square.my_data, my_square.num_rows(), my_square.num_cols());
        }
    }

    Vector<Min_Dim,Precision>& get_diagonal(){ return my_diagonal; }

    Matrix<Min_Dim,Cols,Precision,Reference::RowMajor> get_VT(){
        if(is_vertical()){
            return Matrix<Min_Dim,Cols,Precision,Reference::RowMajor>
                (my_square.my_data, my_square.num_rows(), my_square.num_cols());
        } else {
            return Matrix<Min_Dim,Cols,Precision,Reference::RowMajor>
                (my_copy.my_data,my_copy.num_rows(),my_copy.num_cols());
        }
    }

    void get_inv_diag(Vector<Min_Dim>& inv_diag, const Precision condition){
        for(int i=0; i<min_dim(); i++){
            if(my_diagonal[i] * condition <= my_diagonal[0]){
                inv_diag[i]=0;
            } else {
                inv_diag[i]=static_cast<Precision>(1)/my_diagonal[i];
            }
        }
    }

private:
    Matrix<Rows,Cols,Precision,RowMajor> my_copy;
    Vector<Min_Dim,Precision> my_diagonal;
    Matrix<Min_Dim,Min_Dim,Precision,RowMajor> my_square; // square matrix (U or V' depending on the shape of my_copy)
};






template<int Size, typename Precision>
struct SQSVD : public SVD<Size, Size, Precision> {
    SQSVD() {}
    SQSVD(int size) : SVD<Size,Size,Precision>(size, size) {}
    
    template <int R2, int C2, typename P2, typename B2>
    SQSVD(const Matrix<R2,C2,P2,B2>& m) : SVD<Size,Size,Precision>(m) {}
};


}


#endif