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