SVD.hh

/* 
This file is part of the BIAS library (Basic ImageAlgorithmS).

Copyright (C) 2003-2009    (see file CONTACT for details)
  Multimediale Systeme der Informationsverarbeitung
  Institut fuer Informatik
  Christian-Albrechts-Universitaet Kiel


BIAS is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version.

BIAS 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 Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License
along with BIAS; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/


#ifndef _SVD_hh_
#define _SVD_hh_

#include "bias_config.h" 

#include <Base/Debug/Error.hh>
#include <Base/Debug/Debug.hh>

#include <Base/Math/Matrix.hh>
#include <Base/Math/Vector.hh>

#include <float.h>



// for netlib lapck svd routine tests:
//#include "netlib/dggsvd.c"

// double values below this threshold are treated as zero 
// (for nullspace computation)
//#define DEFAULT_DOUBLE_ZERO_THRESHOLD DBL_EPSILON
#define DEFAULT_DOUBLE_ZERO_THRESHOLD 1E-8

#define SVDSOLVE    0x00000001
#define D_SVD_RANK  0x00000002
#define D_SVD_SOLVE 0x00000004

namespace BIAS {
  /**   
        @class SVD
        @ingroup g_mathalgo
        @brief computes and holds the singular value decomposition 
        of a rectangular (not necessarily quadratic) Matrix A
  
        M = U S V't.
  
        The singular vectors vi and (right) singular values sigma_i satisfy the 
        equation:  A vi = sigma_i vi
  
        The Vector S_ stores the singular values of A in decreasing 
        order (smallest is last)  
    
        The columns of U  are the left singular vectors of A and form 
        an orthonormal basis of A
        which corresponds to the nonzero singular values. 
        while the columns of U corresponding to (nearly) zero values 
        are the nullspace. 
    
        This function is a base class and not a member in Matrix because: 
        - holding the U,S,V matrices allows repeated queries of the 
          same SVD results. 
        -avoid namespace clutter in the Matrix class. 
        There are many other 
        decompositions that might be of interest, 
        and adding them all would make for a very large Matrix class with 
        extra members which are irrelevant for the Matrix itself. 

        You can create the svd object either in relative or in absolute mode:
        if AbsoluteZeroThreshold is true, a singular value is compared
        to the threshold "absolutely", otherwise the ratio of the largest and
        the current singular value is compared against the threshold.
        Usually you want AbsoluteZeroThreshold = false, but for the sake of
        interface stability/downwards compatibility, the default is true
     
        @author Jan Woetzel
  **/
  class BIASMathAlgo_EXPORT SVD : public Debug {
  public:
    /** @brief solve the general eigenproblem of the rectangular matrix M 
        or with other words, make the U*S*V^T decomposition
        @param M matrix to decompose
        @param ZeroThreshold for nullspace determination, if a singular value
               is lower than ZeroThreshold, define the corresponding vector
               in V^T to be a null-vector (also sing.vals are treated this way)
        @param AbsoluteZeroThreshold if true, a singular value is compared
        to the threshold "absoultely", otherwise the ratio of the largest and
        the current singular value is compared against the threshold.
        Usually you want AbsoluteZeroThreshold = false, but for the sake of
        interface stability/downwards compatibility, the default is true
        @author Jan Woetzel **/
    SVD(const Matrix <double> & M, 
        double ZeroThreshold = DEFAULT_DOUBLE_ZERO_THRESHOLD,
        bool AbsoluteZeroThreshold = true);

    /** @brief creates an empty svd without automatically decomposing any
        matrix
        @param AbsoluteZeroThreshold if true, a singular value is compared
        to the threshold "absolutely", otherwise the ratio of the largest and
        the current singular value is compared against the threshold.
        Usually you want AbsoluteZeroThreshold = false, but for the sake of
        interface stability/downwards compatibility, the default is true  */
    inline SVD(bool AbsoluteZeroThreshold = true) { 
      Decomposed_ = false;
      ZeroThreshold_ = DEFAULT_DOUBLE_ZERO_THRESHOLD;
      AbsoluteZeroThreshold_ =  AbsoluteZeroThreshold;
    };

    inline virtual ~SVD(){};
  
    /** @brief set a new matrix and compute its decomposition.
        solves the general eigenproblem of the rectangular matrix M 
        @author Jan Woetzel 05/142002)*/
    int Compute(const Matrix<double>& M, 
                double ZeroThreshold=DEFAULT_DOUBLE_ZERO_THRESHOLD);

    /** use our naming convention */
    inline int compute(const Matrix<double>& M, 
                       double ZeroThreshold=DEFAULT_DOUBLE_ZERO_THRESHOLD)
    { return Compute(M, ZeroThreshold); };

#ifdef BIAS_HAVE_OPENCV
    /** \brief Use OpenCV for decomposition as Lapack frequently leads
               to crashes under windows 
       \attention This does not deliver the same results as Compute(). 
                  The columns of U and VT differ in sign.
       \author Ingo Schiller
    */
    int ComputeOpenCV(const Matrix<double>& M, 
                double ZeroThreshold=DEFAULT_DOUBLE_ZERO_THRESHOLD);
#endif
    /** solve the general (non-special) eigenvalue/eigenvector problem of 
        a general (non-symmetric) matrix M
        calls the extern liblapack routine dgesvd .
        computes the genral singular value decomposition with 
        M = U * Sigma(S) * VT 
        with 
        Sigma(S) is the qudratic, symmetric matrix which contains  
        the singular vales in the diagonal and zero and zero else.
        U, VT as described in the data mambers.
        @author Jan Woetzel (04/04/2002) */
    int General_Eigenproblem_GeneralMatrix_Lapack(const Matrix<double> &M);

    /** @brief return zerothresh currently used
        @author koeser */
    inline double GetZeroThreshold() const 
    { return ZeroThreshold_; }

    /** return S
        which is a vector of the singular values of A in descending order. 
        The diagonalelements of Sigma are the singular values of A.
        @author Jan Woetzel 04/04/2002 */
    inline const Vector<double>& GetS() const 
    { return S_; }
  
    /** return one singular value (which may be zero).
        the index runs from 0.. (size-1)
        @author Jan Woetzel 04/04/2002    **/
    inline const double GetSingularValue(int index ) const;
  
    /** return VT (=transposed(V))
        @author Jan Woetzel (04/04/2002) **/
    inline const Matrix<double>& GetVT() const 
    {  return VT_; };
  
    /** return V
        @author Jan Woetzel (04/04/2002) **/
    inline Matrix<double> GetV() const;

    /** return U
        U is a m x m orthogonal matrix
        @author Jan Woetzel (04/04/2002) **/
    inline const Matrix<double>& GetU() const 
    {  return U_; };

    /** return the length of the singular value vector
        inline const int size() const */ 
    inline int size() const { 
      if (U_.num_cols()>VT_.num_cols())
        return U_.num_cols();
      else 
        return VT_.num_cols();
    };

    /** @brief  return the dim of the nullspace.

        For rectangular matrices this is the number of singular values which
        are (about) zero. The nullspace is spaned by the singular vectors 
        corresponding to the zero singular values.
        woelk 4 2003: changed API, diffentiate between Left and Right Nullspace
        koeser 01/2005: rel/abs handled by flag:
        AbsoluteZeroThreshold_ determines whther we compare absolute or 
        relative singular values against the threshold
        @author Jan Woetzel 04/04/2002    **/
    inline const int NullspaceDim() const {
      if (AbsoluteZeroThreshold_) return AbsNullspaceDim();
      else return RelNullspaceDim();
    }

    /** @brief returns the dim of the matrix's right nullspace
        @author woelk 04 2003 */
    inline const int RightNullspaceDim() const {
      if (AbsoluteZeroThreshold_) return AbsRightNullspaceDim();
      else return RelRightNullspaceDim();
    }

    /** @brief returns the dim of the matrix's left nullspace
        @author woelk 04 2003 */
    inline const int LeftNullspaceDim() const {
      if (AbsoluteZeroThreshold_) return AbsLeftNullspaceDim();
      else return RelLeftNullspaceDim();
    }

    /** @brief returns dim of nullspace using the absolute threshold criterion
        
        For rectangular matrices this is the number of singular values which 
        are (about) zero. The nullspace is spaned by the singular vectors 
        corresponding to the zero singular values.
        woelk 4 2003: changed API, diffentiate between Left and Right Nullspace
        @author Jan Woetzel 04/04/2002    **/
    const int AbsNullspaceDim() const;

    /** @brief returns the dim of the matrix's right nullspace, using absolute 
        threshold criterion
        @author woelk 04 2003 */
    const int AbsRightNullspaceDim() const;

    /** @brief returns the dim of the matrix's left nullspace, using absolute 
        threshold criterion
        @author woelk 04 2003 */
    const int AbsLeftNullspaceDim() const;

    /** @brief compare singular values against greatest, not absolute
        @author Christian Buck
        changed criterion for insignificant singular values according to 
        hartley/zisserman from s_i<ZeroThreshold to s_i/s_0<ZeroThreshold.
        fallback to old criterion iff s_0=0. */
    const int RelNullspaceDim() const;
    const int RelRightNullspaceDim() const;
    const int RelLeftNullspaceDim() const;

    /** return one of the nullvectors.
        if last_offset=0 then the last nullvector 
        (corresponding to the smallest singular value) is returned,
        if last_offset=1 then the last but one nullvector is returned, 
        and so on.
        The last_offset must be in [0..NullspaceDim-1].
        @author Jan Woetzel(08/08/2002), JMF */
    inline Vector<double> GetNullvector(const int last_offset = 0);

    /** Returns one of the nullvectors in argument and true if 
        Nullvector exists.
        if last_offset=0 then the last nullvector (corresponding to the 
        smallest singular value) is returned,
        if last_offset=1 then the last but one nullvector is returned, 
        and so on.
        The last_offset must be in [0..NullspaceDim-1].
        @author Jan Woetzel (08/08/2002), JMF */
    inline bool GetNullvector(Vector<double>& NullVec,
                              const int last_offset = 0 );

    /** Return one of the left nullvectors.
        If last_offset=0 then the last nullvector 
        (corresponding to the smallest
        singular value) is returned,
        if last_offset=1 then the last but one nullvector is returned,
        and so on.
        The last_offset must be in [0..NullspaceDim-1] otherwise 0 is returned
        @author JMF/koeser **/
    inline bool GetLeftNullvector(Vector<double>&nv, const int last_offset=0);

    /** @brief same as above but returning vector */
    inline Vector<double> GetLeftNullvector( const int last_offset = 0 );

    // -- Solve the matrix-vector system M x = y, returning x.
    Vector<double> Solve(const Vector<double>& y)  const;
    /** use our naming convention */
    inline Vector<double> solve(const Vector<double>& y)const
    { return Solve(y); };

    /** returns pseudoinverse of A = U * S * V^T
        A^+ = V *  S^+ * U^T
        @author F Woelk */
    Matrix<double> Invert();
    
    /** returns pseudoinverse of A = U * S * V^T
        A^+ = V *  S^+ * U^T
        the first "rank" elements of S are inverted the others are set to zero 
        @author Christian Beder */
    Matrix<double> Invert(int rank);
    
    /** as above, but compute new svd for a */
    Matrix<double> Invert(Matrix<double> A);
#ifdef BIAS_HAVE_OPENCV
    /** \brief Calculates new inverse with OpenCV cvInvert 
        \attention This does not deliver the same results as Invert(). 
    */
    Matrix<double> InvertOpenCV(const Matrix<double>& A);
#endif
    /** returns the square root of a symmetric positive definite matrix M
        A = sqrt(M) = U * sqrt(S) * V_T, 
        such that A*A = M. 
        The result is only valid when the M *is* symmetric positive definite.
        @author woelk 01/2006 */
    Matrix<double> Sqrt();

    /** as above, but compute new svd for @param A */
    Matrix<double> Sqrt(const Matrix<double>& A);

    /** returns the square root of a symmetric positive definite matrix M
        A = sqrt(M) = U * sqrt(S), 
        such that A*A^T = M. 
        The result is only valid when the M *is* symmetric positive definite.
        @author woelk 01/2006 */
    Matrix<double> SqrtT();

    /** as above, but compute new svd for @param A */
    Matrix<double> SqrtT(const Matrix<double>& A);

    /// returns the rank of A_
    unsigned int Rank();


    /** solves the overdetermined linear system $AX=B$ with the unknown $X$,
        where $A$ is an $m x n$ matrix ($m>n$)and B is a vector of 
        size $m$.
        @author woelk 07/2004 */
    int Solve(Matrix<double>& A, Vector<double>& B, Vector<double>& X);

    /** Call this after Compute(),
        returns the eigenvalues of the matrix in ascending 
        order (smallest first), which are NOT the singular values!
        Works only, if the original Matrix was symmetric.
        @author grest, July 2006 */
    BIAS::Vector<double> GetEigenValues();


  protected:
    /// data members:
    /// original matrix (to be decomposed)
    Matrix<double> A_; 
    
    /// contains the singular values of A_ corresponding to the 
    /// i'th column in descending order.
    Vector<double> S_; 
    
    /// contains columnwise the left singular vectors  
    /// of A_ corresponding to the i'th singular value
    Matrix<double> U_; 
    
    /// contains the right singular vectors of A- in rows 
    /// [because VT is transpose(V) ]
    Matrix<double> VT_;

    /// determines whether we compare singular-values directly (absolute)
    /// to the zero threshold or whether we compare the ratio of the sing-value
    /// under inspection and the largest one to the threshold (relative)
    bool AbsoluteZeroThreshold_;

    /// values below this threshold are treated as zero
    double ZeroThreshold_; 
    
    /// flag for holding decomposed matrix
    bool  Decomposed_;

  }; // class SVD

  /////////////////////////////////////////////////////////////////
  /// inline implementations
  /////////////////////////////////////////////////////////////////

  inline const double SVD::GetSingularValue(int index) const 
  {
    BIASASSERT( index >=0 );
    BIASASSERT( index < size() );
    if (index < S_.size())
      return S_[index];
    else
      return 0.0;
  } 

  inline Matrix<double> SVD::GetV() const
  {
    Matrix<double> matV;
    //matV = transpose( VT_ ); //uncomment if error under UNIX
    matV = VT_.Transpose();
    return matV;
  }

  inline Vector<double> SVD::GetNullvector(const int last_offset) 
  {
    Vector<double> vec(VT_.num_cols(), 0.0);
    if (!GetNullvector(vec, last_offset)) {
      //BIASERR("No Nullvector for "<<last_offset);
    }
    return vec;
  }

  inline bool SVD::GetNullvector(Vector<double>& NullVec,
                                 const int last_offset  ) 
  {
    // create a new vector
    NullVec.newsize(VT_.num_cols());
    // fill the (null-)vector

    for (int col=0; col < VT_.num_cols(); col++) {
      NullVec[col] = VT_[ VT_.num_cols() -1 - last_offset ][col];
    }
    
// check if there is a nullvector with this last_offset_index
    if (last_offset > (NullspaceDim() -1)) {
      // index not OK
      return false;
    } else {
      // index OK      
      return true;
    }
  }

  inline Vector<double> SVD::GetLeftNullvector(const int last_offset) 
  {
    Vector<double> vec(VT_.num_rows(), 0.0);
    if (!GetLeftNullvector(vec, last_offset)) {
      BIASERR("No left nullvector for "<<last_offset);
    }
    return vec;
  }

  inline bool SVD::GetLeftNullvector(Vector<double>& vec, 
                                     const int last_offset) {
    // check if there is a nullvector with this last_offset_index
    if (last_offset > (NullspaceDim() -1)) {
      // index not OK
      BIASERR("for last_offset=" <<last_offset
              <<" with actual NullspaceDim= "<<NullspaceDim()
              <<" can no left nullvector be retrieved. aborting. ");      
      BIASERR(S_ << U_ << VT_);
      return false;
    }
    // index OK
    vec.newsize( U_.num_rows());
    // fill the (null-)vector
    for (int row=0; row < VT_.num_rows(); row++) {
      vec[row] = U_[row][ U_.num_rows() -1 - last_offset ];
    }
    return true;
  }


} // namespace BIAS
#endif // _SVD_hh_