SVDBase.h (14743B)
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@gmail.com> 5 // Copyright (C) 2014 Gael Guennebaud <gael.guennebaud@inria.fr> 6 // 7 // Copyright (C) 2013 Gauthier Brun <brun.gauthier@gmail.com> 8 // Copyright (C) 2013 Nicolas Carre <nicolas.carre@ensimag.fr> 9 // Copyright (C) 2013 Jean Ceccato <jean.ceccato@ensimag.fr> 10 // Copyright (C) 2013 Pierre Zoppitelli <pierre.zoppitelli@ensimag.fr> 11 // 12 // This Source Code Form is subject to the terms of the Mozilla 13 // Public License v. 2.0. If a copy of the MPL was not distributed 14 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 15 16 #ifndef EIGEN_SVDBASE_H 17 #define EIGEN_SVDBASE_H 18 19 namespace Eigen { 20 21 namespace internal { 22 template<typename Derived> struct traits<SVDBase<Derived> > 23 : traits<Derived> 24 { 25 typedef MatrixXpr XprKind; 26 typedef SolverStorage StorageKind; 27 typedef int StorageIndex; 28 enum { Flags = 0 }; 29 }; 30 } 31 32 /** \ingroup SVD_Module 33 * 34 * 35 * \class SVDBase 36 * 37 * \brief Base class of SVD algorithms 38 * 39 * \tparam Derived the type of the actual SVD decomposition 40 * 41 * SVD decomposition consists in decomposing any n-by-p matrix \a A as a product 42 * \f[ A = U S V^* \f] 43 * 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; 44 * 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 45 * and right \em singular \em vectors of \a A respectively. 46 * 47 * Singular values are always sorted in decreasing order. 48 * 49 * 50 * 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 51 * 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 52 * 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, 53 * 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. 54 * 55 * The status of the computation can be retrived using the \a info() method. Unless \a info() returns \a Success, the results should be not 56 * considered well defined. 57 * 58 * If the input matrix has inf or nan coefficients, the result of the computation is undefined, and \a info() will return \a InvalidInput, but the computation is guaranteed to 59 * terminate in finite (and reasonable) time. 60 * \sa class BDCSVD, class JacobiSVD 61 */ 62 template<typename Derived> class SVDBase 63 : public SolverBase<SVDBase<Derived> > 64 { 65 public: 66 67 template<typename Derived_> 68 friend struct internal::solve_assertion; 69 70 typedef typename internal::traits<Derived>::MatrixType MatrixType; 71 typedef typename MatrixType::Scalar Scalar; 72 typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; 73 typedef typename Eigen::internal::traits<SVDBase>::StorageIndex StorageIndex; 74 typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3 75 enum { 76 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 77 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 78 DiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_DYNAMIC(RowsAtCompileTime,ColsAtCompileTime), 79 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 80 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, 81 MaxDiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(MaxRowsAtCompileTime,MaxColsAtCompileTime), 82 MatrixOptions = MatrixType::Options 83 }; 84 85 typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, MatrixOptions, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixUType; 86 typedef Matrix<Scalar, ColsAtCompileTime, ColsAtCompileTime, MatrixOptions, MaxColsAtCompileTime, MaxColsAtCompileTime> MatrixVType; 87 typedef typename internal::plain_diag_type<MatrixType, RealScalar>::type SingularValuesType; 88 89 Derived& derived() { return *static_cast<Derived*>(this); } 90 const Derived& derived() const { return *static_cast<const Derived*>(this); } 91 92 /** \returns the \a U matrix. 93 * 94 * For the SVD decomposition of a n-by-p matrix, letting \a m be the minimum of \a n and \a p, 95 * the U matrix is n-by-n if you asked for \link Eigen::ComputeFullU ComputeFullU \endlink, and is n-by-m if you asked for \link Eigen::ComputeThinU ComputeThinU \endlink. 96 * 97 * The \a m first columns of \a U are the left singular vectors of the matrix being decomposed. 98 * 99 * This method asserts that you asked for \a U to be computed. 100 */ 101 const MatrixUType& matrixU() const 102 { 103 _check_compute_assertions(); 104 eigen_assert(computeU() && "This SVD decomposition didn't compute U. Did you ask for it?"); 105 return m_matrixU; 106 } 107 108 /** \returns the \a V matrix. 109 * 110 * For the SVD decomposition of a n-by-p matrix, letting \a m be the minimum of \a n and \a p, 111 * the V matrix is p-by-p if you asked for \link Eigen::ComputeFullV ComputeFullV \endlink, and is p-by-m if you asked for \link Eigen::ComputeThinV ComputeThinV \endlink. 112 * 113 * The \a m first columns of \a V are the right singular vectors of the matrix being decomposed. 114 * 115 * This method asserts that you asked for \a V to be computed. 116 */ 117 const MatrixVType& matrixV() const 118 { 119 _check_compute_assertions(); 120 eigen_assert(computeV() && "This SVD decomposition didn't compute V. Did you ask for it?"); 121 return m_matrixV; 122 } 123 124 /** \returns the vector of singular values. 125 * 126 * For the SVD decomposition of a n-by-p matrix, letting \a m be the minimum of \a n and \a p, the 127 * returned vector has size \a m. Singular values are always sorted in decreasing order. 128 */ 129 const SingularValuesType& singularValues() const 130 { 131 _check_compute_assertions(); 132 return m_singularValues; 133 } 134 135 /** \returns the number of singular values that are not exactly 0 */ 136 Index nonzeroSingularValues() const 137 { 138 _check_compute_assertions(); 139 return m_nonzeroSingularValues; 140 } 141 142 /** \returns the rank of the matrix of which \c *this is the SVD. 143 * 144 * \note This method has to determine which singular values should be considered nonzero. 145 * For that, it uses the threshold value that you can control by calling 146 * setThreshold(const RealScalar&). 147 */ 148 inline Index rank() const 149 { 150 using std::abs; 151 _check_compute_assertions(); 152 if(m_singularValues.size()==0) return 0; 153 RealScalar premultiplied_threshold = numext::maxi<RealScalar>(m_singularValues.coeff(0) * threshold(), (std::numeric_limits<RealScalar>::min)()); 154 Index i = m_nonzeroSingularValues-1; 155 while(i>=0 && m_singularValues.coeff(i) < premultiplied_threshold) --i; 156 return i+1; 157 } 158 159 /** Allows to prescribe a threshold to be used by certain methods, such as rank() and solve(), 160 * which need to determine when singular values are to be considered nonzero. 161 * This is not used for the SVD decomposition itself. 162 * 163 * When it needs to get the threshold value, Eigen calls threshold(). 164 * The default is \c NumTraits<Scalar>::epsilon() 165 * 166 * \param threshold The new value to use as the threshold. 167 * 168 * A singular value will be considered nonzero if its value is strictly greater than 169 * \f$ \vert singular value \vert \leqslant threshold \times \vert max singular value \vert \f$. 170 * 171 * If you want to come back to the default behavior, call setThreshold(Default_t) 172 */ 173 Derived& setThreshold(const RealScalar& threshold) 174 { 175 m_usePrescribedThreshold = true; 176 m_prescribedThreshold = threshold; 177 return derived(); 178 } 179 180 /** Allows to come back to the default behavior, letting Eigen use its default formula for 181 * determining the threshold. 182 * 183 * You should pass the special object Eigen::Default as parameter here. 184 * \code svd.setThreshold(Eigen::Default); \endcode 185 * 186 * See the documentation of setThreshold(const RealScalar&). 187 */ 188 Derived& setThreshold(Default_t) 189 { 190 m_usePrescribedThreshold = false; 191 return derived(); 192 } 193 194 /** Returns the threshold that will be used by certain methods such as rank(). 195 * 196 * See the documentation of setThreshold(const RealScalar&). 197 */ 198 RealScalar threshold() const 199 { 200 eigen_assert(m_isInitialized || m_usePrescribedThreshold); 201 // this temporary is needed to workaround a MSVC issue 202 Index diagSize = (std::max<Index>)(1,m_diagSize); 203 return m_usePrescribedThreshold ? m_prescribedThreshold 204 : RealScalar(diagSize)*NumTraits<Scalar>::epsilon(); 205 } 206 207 /** \returns true if \a U (full or thin) is asked for in this SVD decomposition */ 208 inline bool computeU() const { return m_computeFullU || m_computeThinU; } 209 /** \returns true if \a V (full or thin) is asked for in this SVD decomposition */ 210 inline bool computeV() const { return m_computeFullV || m_computeThinV; } 211 212 inline Index rows() const { return m_rows; } 213 inline Index cols() const { return m_cols; } 214 215 #ifdef EIGEN_PARSED_BY_DOXYGEN 216 /** \returns a (least squares) solution of \f$ A x = b \f$ using the current SVD decomposition of A. 217 * 218 * \param b the right-hand-side of the equation to solve. 219 * 220 * \note Solving requires both U and V to be computed. Thin U and V are enough, there is no need for full U or V. 221 * 222 * \note SVD solving is implicitly least-squares. Thus, this method serves both purposes of exact solving and least-squares solving. 223 * In other words, the returned solution is guaranteed to minimize the Euclidean norm \f$ \Vert A x - b \Vert \f$. 224 */ 225 template<typename Rhs> 226 inline const Solve<Derived, Rhs> 227 solve(const MatrixBase<Rhs>& b) const; 228 #endif 229 230 231 /** \brief Reports whether previous computation was successful. 232 * 233 * \returns \c Success if computation was successful. 234 */ 235 EIGEN_DEVICE_FUNC 236 ComputationInfo info() const 237 { 238 eigen_assert(m_isInitialized && "SVD is not initialized."); 239 return m_info; 240 } 241 242 #ifndef EIGEN_PARSED_BY_DOXYGEN 243 template<typename RhsType, typename DstType> 244 void _solve_impl(const RhsType &rhs, DstType &dst) const; 245 246 template<bool Conjugate, typename RhsType, typename DstType> 247 void _solve_impl_transposed(const RhsType &rhs, DstType &dst) const; 248 #endif 249 250 protected: 251 252 static void check_template_parameters() 253 { 254 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar); 255 } 256 257 void _check_compute_assertions() const { 258 eigen_assert(m_isInitialized && "SVD is not initialized."); 259 } 260 261 template<bool Transpose_, typename Rhs> 262 void _check_solve_assertion(const Rhs& b) const { 263 EIGEN_ONLY_USED_FOR_DEBUG(b); 264 _check_compute_assertions(); 265 eigen_assert(computeU() && computeV() && "SVDBase::solve(): Both unitaries U and V are required to be computed (thin unitaries suffice)."); 266 eigen_assert((Transpose_?cols():rows())==b.rows() && "SVDBase::solve(): invalid number of rows of the right hand side matrix b"); 267 } 268 269 // return true if already allocated 270 bool allocate(Index rows, Index cols, unsigned int computationOptions) ; 271 272 MatrixUType m_matrixU; 273 MatrixVType m_matrixV; 274 SingularValuesType m_singularValues; 275 ComputationInfo m_info; 276 bool m_isInitialized, m_isAllocated, m_usePrescribedThreshold; 277 bool m_computeFullU, m_computeThinU; 278 bool m_computeFullV, m_computeThinV; 279 unsigned int m_computationOptions; 280 Index m_nonzeroSingularValues, m_rows, m_cols, m_diagSize; 281 RealScalar m_prescribedThreshold; 282 283 /** \brief Default Constructor. 284 * 285 * Default constructor of SVDBase 286 */ 287 SVDBase() 288 : m_info(Success), 289 m_isInitialized(false), 290 m_isAllocated(false), 291 m_usePrescribedThreshold(false), 292 m_computeFullU(false), 293 m_computeThinU(false), 294 m_computeFullV(false), 295 m_computeThinV(false), 296 m_computationOptions(0), 297 m_rows(-1), m_cols(-1), m_diagSize(0) 298 { 299 check_template_parameters(); 300 } 301 302 303 }; 304 305 #ifndef EIGEN_PARSED_BY_DOXYGEN 306 template<typename Derived> 307 template<typename RhsType, typename DstType> 308 void SVDBase<Derived>::_solve_impl(const RhsType &rhs, DstType &dst) const 309 { 310 // A = U S V^* 311 // So A^{-1} = V S^{-1} U^* 312 313 Matrix<typename RhsType::Scalar, Dynamic, RhsType::ColsAtCompileTime, 0, MatrixType::MaxRowsAtCompileTime, RhsType::MaxColsAtCompileTime> tmp; 314 Index l_rank = rank(); 315 tmp.noalias() = m_matrixU.leftCols(l_rank).adjoint() * rhs; 316 tmp = m_singularValues.head(l_rank).asDiagonal().inverse() * tmp; 317 dst = m_matrixV.leftCols(l_rank) * tmp; 318 } 319 320 template<typename Derived> 321 template<bool Conjugate, typename RhsType, typename DstType> 322 void SVDBase<Derived>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const 323 { 324 // A = U S V^* 325 // So A^{-*} = U S^{-1} V^* 326 // And A^{-T} = U_conj S^{-1} V^T 327 Matrix<typename RhsType::Scalar, Dynamic, RhsType::ColsAtCompileTime, 0, MatrixType::MaxRowsAtCompileTime, RhsType::MaxColsAtCompileTime> tmp; 328 Index l_rank = rank(); 329 330 tmp.noalias() = m_matrixV.leftCols(l_rank).transpose().template conjugateIf<Conjugate>() * rhs; 331 tmp = m_singularValues.head(l_rank).asDiagonal().inverse() * tmp; 332 dst = m_matrixU.template conjugateIf<!Conjugate>().leftCols(l_rank) * tmp; 333 } 334 #endif 335 336 template<typename MatrixType> 337 bool SVDBase<MatrixType>::allocate(Index rows, Index cols, unsigned int computationOptions) 338 { 339 eigen_assert(rows >= 0 && cols >= 0); 340 341 if (m_isAllocated && 342 rows == m_rows && 343 cols == m_cols && 344 computationOptions == m_computationOptions) 345 { 346 return true; 347 } 348 349 m_rows = rows; 350 m_cols = cols; 351 m_info = Success; 352 m_isInitialized = false; 353 m_isAllocated = true; 354 m_computationOptions = computationOptions; 355 m_computeFullU = (computationOptions & ComputeFullU) != 0; 356 m_computeThinU = (computationOptions & ComputeThinU) != 0; 357 m_computeFullV = (computationOptions & ComputeFullV) != 0; 358 m_computeThinV = (computationOptions & ComputeThinV) != 0; 359 eigen_assert(!(m_computeFullU && m_computeThinU) && "SVDBase: you can't ask for both full and thin U"); 360 eigen_assert(!(m_computeFullV && m_computeThinV) && "SVDBase: you can't ask for both full and thin V"); 361 eigen_assert(EIGEN_IMPLIES(m_computeThinU || m_computeThinV, MatrixType::ColsAtCompileTime==Dynamic) && 362 "SVDBase: thin U and V are only available when your matrix has a dynamic number of columns."); 363 364 m_diagSize = (std::min)(m_rows, m_cols); 365 m_singularValues.resize(m_diagSize); 366 if(RowsAtCompileTime==Dynamic) 367 m_matrixU.resize(m_rows, m_computeFullU ? m_rows : m_computeThinU ? m_diagSize : 0); 368 if(ColsAtCompileTime==Dynamic) 369 m_matrixV.resize(m_cols, m_computeFullV ? m_cols : m_computeThinV ? m_diagSize : 0); 370 371 return false; 372 } 373 374 }// end namespace 375 376 #endif // EIGEN_SVDBASE_H