SVD.h

// -*- c++ -*-

// Copyright (C) 2005,2009 Tom Drummond (twd20@cam.ac.uk)

//All rights reserved.
//
//Redistribution and use in source and binary forms, with or without
//modification, are permitted provided that the following conditions
//are met:
//1. Redistributions of source code must retain the above copyright
//    notice, this list of conditions and the following disclaimer.
//2. Redistributions in binary form must reproduce the above copyright
//   notice, this list of conditions and the following disclaimer in the
//   documentation and/or other materials provided with the distribution.
//
//THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND OTHER CONTRIBUTORS ``AS IS''
//AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
//IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
//ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR OTHER CONTRIBUTORS BE
//LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
//CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
//SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
//INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
//CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
//ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
//POSSIBILITY OF SUCH DAMAGE.

#ifndef __SVD_H
#define __SVD_H

#include <TooN/TooN.h>
#include <TooN/lapack.h>
#include <algorithm>

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.
        gesvd_( &JOBVT, &JOBU, &m, &n, a, &lda, s, uorvt,
                 &ldvt, uorvt, &ldu, &size, &LWORK, &INFO);
    
        LWORK = (long int)(size);
        wk = new Precision[LWORK];

        gesvd_( &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