Home | History | Annotate | Download | only in SVD 
      1 // This file is part of Eigen, a lightweight C++ template library
      2 // for linear algebra.
      3 //
      4 // Copyright (C) 2009-2010 Benoit Jacob <jacob.benoit.1 (at) gmail.com>
      5 //
      6 // Copyright (C) 2013 Gauthier Brun <brun.gauthier (at) gmail.com>
      7 // Copyright (C) 2013 Nicolas Carre <nicolas.carre (at) ensimag.fr>
      8 // Copyright (C) 2013 Jean Ceccato <jean.ceccato (at) ensimag.fr>
      9 // Copyright (C) 2013 Pierre Zoppitelli <pierre.zoppitelli (at) ensimag.fr>
     10 //
     11 // This Source Code Form is subject to the terms of the Mozilla
     12 // Public License v. 2.0. If a copy of the MPL was not distributed
     13 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
     14 
     15 #ifndef EIGEN_SVD_H
     16 #define EIGEN_SVD_H
     17 
     18 namespace Eigen {
     19 /** \ingroup SVD_Module
     20  *
     21  *
     22  * \class SVDBase
     23  *
     24  * \brief Mother class of SVD classes algorithms
     25  *
     26  * \param MatrixType the type of the matrix of which we are computing the SVD decomposition
     27  * SVD decomposition consists in decomposing any n-by-p matrix \a A as a product
     28  *   \f[ A = U S V^* \f]
     29  * where \a U is a n-by-n unitary, \a V is a p-by-p unitary, and \a S is a n-by-p real positive matrix which is zero outside of its main diagonal;
     30  * the diagonal entries of S are known as the \em singular \em values of \a A and the columns of \a U and \a V are known as the left
     31  * and right \em singular \em vectors of \a A respectively.
     32  *
     33  * Singular values are always sorted in decreasing order.
     34  *
     35  *
     36  * You can ask for only \em thin \a U or \a V to be computed, meaning the following. In case of a rectangular n-by-p matrix, letting \a m be the
     37  * smaller value among \a n and \a p, there are only \a m singular vectors; the remaining columns of \a U and \a V do not correspond to actual
     38  * singular vectors. Asking for \em thin \a U or \a V means asking for only their \a m first columns to be formed. So \a U is then a n-by-m matrix,
     39  * and \a V is then a p-by-m matrix. Notice that thin \a U and \a V are all you need for (least squares) solving.
     40  *
     41  * If the input matrix has inf or nan coefficients, the result of the computation is undefined, but the computation is guaranteed to
     42  * terminate in finite (and reasonable) time.
     43  * \sa MatrixBase::genericSvd()
     44  */
     45 template<typename _MatrixType>
     46 class SVDBase
     47 {
     48 
     49 public:
     50   typedef _MatrixType MatrixType;
     51   typedef typename MatrixType::Scalar Scalar;
     52   typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
     53   typedef typename MatrixType::Index Index;
     54   enum {
     55     RowsAtCompileTime = MatrixType::RowsAtCompileTime,
     56     ColsAtCompileTime = MatrixType::ColsAtCompileTime,
     57     DiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_DYNAMIC(RowsAtCompileTime,ColsAtCompileTime),
     58     MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
     59     MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
     60     MaxDiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(MaxRowsAtCompileTime,MaxColsAtCompileTime),
     61     MatrixOptions = MatrixType::Options
     62   };
     63 
     64   typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime,
     65 		 MatrixOptions, MaxRowsAtCompileTime, MaxRowsAtCompileTime>
     66   MatrixUType;
     67   typedef Matrix<Scalar, ColsAtCompileTime, ColsAtCompileTime,
     68 		 MatrixOptions, MaxColsAtCompileTime, MaxColsAtCompileTime>
     69   MatrixVType;
     70   typedef typename internal::plain_diag_type<MatrixType, RealScalar>::type SingularValuesType;
     71   typedef typename internal::plain_row_type<MatrixType>::type RowType;
     72   typedef typename internal::plain_col_type<MatrixType>::type ColType;
     73   typedef Matrix<Scalar, DiagSizeAtCompileTime, DiagSizeAtCompileTime,
     74 		 MatrixOptions, MaxDiagSizeAtCompileTime, MaxDiagSizeAtCompileTime>
     75   WorkMatrixType;
     76 
     77 
     78 
     79 
     80   /** \brief Method performing the decomposition of given matrix using custom options.
     81    *
     82    * \param matrix the matrix to decompose
     83    * \param computationOptions optional parameter allowing to specify if you want full or thin U or V unitaries to be computed.
     84    *                           By default, none is computed. This is a bit-field, the possible bits are #ComputeFullU, #ComputeThinU,
     85    *                           #ComputeFullV, #ComputeThinV.
     86    *
     87    * Thin unitaries are only available if your matrix type has a Dynamic number of columns (for example MatrixXf). They also are not
     88    * available with the (non-default) FullPivHouseholderQR preconditioner.
     89    */
     90   SVDBase& compute(const MatrixType& matrix, unsigned int computationOptions);
     91 
     92   /** \brief Method performing the decomposition of given matrix using current options.
     93    *
     94    * \param matrix the matrix to decompose
     95    *
     96    * This method uses the current \a computationOptions, as already passed to the constructor or to compute(const MatrixType&, unsigned int).
     97    */
     98   //virtual SVDBase& compute(const MatrixType& matrix) = 0;
     99   SVDBase& compute(const MatrixType& matrix);
    100 
    101   /** \returns the \a U matrix.
    102    *
    103    * For the SVDBase decomposition of a n-by-p matrix, letting \a m be the minimum of \a n and \a p,
    104    * the U matrix is n-by-n if you asked for #ComputeFullU, and is n-by-m if you asked for #ComputeThinU.
    105    *
    106    * The \a m first columns of \a U are the left singular vectors of the matrix being decomposed.
    107    *
    108    * This method asserts that you asked for \a U to be computed.
    109    */
    110   const MatrixUType& matrixU() const
    111   {
    112     eigen_assert(m_isInitialized && "SVD is not initialized.");
    113     eigen_assert(computeU() && "This SVD decomposition didn't compute U. Did you ask for it?");
    114     return m_matrixU;
    115   }
    116 
    117   /** \returns the \a V matrix.
    118    *
    119    * For the SVD decomposition of a n-by-p matrix, letting \a m be the minimum of \a n and \a p,
    120    * the V matrix is p-by-p if you asked for #ComputeFullV, and is p-by-m if you asked for ComputeThinV.
    121    *
    122    * The \a m first columns of \a V are the right singular vectors of the matrix being decomposed.
    123    *
    124    * This method asserts that you asked for \a V to be computed.
    125    */
    126   const MatrixVType& matrixV() const
    127   {
    128     eigen_assert(m_isInitialized && "SVD is not initialized.");
    129     eigen_assert(computeV() && "This SVD decomposition didn't compute V. Did you ask for it?");
    130     return m_matrixV;
    131   }
    132 
    133   /** \returns the vector of singular values.
    134    *
    135    * For the SVD decomposition of a n-by-p matrix, letting \a m be the minimum of \a n and \a p, the
    136    * returned vector has size \a m.  Singular values are always sorted in decreasing order.
    137    */
    138   const SingularValuesType& singularValues() const
    139   {
    140     eigen_assert(m_isInitialized && "SVD is not initialized.");
    141     return m_singularValues;
    142   }
    143 
    144 
    145 
    146   /** \returns the number of singular values that are not exactly 0 */
    147   Index nonzeroSingularValues() const
    148   {
    149     eigen_assert(m_isInitialized && "SVD is not initialized.");
    150     return m_nonzeroSingularValues;
    151   }
    152 
    153 
    154   /** \returns true if \a U (full or thin) is asked for in this SVD decomposition */
    155   inline bool computeU() const { return m_computeFullU || m_computeThinU; }
    156   /** \returns true if \a V (full or thin) is asked for in this SVD decomposition */
    157   inline bool computeV() const { return m_computeFullV || m_computeThinV; }
    158 
    159 
    160   inline Index rows() const { return m_rows; }
    161   inline Index cols() const { return m_cols; }
    162 
    163 
    164 protected:
    165   // return true if already allocated
    166   bool allocate(Index rows, Index cols, unsigned int computationOptions) ;
    167 
    168   MatrixUType m_matrixU;
    169   MatrixVType m_matrixV;
    170   SingularValuesType m_singularValues;
    171   bool m_isInitialized, m_isAllocated;
    172   bool m_computeFullU, m_computeThinU;
    173   bool m_computeFullV, m_computeThinV;
    174   unsigned int m_computationOptions;
    175   Index m_nonzeroSingularValues, m_rows, m_cols, m_diagSize;
    176 
    177 
    178   /** \brief Default Constructor.
    179    *
    180    * Default constructor of SVDBase
    181    */
    182   SVDBase()
    183     : m_isInitialized(false),
    184       m_isAllocated(false),
    185       m_computationOptions(0),
    186       m_rows(-1), m_cols(-1)
    187   {}
    188 
    189 
    190 };
    191 
    192 
    193 template<typename MatrixType>
    194 bool SVDBase<MatrixType>::allocate(Index rows, Index cols, unsigned int computationOptions)
    195 {
    196   eigen_assert(rows >= 0 && cols >= 0);
    197 
    198   if (m_isAllocated &&
    199       rows == m_rows &&
    200       cols == m_cols &&
    201       computationOptions == m_computationOptions)
    202   {
    203     return true;
    204   }
    205 
    206   m_rows = rows;
    207   m_cols = cols;
    208   m_isInitialized = false;
    209   m_isAllocated = true;
    210   m_computationOptions = computationOptions;
    211   m_computeFullU = (computationOptions & ComputeFullU) != 0;
    212   m_computeThinU = (computationOptions & ComputeThinU) != 0;
    213   m_computeFullV = (computationOptions & ComputeFullV) != 0;
    214   m_computeThinV = (computationOptions & ComputeThinV) != 0;
    215   eigen_assert(!(m_computeFullU && m_computeThinU) && "SVDBase: you can't ask for both full and thin U");
    216   eigen_assert(!(m_computeFullV && m_computeThinV) && "SVDBase: you can't ask for both full and thin V");
    217   eigen_assert(EIGEN_IMPLIES(m_computeThinU || m_computeThinV, MatrixType::ColsAtCompileTime==Dynamic) &&
    218 	       "SVDBase: thin U and V are only available when your matrix has a dynamic number of columns.");
    219 
    220   m_diagSize = (std::min)(m_rows, m_cols);
    221   m_singularValues.resize(m_diagSize);
    222   if(RowsAtCompileTime==Dynamic)
    223     m_matrixU.resize(m_rows, m_computeFullU ? m_rows
    224 		     : m_computeThinU ? m_diagSize
    225 		     : 0);
    226   if(ColsAtCompileTime==Dynamic)
    227     m_matrixV.resize(m_cols, m_computeFullV ? m_cols
    228 		     : m_computeThinV ? m_diagSize
    229 		     : 0);
    230 
    231   return false;
    232 }
    233 
    234 }// end namespace
    235 
    236 #endif // EIGEN_SVD_H
    237