![]() |
Eigen
3.3.3
|
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2009-2010 Benoit Jacob <jacob.benoit.1@gmail.com> 00005 // Copyright (C) 2013-2014 Gael Guennebaud <gael.guennebaud@inria.fr> 00006 // 00007 // This Source Code Form is subject to the terms of the Mozilla 00008 // Public License v. 2.0. If a copy of the MPL was not distributed 00009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00010 00011 #ifndef EIGEN_JACOBISVD_H 00012 #define EIGEN_JACOBISVD_H 00013 00014 namespace Eigen { 00015 00016 namespace internal { 00017 // forward declaration (needed by ICC) 00018 // the empty body is required by MSVC 00019 template<typename MatrixType, int QRPreconditioner, 00020 bool IsComplex = NumTraits<typename MatrixType::Scalar>::IsComplex> 00021 struct svd_precondition_2x2_block_to_be_real {}; 00022 00023 /*** QR preconditioners (R-SVD) 00024 *** 00025 *** Their role is to reduce the problem of computing the SVD to the case of a square matrix. 00026 *** This approach, known as R-SVD, is an optimization for rectangular-enough matrices, and is a requirement for 00027 *** JacobiSVD which by itself is only able to work on square matrices. 00028 ***/ 00029 00030 enum { PreconditionIfMoreColsThanRows, PreconditionIfMoreRowsThanCols }; 00031 00032 template<typename MatrixType, int QRPreconditioner, int Case> 00033 struct qr_preconditioner_should_do_anything 00034 { 00035 enum { a = MatrixType::RowsAtCompileTime != Dynamic && 00036 MatrixType::ColsAtCompileTime != Dynamic && 00037 MatrixType::ColsAtCompileTime <= MatrixType::RowsAtCompileTime, 00038 b = MatrixType::RowsAtCompileTime != Dynamic && 00039 MatrixType::ColsAtCompileTime != Dynamic && 00040 MatrixType::RowsAtCompileTime <= MatrixType::ColsAtCompileTime, 00041 ret = !( (QRPreconditioner == NoQRPreconditioner) || 00042 (Case == PreconditionIfMoreColsThanRows && bool(a)) || 00043 (Case == PreconditionIfMoreRowsThanCols && bool(b)) ) 00044 }; 00045 }; 00046 00047 template<typename MatrixType, int QRPreconditioner, int Case, 00048 bool DoAnything = qr_preconditioner_should_do_anything<MatrixType, QRPreconditioner, Case>::ret 00049 > struct qr_preconditioner_impl {}; 00050 00051 template<typename MatrixType, int QRPreconditioner, int Case> 00052 class qr_preconditioner_impl<MatrixType, QRPreconditioner, Case, false> 00053 { 00054 public: 00055 void allocate(const JacobiSVD<MatrixType, QRPreconditioner>&) {} 00056 bool run(JacobiSVD<MatrixType, QRPreconditioner>&, const MatrixType&) 00057 { 00058 return false; 00059 } 00060 }; 00061 00062 /*** preconditioner using FullPivHouseholderQR ***/ 00063 00064 template<typename MatrixType> 00065 class qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true> 00066 { 00067 public: 00068 typedef typename MatrixType::Scalar Scalar; 00069 enum 00070 { 00071 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 00072 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime 00073 }; 00074 typedef Matrix<Scalar, 1, RowsAtCompileTime, RowMajor, 1, MaxRowsAtCompileTime> WorkspaceType; 00075 00076 void allocate(const JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd) 00077 { 00078 if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols()) 00079 { 00080 m_qr.~QRType(); 00081 ::new (&m_qr) QRType(svd.rows(), svd.cols()); 00082 } 00083 if (svd.m_computeFullU) m_workspace.resize(svd.rows()); 00084 } 00085 00086 bool run(JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix) 00087 { 00088 if(matrix.rows() > matrix.cols()) 00089 { 00090 m_qr.compute(matrix); 00091 svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>(); 00092 if(svd.m_computeFullU) m_qr.matrixQ().evalTo(svd.m_matrixU, m_workspace); 00093 if(svd.computeV()) svd.m_matrixV = m_qr.colsPermutation(); 00094 return true; 00095 } 00096 return false; 00097 } 00098 private: 00099 typedef FullPivHouseholderQR<MatrixType> QRType; 00100 QRType m_qr; 00101 WorkspaceType m_workspace; 00102 }; 00103 00104 template<typename MatrixType> 00105 class qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true> 00106 { 00107 public: 00108 typedef typename MatrixType::Scalar Scalar; 00109 enum 00110 { 00111 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 00112 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 00113 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 00114 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, 00115 TrOptions = RowsAtCompileTime==1 ? (MatrixType::Options & ~(RowMajor)) 00116 : ColsAtCompileTime==1 ? (MatrixType::Options | RowMajor) 00117 : MatrixType::Options 00118 }; 00119 typedef Matrix<Scalar, ColsAtCompileTime, RowsAtCompileTime, TrOptions, MaxColsAtCompileTime, MaxRowsAtCompileTime> 00120 TransposeTypeWithSameStorageOrder; 00121 00122 void allocate(const JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd) 00123 { 00124 if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols()) 00125 { 00126 m_qr.~QRType(); 00127 ::new (&m_qr) QRType(svd.cols(), svd.rows()); 00128 } 00129 m_adjoint.resize(svd.cols(), svd.rows()); 00130 if (svd.m_computeFullV) m_workspace.resize(svd.cols()); 00131 } 00132 00133 bool run(JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix) 00134 { 00135 if(matrix.cols() > matrix.rows()) 00136 { 00137 m_adjoint = matrix.adjoint(); 00138 m_qr.compute(m_adjoint); 00139 svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint(); 00140 if(svd.m_computeFullV) m_qr.matrixQ().evalTo(svd.m_matrixV, m_workspace); 00141 if(svd.computeU()) svd.m_matrixU = m_qr.colsPermutation(); 00142 return true; 00143 } 00144 else return false; 00145 } 00146 private: 00147 typedef FullPivHouseholderQR<TransposeTypeWithSameStorageOrder> QRType; 00148 QRType m_qr; 00149 TransposeTypeWithSameStorageOrder m_adjoint; 00150 typename internal::plain_row_type<MatrixType>::type m_workspace; 00151 }; 00152 00153 /*** preconditioner using ColPivHouseholderQR ***/ 00154 00155 template<typename MatrixType> 00156 class qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true> 00157 { 00158 public: 00159 void allocate(const JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd) 00160 { 00161 if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols()) 00162 { 00163 m_qr.~QRType(); 00164 ::new (&m_qr) QRType(svd.rows(), svd.cols()); 00165 } 00166 if (svd.m_computeFullU) m_workspace.resize(svd.rows()); 00167 else if (svd.m_computeThinU) m_workspace.resize(svd.cols()); 00168 } 00169 00170 bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix) 00171 { 00172 if(matrix.rows() > matrix.cols()) 00173 { 00174 m_qr.compute(matrix); 00175 svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>(); 00176 if(svd.m_computeFullU) m_qr.householderQ().evalTo(svd.m_matrixU, m_workspace); 00177 else if(svd.m_computeThinU) 00178 { 00179 svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols()); 00180 m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_workspace); 00181 } 00182 if(svd.computeV()) svd.m_matrixV = m_qr.colsPermutation(); 00183 return true; 00184 } 00185 return false; 00186 } 00187 00188 private: 00189 typedef ColPivHouseholderQR<MatrixType> QRType; 00190 QRType m_qr; 00191 typename internal::plain_col_type<MatrixType>::type m_workspace; 00192 }; 00193 00194 template<typename MatrixType> 00195 class qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true> 00196 { 00197 public: 00198 typedef typename MatrixType::Scalar Scalar; 00199 enum 00200 { 00201 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 00202 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 00203 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 00204 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, 00205 TrOptions = RowsAtCompileTime==1 ? (MatrixType::Options & ~(RowMajor)) 00206 : ColsAtCompileTime==1 ? (MatrixType::Options | RowMajor) 00207 : MatrixType::Options 00208 }; 00209 00210 typedef Matrix<Scalar, ColsAtCompileTime, RowsAtCompileTime, TrOptions, MaxColsAtCompileTime, MaxRowsAtCompileTime> 00211 TransposeTypeWithSameStorageOrder; 00212 00213 void allocate(const JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd) 00214 { 00215 if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols()) 00216 { 00217 m_qr.~QRType(); 00218 ::new (&m_qr) QRType(svd.cols(), svd.rows()); 00219 } 00220 if (svd.m_computeFullV) m_workspace.resize(svd.cols()); 00221 else if (svd.m_computeThinV) m_workspace.resize(svd.rows()); 00222 m_adjoint.resize(svd.cols(), svd.rows()); 00223 } 00224 00225 bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix) 00226 { 00227 if(matrix.cols() > matrix.rows()) 00228 { 00229 m_adjoint = matrix.adjoint(); 00230 m_qr.compute(m_adjoint); 00231 00232 svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint(); 00233 if(svd.m_computeFullV) m_qr.householderQ().evalTo(svd.m_matrixV, m_workspace); 00234 else if(svd.m_computeThinV) 00235 { 00236 svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows()); 00237 m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_workspace); 00238 } 00239 if(svd.computeU()) svd.m_matrixU = m_qr.colsPermutation(); 00240 return true; 00241 } 00242 else return false; 00243 } 00244 00245 private: 00246 typedef ColPivHouseholderQR<TransposeTypeWithSameStorageOrder> QRType; 00247 QRType m_qr; 00248 TransposeTypeWithSameStorageOrder m_adjoint; 00249 typename internal::plain_row_type<MatrixType>::type m_workspace; 00250 }; 00251 00252 /*** preconditioner using HouseholderQR ***/ 00253 00254 template<typename MatrixType> 00255 class qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true> 00256 { 00257 public: 00258 void allocate(const JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd) 00259 { 00260 if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols()) 00261 { 00262 m_qr.~QRType(); 00263 ::new (&m_qr) QRType(svd.rows(), svd.cols()); 00264 } 00265 if (svd.m_computeFullU) m_workspace.resize(svd.rows()); 00266 else if (svd.m_computeThinU) m_workspace.resize(svd.cols()); 00267 } 00268 00269 bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix) 00270 { 00271 if(matrix.rows() > matrix.cols()) 00272 { 00273 m_qr.compute(matrix); 00274 svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>(); 00275 if(svd.m_computeFullU) m_qr.householderQ().evalTo(svd.m_matrixU, m_workspace); 00276 else if(svd.m_computeThinU) 00277 { 00278 svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols()); 00279 m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_workspace); 00280 } 00281 if(svd.computeV()) svd.m_matrixV.setIdentity(matrix.cols(), matrix.cols()); 00282 return true; 00283 } 00284 return false; 00285 } 00286 private: 00287 typedef HouseholderQR<MatrixType> QRType; 00288 QRType m_qr; 00289 typename internal::plain_col_type<MatrixType>::type m_workspace; 00290 }; 00291 00292 template<typename MatrixType> 00293 class qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true> 00294 { 00295 public: 00296 typedef typename MatrixType::Scalar Scalar; 00297 enum 00298 { 00299 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 00300 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 00301 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 00302 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, 00303 Options = MatrixType::Options 00304 }; 00305 00306 typedef Matrix<Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime> 00307 TransposeTypeWithSameStorageOrder; 00308 00309 void allocate(const JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd) 00310 { 00311 if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols()) 00312 { 00313 m_qr.~QRType(); 00314 ::new (&m_qr) QRType(svd.cols(), svd.rows()); 00315 } 00316 if (svd.m_computeFullV) m_workspace.resize(svd.cols()); 00317 else if (svd.m_computeThinV) m_workspace.resize(svd.rows()); 00318 m_adjoint.resize(svd.cols(), svd.rows()); 00319 } 00320 00321 bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix) 00322 { 00323 if(matrix.cols() > matrix.rows()) 00324 { 00325 m_adjoint = matrix.adjoint(); 00326 m_qr.compute(m_adjoint); 00327 00328 svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint(); 00329 if(svd.m_computeFullV) m_qr.householderQ().evalTo(svd.m_matrixV, m_workspace); 00330 else if(svd.m_computeThinV) 00331 { 00332 svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows()); 00333 m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_workspace); 00334 } 00335 if(svd.computeU()) svd.m_matrixU.setIdentity(matrix.rows(), matrix.rows()); 00336 return true; 00337 } 00338 else return false; 00339 } 00340 00341 private: 00342 typedef HouseholderQR<TransposeTypeWithSameStorageOrder> QRType; 00343 QRType m_qr; 00344 TransposeTypeWithSameStorageOrder m_adjoint; 00345 typename internal::plain_row_type<MatrixType>::type m_workspace; 00346 }; 00347 00348 /*** 2x2 SVD implementation 00349 *** 00350 *** JacobiSVD consists in performing a series of 2x2 SVD subproblems 00351 ***/ 00352 00353 template<typename MatrixType, int QRPreconditioner> 00354 struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, false> 00355 { 00356 typedef JacobiSVD<MatrixType, QRPreconditioner> SVD; 00357 typedef typename MatrixType::RealScalar RealScalar; 00358 static bool run(typename SVD::WorkMatrixType&, SVD&, Index, Index, RealScalar&) { return true; } 00359 }; 00360 00361 template<typename MatrixType, int QRPreconditioner> 00362 struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, true> 00363 { 00364 typedef JacobiSVD<MatrixType, QRPreconditioner> SVD; 00365 typedef typename MatrixType::Scalar Scalar; 00366 typedef typename MatrixType::RealScalar RealScalar; 00367 static bool run(typename SVD::WorkMatrixType& work_matrix, SVD& svd, Index p, Index q, RealScalar& maxDiagEntry) 00368 { 00369 using std::sqrt; 00370 using std::abs; 00371 Scalar z; 00372 JacobiRotation<Scalar> rot; 00373 RealScalar n = sqrt(numext::abs2(work_matrix.coeff(p,p)) + numext::abs2(work_matrix.coeff(q,p))); 00374 00375 const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)(); 00376 const RealScalar precision = NumTraits<Scalar>::epsilon(); 00377 00378 if(n==0) 00379 { 00380 // make sure first column is zero 00381 work_matrix.coeffRef(p,p) = work_matrix.coeffRef(q,p) = Scalar(0); 00382 00383 if(abs(numext::imag(work_matrix.coeff(p,q)))>considerAsZero) 00384 { 00385 // work_matrix.coeff(p,q) can be zero if work_matrix.coeff(q,p) is not zero but small enough to underflow when computing n 00386 z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q); 00387 work_matrix.row(p) *= z; 00388 if(svd.computeU()) svd.m_matrixU.col(p) *= conj(z); 00389 } 00390 if(abs(numext::imag(work_matrix.coeff(q,q)))>considerAsZero) 00391 { 00392 z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q); 00393 work_matrix.row(q) *= z; 00394 if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z); 00395 } 00396 // otherwise the second row is already zero, so we have nothing to do. 00397 } 00398 else 00399 { 00400 rot.c() = conj(work_matrix.coeff(p,p)) / n; 00401 rot.s() = work_matrix.coeff(q,p) / n; 00402 work_matrix.applyOnTheLeft(p,q,rot); 00403 if(svd.computeU()) svd.m_matrixU.applyOnTheRight(p,q,rot.adjoint()); 00404 if(abs(numext::imag(work_matrix.coeff(p,q)))>considerAsZero) 00405 { 00406 z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q); 00407 work_matrix.col(q) *= z; 00408 if(svd.computeV()) svd.m_matrixV.col(q) *= z; 00409 } 00410 if(abs(numext::imag(work_matrix.coeff(q,q)))>considerAsZero) 00411 { 00412 z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q); 00413 work_matrix.row(q) *= z; 00414 if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z); 00415 } 00416 } 00417 00418 // update largest diagonal entry 00419 maxDiagEntry = numext::maxi<RealScalar>(maxDiagEntry,numext::maxi<RealScalar>(abs(work_matrix.coeff(p,p)), abs(work_matrix.coeff(q,q)))); 00420 // and check whether the 2x2 block is already diagonal 00421 RealScalar threshold = numext::maxi<RealScalar>(considerAsZero, precision * maxDiagEntry); 00422 return abs(work_matrix.coeff(p,q))>threshold || abs(work_matrix.coeff(q,p)) > threshold; 00423 } 00424 }; 00425 00426 template<typename _MatrixType, int QRPreconditioner> 00427 struct traits<JacobiSVD<_MatrixType,QRPreconditioner> > 00428 { 00429 typedef _MatrixType MatrixType; 00430 }; 00431 00432 } // end namespace internal 00433 00487 template<typename _MatrixType, int QRPreconditioner> class JacobiSVD 00488 : public SVDBase<JacobiSVD<_MatrixType,QRPreconditioner> > 00489 { 00490 typedef SVDBase<JacobiSVD> Base; 00491 public: 00492 00493 typedef _MatrixType MatrixType; 00494 typedef typename MatrixType::Scalar Scalar; 00495 typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; 00496 enum { 00497 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 00498 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 00499 DiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_DYNAMIC(RowsAtCompileTime,ColsAtCompileTime), 00500 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 00501 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, 00502 MaxDiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(MaxRowsAtCompileTime,MaxColsAtCompileTime), 00503 MatrixOptions = MatrixType::Options 00504 }; 00505 00506 typedef typename Base::MatrixUType MatrixUType; 00507 typedef typename Base::MatrixVType MatrixVType; 00508 typedef typename Base::SingularValuesType SingularValuesType; 00509 00510 typedef typename internal::plain_row_type<MatrixType>::type RowType; 00511 typedef typename internal::plain_col_type<MatrixType>::type ColType; 00512 typedef Matrix<Scalar, DiagSizeAtCompileTime, DiagSizeAtCompileTime, 00513 MatrixOptions, MaxDiagSizeAtCompileTime, MaxDiagSizeAtCompileTime> 00514 WorkMatrixType; 00515 00521 JacobiSVD() 00522 {} 00523 00524 00531 JacobiSVD(Index rows, Index cols, unsigned int computationOptions = 0) 00532 { 00533 allocate(rows, cols, computationOptions); 00534 } 00535 00546 explicit JacobiSVD(const MatrixType& matrix, unsigned int computationOptions = 0) 00547 { 00548 compute(matrix, computationOptions); 00549 } 00550 00561 JacobiSVD& compute(const MatrixType& matrix, unsigned int computationOptions); 00562 00569 JacobiSVD& compute(const MatrixType& matrix) 00570 { 00571 return compute(matrix, m_computationOptions); 00572 } 00573 00574 using Base::computeU; 00575 using Base::computeV; 00576 using Base::rows; 00577 using Base::cols; 00578 using Base::rank; 00579 00580 private: 00581 void allocate(Index rows, Index cols, unsigned int computationOptions); 00582 00583 protected: 00584 using Base::m_matrixU; 00585 using Base::m_matrixV; 00586 using Base::m_singularValues; 00587 using Base::m_isInitialized; 00588 using Base::m_isAllocated; 00589 using Base::m_usePrescribedThreshold; 00590 using Base::m_computeFullU; 00591 using Base::m_computeThinU; 00592 using Base::m_computeFullV; 00593 using Base::m_computeThinV; 00594 using Base::m_computationOptions; 00595 using Base::m_nonzeroSingularValues; 00596 using Base::m_rows; 00597 using Base::m_cols; 00598 using Base::m_diagSize; 00599 using Base::m_prescribedThreshold; 00600 WorkMatrixType m_workMatrix; 00601 00602 template<typename __MatrixType, int _QRPreconditioner, bool _IsComplex> 00603 friend struct internal::svd_precondition_2x2_block_to_be_real; 00604 template<typename __MatrixType, int _QRPreconditioner, int _Case, bool _DoAnything> 00605 friend struct internal::qr_preconditioner_impl; 00606 00607 internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreColsThanRows> m_qr_precond_morecols; 00608 internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreRowsThanCols> m_qr_precond_morerows; 00609 MatrixType m_scaledMatrix; 00610 }; 00611 00612 template<typename MatrixType, int QRPreconditioner> 00613 void JacobiSVD<MatrixType, QRPreconditioner>::allocate(Index rows, Index cols, unsigned int computationOptions) 00614 { 00615 eigen_assert(rows >= 0 && cols >= 0); 00616 00617 if (m_isAllocated && 00618 rows == m_rows && 00619 cols == m_cols && 00620 computationOptions == m_computationOptions) 00621 { 00622 return; 00623 } 00624 00625 m_rows = rows; 00626 m_cols = cols; 00627 m_isInitialized = false; 00628 m_isAllocated = true; 00629 m_computationOptions = computationOptions; 00630 m_computeFullU = (computationOptions & ComputeFullU) != 0; 00631 m_computeThinU = (computationOptions & ComputeThinU) != 0; 00632 m_computeFullV = (computationOptions & ComputeFullV) != 0; 00633 m_computeThinV = (computationOptions & ComputeThinV) != 0; 00634 eigen_assert(!(m_computeFullU && m_computeThinU) && "JacobiSVD: you can't ask for both full and thin U"); 00635 eigen_assert(!(m_computeFullV && m_computeThinV) && "JacobiSVD: you can't ask for both full and thin V"); 00636 eigen_assert(EIGEN_IMPLIES(m_computeThinU || m_computeThinV, MatrixType::ColsAtCompileTime==Dynamic) && 00637 "JacobiSVD: thin U and V are only available when your matrix has a dynamic number of columns."); 00638 if (QRPreconditioner == FullPivHouseholderQRPreconditioner) 00639 { 00640 eigen_assert(!(m_computeThinU || m_computeThinV) && 00641 "JacobiSVD: can't compute thin U or thin V with the FullPivHouseholderQR preconditioner. " 00642 "Use the ColPivHouseholderQR preconditioner instead."); 00643 } 00644 m_diagSize = (std::min)(m_rows, m_cols); 00645 m_singularValues.resize(m_diagSize); 00646 if(RowsAtCompileTime==Dynamic) 00647 m_matrixU.resize(m_rows, m_computeFullU ? m_rows 00648 : m_computeThinU ? m_diagSize 00649 : 0); 00650 if(ColsAtCompileTime==Dynamic) 00651 m_matrixV.resize(m_cols, m_computeFullV ? m_cols 00652 : m_computeThinV ? m_diagSize 00653 : 0); 00654 m_workMatrix.resize(m_diagSize, m_diagSize); 00655 00656 if(m_cols>m_rows) m_qr_precond_morecols.allocate(*this); 00657 if(m_rows>m_cols) m_qr_precond_morerows.allocate(*this); 00658 if(m_rows!=m_cols) m_scaledMatrix.resize(rows,cols); 00659 } 00660 00661 template<typename MatrixType, int QRPreconditioner> 00662 JacobiSVD<MatrixType, QRPreconditioner>& 00663 JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsigned int computationOptions) 00664 { 00665 using std::abs; 00666 allocate(matrix.rows(), matrix.cols(), computationOptions); 00667 00668 // currently we stop when we reach precision 2*epsilon as the last bit of precision can require an unreasonable number of iterations, 00669 // only worsening the precision of U and V as we accumulate more rotations 00670 const RealScalar precision = RealScalar(2) * NumTraits<Scalar>::epsilon(); 00671 00672 // limit for denormal numbers to be considered zero in order to avoid infinite loops (see bug 286) 00673 const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)(); 00674 00675 // Scaling factor to reduce over/under-flows 00676 RealScalar scale = matrix.cwiseAbs().maxCoeff(); 00677 if(scale==RealScalar(0)) scale = RealScalar(1); 00678 00679 /*** step 1. The R-SVD step: we use a QR decomposition to reduce to the case of a square matrix */ 00680 00681 if(m_rows!=m_cols) 00682 { 00683 m_scaledMatrix = matrix / scale; 00684 m_qr_precond_morecols.run(*this, m_scaledMatrix); 00685 m_qr_precond_morerows.run(*this, m_scaledMatrix); 00686 } 00687 else 00688 { 00689 m_workMatrix = matrix.block(0,0,m_diagSize,m_diagSize) / scale; 00690 if(m_computeFullU) m_matrixU.setIdentity(m_rows,m_rows); 00691 if(m_computeThinU) m_matrixU.setIdentity(m_rows,m_diagSize); 00692 if(m_computeFullV) m_matrixV.setIdentity(m_cols,m_cols); 00693 if(m_computeThinV) m_matrixV.setIdentity(m_cols, m_diagSize); 00694 } 00695 00696 /*** step 2. The main Jacobi SVD iteration. ***/ 00697 RealScalar maxDiagEntry = m_workMatrix.cwiseAbs().diagonal().maxCoeff(); 00698 00699 bool finished = false; 00700 while(!finished) 00701 { 00702 finished = true; 00703 00704 // do a sweep: for all index pairs (p,q), perform SVD of the corresponding 2x2 sub-matrix 00705 00706 for(Index p = 1; p < m_diagSize; ++p) 00707 { 00708 for(Index q = 0; q < p; ++q) 00709 { 00710 // if this 2x2 sub-matrix is not diagonal already... 00711 // notice that this comparison will evaluate to false if any NaN is involved, ensuring that NaN's don't 00712 // keep us iterating forever. Similarly, small denormal numbers are considered zero. 00713 RealScalar threshold = numext::maxi<RealScalar>(considerAsZero, precision * maxDiagEntry); 00714 if(abs(m_workMatrix.coeff(p,q))>threshold || abs(m_workMatrix.coeff(q,p)) > threshold) 00715 { 00716 finished = false; 00717 // perform SVD decomposition of 2x2 sub-matrix corresponding to indices p,q to make it diagonal 00718 // the complex to real operation returns true if the updated 2x2 block is not already diagonal 00719 if(internal::svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner>::run(m_workMatrix, *this, p, q, maxDiagEntry)) 00720 { 00721 JacobiRotation<RealScalar> j_left, j_right; 00722 internal::real_2x2_jacobi_svd(m_workMatrix, p, q, &j_left, &j_right); 00723 00724 // accumulate resulting Jacobi rotations 00725 m_workMatrix.applyOnTheLeft(p,q,j_left); 00726 if(computeU()) m_matrixU.applyOnTheRight(p,q,j_left.transpose()); 00727 00728 m_workMatrix.applyOnTheRight(p,q,j_right); 00729 if(computeV()) m_matrixV.applyOnTheRight(p,q,j_right); 00730 00731 // keep track of the largest diagonal coefficient 00732 maxDiagEntry = numext::maxi<RealScalar>(maxDiagEntry,numext::maxi<RealScalar>(abs(m_workMatrix.coeff(p,p)), abs(m_workMatrix.coeff(q,q)))); 00733 } 00734 } 00735 } 00736 } 00737 } 00738 00739 /*** step 3. The work matrix is now diagonal, so ensure it's positive so its diagonal entries are the singular values ***/ 00740 00741 for(Index i = 0; i < m_diagSize; ++i) 00742 { 00743 // For a complex matrix, some diagonal coefficients might note have been 00744 // treated by svd_precondition_2x2_block_to_be_real, and the imaginary part 00745 // of some diagonal entry might not be null. 00746 if(NumTraits<Scalar>::IsComplex && abs(numext::imag(m_workMatrix.coeff(i,i)))>considerAsZero) 00747 { 00748 RealScalar a = abs(m_workMatrix.coeff(i,i)); 00749 m_singularValues.coeffRef(i) = abs(a); 00750 if(computeU()) m_matrixU.col(i) *= m_workMatrix.coeff(i,i)/a; 00751 } 00752 else 00753 { 00754 // m_workMatrix.coeff(i,i) is already real, no difficulty: 00755 RealScalar a = numext::real(m_workMatrix.coeff(i,i)); 00756 m_singularValues.coeffRef(i) = abs(a); 00757 if(computeU() && (a<RealScalar(0))) m_matrixU.col(i) = -m_matrixU.col(i); 00758 } 00759 } 00760 00761 m_singularValues *= scale; 00762 00763 /*** step 4. Sort singular values in descending order and compute the number of nonzero singular values ***/ 00764 00765 m_nonzeroSingularValues = m_diagSize; 00766 for(Index i = 0; i < m_diagSize; i++) 00767 { 00768 Index pos; 00769 RealScalar maxRemainingSingularValue = m_singularValues.tail(m_diagSize-i).maxCoeff(&pos); 00770 if(maxRemainingSingularValue == RealScalar(0)) 00771 { 00772 m_nonzeroSingularValues = i; 00773 break; 00774 } 00775 if(pos) 00776 { 00777 pos += i; 00778 std::swap(m_singularValues.coeffRef(i), m_singularValues.coeffRef(pos)); 00779 if(computeU()) m_matrixU.col(pos).swap(m_matrixU.col(i)); 00780 if(computeV()) m_matrixV.col(pos).swap(m_matrixV.col(i)); 00781 } 00782 } 00783 00784 m_isInitialized = true; 00785 return *this; 00786 } 00787 00795 template<typename Derived> 00796 JacobiSVD<typename MatrixBase<Derived>::PlainObject> 00797 MatrixBase<Derived>::jacobiSvd(unsigned int computationOptions) const 00798 { 00799 return JacobiSVD<PlainObject>(*this, computationOptions); 00800 } 00801 00802 } // end namespace Eigen 00803 00804 #endif // EIGEN_JACOBISVD_H