$treeview $search $mathjax
Eigen-unsupported
3.2.5
$projectbrief
|
$projectbrief
|
$searchbox |
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 // 00006 // This Source Code Form is subject to the terms of the Mozilla 00007 // Public License v. 2.0. If a copy of the MPL was not distributed 00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00009 00010 #ifndef EIGEN_JACOBISVD_H 00011 #define EIGEN_JACOBISVD_H 00012 00013 namespace Eigen { 00014 00015 namespace internal { 00016 // forward declaration (needed by ICC) 00017 // the empty body is required by MSVC 00018 template<typename MatrixType, int QRPreconditioner, 00019 bool IsComplex = NumTraits<typename MatrixType::Scalar>::IsComplex> 00020 struct svd_precondition_2x2_block_to_be_real {}; 00021 00022 /*** QR preconditioners (R-SVD) 00023 *** 00024 *** Their role is to reduce the problem of computing the SVD to the case of a square matrix. 00025 *** This approach, known as R-SVD, is an optimization for rectangular-enough matrices, and is a requirement for 00026 *** JacobiSVD which by itself is only able to work on square matrices. 00027 ***/ 00028 00029 enum { PreconditionIfMoreColsThanRows, PreconditionIfMoreRowsThanCols }; 00030 00031 template<typename MatrixType, int QRPreconditioner, int Case> 00032 struct qr_preconditioner_should_do_anything 00033 { 00034 enum { a = MatrixType::RowsAtCompileTime != Dynamic && 00035 MatrixType::ColsAtCompileTime != Dynamic && 00036 MatrixType::ColsAtCompileTime <= MatrixType::RowsAtCompileTime, 00037 b = MatrixType::RowsAtCompileTime != Dynamic && 00038 MatrixType::ColsAtCompileTime != Dynamic && 00039 MatrixType::RowsAtCompileTime <= MatrixType::ColsAtCompileTime, 00040 ret = !( (QRPreconditioner == NoQRPreconditioner) || 00041 (Case == PreconditionIfMoreColsThanRows && bool(a)) || 00042 (Case == PreconditionIfMoreRowsThanCols && bool(b)) ) 00043 }; 00044 }; 00045 00046 template<typename MatrixType, int QRPreconditioner, int Case, 00047 bool DoAnything = qr_preconditioner_should_do_anything<MatrixType, QRPreconditioner, Case>::ret 00048 > struct qr_preconditioner_impl {}; 00049 00050 template<typename MatrixType, int QRPreconditioner, int Case> 00051 class qr_preconditioner_impl<MatrixType, QRPreconditioner, Case, false> 00052 { 00053 public: 00054 typedef typename MatrixType::Index Index; 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::Index Index; 00069 typedef typename MatrixType::Scalar Scalar; 00070 enum 00071 { 00072 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 00073 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime 00074 }; 00075 typedef Matrix<Scalar, 1, RowsAtCompileTime, RowMajor, 1, MaxRowsAtCompileTime> WorkspaceType; 00076 00077 void allocate(const JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd) 00078 { 00079 if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols()) 00080 { 00081 m_qr.~QRType(); 00082 ::new (&m_qr) QRType(svd.rows(), svd.cols()); 00083 } 00084 if (svd.m_computeFullU) m_workspace.resize(svd.rows()); 00085 } 00086 00087 bool run(JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix) 00088 { 00089 if(matrix.rows() > matrix.cols()) 00090 { 00091 m_qr.compute(matrix); 00092 svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>(); 00093 if(svd.m_computeFullU) m_qr.matrixQ().evalTo(svd.m_matrixU, m_workspace); 00094 if(svd.computeV()) svd.m_matrixV = m_qr.colsPermutation(); 00095 return true; 00096 } 00097 return false; 00098 } 00099 private: 00100 typedef FullPivHouseholderQR<MatrixType> QRType; 00101 QRType m_qr; 00102 WorkspaceType m_workspace; 00103 }; 00104 00105 template<typename MatrixType> 00106 class qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true> 00107 { 00108 public: 00109 typedef typename MatrixType::Index Index; 00110 typedef typename MatrixType::Scalar Scalar; 00111 enum 00112 { 00113 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 00114 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 00115 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 00116 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, 00117 Options = MatrixType::Options 00118 }; 00119 typedef Matrix<Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, 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 typedef typename MatrixType::Index Index; 00160 00161 void allocate(const JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd) 00162 { 00163 if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols()) 00164 { 00165 m_qr.~QRType(); 00166 ::new (&m_qr) QRType(svd.rows(), svd.cols()); 00167 } 00168 if (svd.m_computeFullU) m_workspace.resize(svd.rows()); 00169 else if (svd.m_computeThinU) m_workspace.resize(svd.cols()); 00170 } 00171 00172 bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix) 00173 { 00174 if(matrix.rows() > matrix.cols()) 00175 { 00176 m_qr.compute(matrix); 00177 svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>(); 00178 if(svd.m_computeFullU) m_qr.householderQ().evalTo(svd.m_matrixU, m_workspace); 00179 else if(svd.m_computeThinU) 00180 { 00181 svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols()); 00182 m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_workspace); 00183 } 00184 if(svd.computeV()) svd.m_matrixV = m_qr.colsPermutation(); 00185 return true; 00186 } 00187 return false; 00188 } 00189 00190 private: 00191 typedef ColPivHouseholderQR<MatrixType> QRType; 00192 QRType m_qr; 00193 typename internal::plain_col_type<MatrixType>::type m_workspace; 00194 }; 00195 00196 template<typename MatrixType> 00197 class qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true> 00198 { 00199 public: 00200 typedef typename MatrixType::Index Index; 00201 typedef typename MatrixType::Scalar Scalar; 00202 enum 00203 { 00204 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 00205 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 00206 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 00207 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, 00208 Options = MatrixType::Options 00209 }; 00210 00211 typedef Matrix<Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime> 00212 TransposeTypeWithSameStorageOrder; 00213 00214 void allocate(const JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd) 00215 { 00216 if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols()) 00217 { 00218 m_qr.~QRType(); 00219 ::new (&m_qr) QRType(svd.cols(), svd.rows()); 00220 } 00221 if (svd.m_computeFullV) m_workspace.resize(svd.cols()); 00222 else if (svd.m_computeThinV) m_workspace.resize(svd.rows()); 00223 m_adjoint.resize(svd.cols(), svd.rows()); 00224 } 00225 00226 bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix) 00227 { 00228 if(matrix.cols() > matrix.rows()) 00229 { 00230 m_adjoint = matrix.adjoint(); 00231 m_qr.compute(m_adjoint); 00232 00233 svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint(); 00234 if(svd.m_computeFullV) m_qr.householderQ().evalTo(svd.m_matrixV, m_workspace); 00235 else if(svd.m_computeThinV) 00236 { 00237 svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows()); 00238 m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_workspace); 00239 } 00240 if(svd.computeU()) svd.m_matrixU = m_qr.colsPermutation(); 00241 return true; 00242 } 00243 else return false; 00244 } 00245 00246 private: 00247 typedef ColPivHouseholderQR<TransposeTypeWithSameStorageOrder> QRType; 00248 QRType m_qr; 00249 TransposeTypeWithSameStorageOrder m_adjoint; 00250 typename internal::plain_row_type<MatrixType>::type m_workspace; 00251 }; 00252 00253 /*** preconditioner using HouseholderQR ***/ 00254 00255 template<typename MatrixType> 00256 class qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true> 00257 { 00258 public: 00259 typedef typename MatrixType::Index Index; 00260 00261 void allocate(const JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd) 00262 { 00263 if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols()) 00264 { 00265 m_qr.~QRType(); 00266 ::new (&m_qr) QRType(svd.rows(), svd.cols()); 00267 } 00268 if (svd.m_computeFullU) m_workspace.resize(svd.rows()); 00269 else if (svd.m_computeThinU) m_workspace.resize(svd.cols()); 00270 } 00271 00272 bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix) 00273 { 00274 if(matrix.rows() > matrix.cols()) 00275 { 00276 m_qr.compute(matrix); 00277 svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>(); 00278 if(svd.m_computeFullU) m_qr.householderQ().evalTo(svd.m_matrixU, m_workspace); 00279 else if(svd.m_computeThinU) 00280 { 00281 svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols()); 00282 m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_workspace); 00283 } 00284 if(svd.computeV()) svd.m_matrixV.setIdentity(matrix.cols(), matrix.cols()); 00285 return true; 00286 } 00287 return false; 00288 } 00289 private: 00290 typedef HouseholderQR<MatrixType> QRType; 00291 QRType m_qr; 00292 typename internal::plain_col_type<MatrixType>::type m_workspace; 00293 }; 00294 00295 template<typename MatrixType> 00296 class qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true> 00297 { 00298 public: 00299 typedef typename MatrixType::Index Index; 00300 typedef typename MatrixType::Scalar Scalar; 00301 enum 00302 { 00303 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 00304 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 00305 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 00306 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, 00307 Options = MatrixType::Options 00308 }; 00309 00310 typedef Matrix<Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime> 00311 TransposeTypeWithSameStorageOrder; 00312 00313 void allocate(const JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd) 00314 { 00315 if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols()) 00316 { 00317 m_qr.~QRType(); 00318 ::new (&m_qr) QRType(svd.cols(), svd.rows()); 00319 } 00320 if (svd.m_computeFullV) m_workspace.resize(svd.cols()); 00321 else if (svd.m_computeThinV) m_workspace.resize(svd.rows()); 00322 m_adjoint.resize(svd.cols(), svd.rows()); 00323 } 00324 00325 bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix) 00326 { 00327 if(matrix.cols() > matrix.rows()) 00328 { 00329 m_adjoint = matrix.adjoint(); 00330 m_qr.compute(m_adjoint); 00331 00332 svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint(); 00333 if(svd.m_computeFullV) m_qr.householderQ().evalTo(svd.m_matrixV, m_workspace); 00334 else if(svd.m_computeThinV) 00335 { 00336 svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows()); 00337 m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_workspace); 00338 } 00339 if(svd.computeU()) svd.m_matrixU.setIdentity(matrix.rows(), matrix.rows()); 00340 return true; 00341 } 00342 else return false; 00343 } 00344 00345 private: 00346 typedef HouseholderQR<TransposeTypeWithSameStorageOrder> QRType; 00347 QRType m_qr; 00348 TransposeTypeWithSameStorageOrder m_adjoint; 00349 typename internal::plain_row_type<MatrixType>::type m_workspace; 00350 }; 00351 00352 /*** 2x2 SVD implementation 00353 *** 00354 *** JacobiSVD consists in performing a series of 2x2 SVD subproblems 00355 ***/ 00356 00357 template<typename MatrixType, int QRPreconditioner> 00358 struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, false> 00359 { 00360 typedef JacobiSVD<MatrixType, QRPreconditioner> SVD; 00361 typedef typename SVD::Index Index; 00362 static void run(typename SVD::WorkMatrixType&, SVD&, Index, Index) {} 00363 }; 00364 00365 template<typename MatrixType, int QRPreconditioner> 00366 struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, true> 00367 { 00368 typedef JacobiSVD<MatrixType, QRPreconditioner> SVD; 00369 typedef typename MatrixType::Scalar Scalar; 00370 typedef typename MatrixType::RealScalar RealScalar; 00371 typedef typename SVD::Index Index; 00372 static void run(typename SVD::WorkMatrixType& work_matrix, SVD& svd, Index p, Index q) 00373 { 00374 using std::sqrt; 00375 Scalar z; 00376 JacobiRotation<Scalar> rot; 00377 RealScalar n = sqrt(numext::abs2(work_matrix.coeff(p,p)) + numext::abs2(work_matrix.coeff(q,p))); 00378 if(n==0) 00379 { 00380 z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q); 00381 work_matrix.row(p) *= z; 00382 if(svd.computeU()) svd.m_matrixU.col(p) *= conj(z); 00383 z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q); 00384 work_matrix.row(q) *= z; 00385 if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z); 00386 } 00387 else 00388 { 00389 rot.c() = conj(work_matrix.coeff(p,p)) / n; 00390 rot.s() = work_matrix.coeff(q,p) / n; 00391 work_matrix.applyOnTheLeft(p,q,rot); 00392 if(svd.computeU()) svd.m_matrixU.applyOnTheRight(p,q,rot.adjoint()); 00393 if(work_matrix.coeff(p,q) != Scalar(0)) 00394 { 00395 Scalar z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q); 00396 work_matrix.col(q) *= z; 00397 if(svd.computeV()) svd.m_matrixV.col(q) *= z; 00398 } 00399 if(work_matrix.coeff(q,q) != Scalar(0)) 00400 { 00401 z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q); 00402 work_matrix.row(q) *= z; 00403 if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z); 00404 } 00405 } 00406 } 00407 }; 00408 00409 template<typename MatrixType, typename RealScalar, typename Index> 00410 void real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q, 00411 JacobiRotation<RealScalar> *j_left, 00412 JacobiRotation<RealScalar> *j_right) 00413 { 00414 using std::sqrt; 00415 Matrix<RealScalar,2,2> m; 00416 m << numext::real(matrix.coeff(p,p)), numext::real(matrix.coeff(p,q)), 00417 numext::real(matrix.coeff(q,p)), numext::real(matrix.coeff(q,q)); 00418 JacobiRotation<RealScalar> rot1; 00419 RealScalar t = m.coeff(0,0) + m.coeff(1,1); 00420 RealScalar d = m.coeff(1,0) - m.coeff(0,1); 00421 if(t == RealScalar(0)) 00422 { 00423 rot1.c() = RealScalar(0); 00424 rot1.s() = d > RealScalar(0) ? RealScalar(1) : RealScalar(-1); 00425 } 00426 else 00427 { 00428 RealScalar u = d / t; 00429 rot1.c() = RealScalar(1) / sqrt(RealScalar(1) + numext::abs2(u)); 00430 rot1.s() = rot1.c() * u; 00431 } 00432 m.applyOnTheLeft(0,1,rot1); 00433 j_right->makeJacobi(m,0,1); 00434 *j_left = rot1 * j_right->transpose(); 00435 } 00436 00437 } // end namespace internal 00438 00492 template<typename _MatrixType, int QRPreconditioner> 00493 class JacobiSVD : public SVDBase<_MatrixType> 00494 { 00495 public: 00496 00497 typedef _MatrixType MatrixType; 00498 typedef typename MatrixType::Scalar Scalar; 00499 typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; 00500 typedef typename MatrixType::Index Index; 00501 enum { 00502 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 00503 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 00504 DiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_DYNAMIC(RowsAtCompileTime,ColsAtCompileTime), 00505 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 00506 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, 00507 MaxDiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(MaxRowsAtCompileTime,MaxColsAtCompileTime), 00508 MatrixOptions = MatrixType::Options 00509 }; 00510 00511 typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, 00512 MatrixOptions, MaxRowsAtCompileTime, MaxRowsAtCompileTime> 00513 MatrixUType; 00514 typedef Matrix<Scalar, ColsAtCompileTime, ColsAtCompileTime, 00515 MatrixOptions, MaxColsAtCompileTime, MaxColsAtCompileTime> 00516 MatrixVType; 00517 typedef typename internal::plain_diag_type<MatrixType, RealScalar>::type SingularValuesType; 00518 typedef typename internal::plain_row_type<MatrixType>::type RowType; 00519 typedef typename internal::plain_col_type<MatrixType>::type ColType; 00520 typedef Matrix<Scalar, DiagSizeAtCompileTime, DiagSizeAtCompileTime, 00521 MatrixOptions, MaxDiagSizeAtCompileTime, MaxDiagSizeAtCompileTime> 00522 WorkMatrixType; 00523 00529 JacobiSVD() 00530 : SVDBase<_MatrixType>::SVDBase() 00531 {} 00532 00533 00540 JacobiSVD(Index rows, Index cols, unsigned int computationOptions = 0) 00541 : SVDBase<_MatrixType>::SVDBase() 00542 { 00543 allocate(rows, cols, computationOptions); 00544 } 00545 00556 JacobiSVD(const MatrixType& matrix, unsigned int computationOptions = 0) 00557 : SVDBase<_MatrixType>::SVDBase() 00558 { 00559 compute(matrix, computationOptions); 00560 } 00561 00572 SVDBase<MatrixType>& compute(const MatrixType& matrix, unsigned int computationOptions); 00573 00580 SVDBase<MatrixType>& compute(const MatrixType& matrix) 00581 { 00582 return compute(matrix, this->m_computationOptions); 00583 } 00584 00594 template<typename Rhs> 00595 inline const internal::solve_retval<JacobiSVD, Rhs> 00596 solve(const MatrixBase<Rhs>& b) const 00597 { 00598 eigen_assert(this->m_isInitialized && "JacobiSVD is not initialized."); 00599 eigen_assert(SVDBase<MatrixType>::computeU() && SVDBase<MatrixType>::computeV() && "JacobiSVD::solve() requires both unitaries U and V to be computed (thin unitaries suffice)."); 00600 return internal::solve_retval<JacobiSVD, Rhs>(*this, b.derived()); 00601 } 00602 00603 00604 00605 private: 00606 void allocate(Index rows, Index cols, unsigned int computationOptions); 00607 00608 protected: 00609 WorkMatrixType m_workMatrix; 00610 00611 template<typename __MatrixType, int _QRPreconditioner, bool _IsComplex> 00612 friend struct internal::svd_precondition_2x2_block_to_be_real; 00613 template<typename __MatrixType, int _QRPreconditioner, int _Case, bool _DoAnything> 00614 friend struct internal::qr_preconditioner_impl; 00615 00616 internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreColsThanRows> m_qr_precond_morecols; 00617 internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreRowsThanCols> m_qr_precond_morerows; 00618 }; 00619 00620 template<typename MatrixType, int QRPreconditioner> 00621 void JacobiSVD<MatrixType, QRPreconditioner>::allocate(Index rows, Index cols, unsigned int computationOptions) 00622 { 00623 if (SVDBase<MatrixType>::allocate(rows, cols, computationOptions)) return; 00624 00625 if (QRPreconditioner == FullPivHouseholderQRPreconditioner) 00626 { 00627 eigen_assert(!(this->m_computeThinU || this->m_computeThinV) && 00628 "JacobiSVD: can't compute thin U or thin V with the FullPivHouseholderQR preconditioner. " 00629 "Use the ColPivHouseholderQR preconditioner instead."); 00630 } 00631 00632 m_workMatrix.resize(this->m_diagSize, this->m_diagSize); 00633 00634 if(this->m_cols>this->m_rows) m_qr_precond_morecols.allocate(*this); 00635 if(this->m_rows>this->m_cols) m_qr_precond_morerows.allocate(*this); 00636 } 00637 00638 template<typename MatrixType, int QRPreconditioner> 00639 SVDBase<MatrixType>& 00640 JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsigned int computationOptions) 00641 { 00642 using std::abs; 00643 allocate(matrix.rows(), matrix.cols(), computationOptions); 00644 00645 // currently we stop when we reach precision 2*epsilon as the last bit of precision can require an unreasonable number of iterations, 00646 // only worsening the precision of U and V as we accumulate more rotations 00647 const RealScalar precision = RealScalar(2) * NumTraits<Scalar>::epsilon(); 00648 00649 // limit for very small denormal numbers to be considered zero in order to avoid infinite loops (see bug 286) 00650 const RealScalar considerAsZero = RealScalar(2) * std::numeric_limits<RealScalar>::denorm_min(); 00651 00652 /*** step 1. The R-SVD step: we use a QR decomposition to reduce to the case of a square matrix */ 00653 00654 if(!m_qr_precond_morecols.run(*this, matrix) && !m_qr_precond_morerows.run(*this, matrix)) 00655 { 00656 m_workMatrix = matrix.block(0,0,this->m_diagSize,this->m_diagSize); 00657 if(this->m_computeFullU) this->m_matrixU.setIdentity(this->m_rows,this->m_rows); 00658 if(this->m_computeThinU) this->m_matrixU.setIdentity(this->m_rows,this->m_diagSize); 00659 if(this->m_computeFullV) this->m_matrixV.setIdentity(this->m_cols,this->m_cols); 00660 if(this->m_computeThinV) this->m_matrixV.setIdentity(this->m_cols, this->m_diagSize); 00661 } 00662 00663 /*** step 2. The main Jacobi SVD iteration. ***/ 00664 00665 bool finished = false; 00666 while(!finished) 00667 { 00668 finished = true; 00669 00670 // do a sweep: for all index pairs (p,q), perform SVD of the corresponding 2x2 sub-matrix 00671 00672 for(Index p = 1; p < this->m_diagSize; ++p) 00673 { 00674 for(Index q = 0; q < p; ++q) 00675 { 00676 // if this 2x2 sub-matrix is not diagonal already... 00677 // notice that this comparison will evaluate to false if any NaN is involved, ensuring that NaN's don't 00678 // keep us iterating forever. Similarly, small denormal numbers are considered zero. 00679 using std::max; 00680 RealScalar threshold = (max)(considerAsZero, precision * (max)(abs(m_workMatrix.coeff(p,p)), 00681 abs(m_workMatrix.coeff(q,q)))); 00682 if((max)(abs(m_workMatrix.coeff(p,q)),abs(m_workMatrix.coeff(q,p))) > threshold) 00683 { 00684 finished = false; 00685 00686 // perform SVD decomposition of 2x2 sub-matrix corresponding to indices p,q to make it diagonal 00687 internal::svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner>::run(m_workMatrix, *this, p, q); 00688 JacobiRotation<RealScalar> j_left, j_right; 00689 internal::real_2x2_jacobi_svd(m_workMatrix, p, q, &j_left, &j_right); 00690 00691 // accumulate resulting Jacobi rotations 00692 m_workMatrix.applyOnTheLeft(p,q,j_left); 00693 if(SVDBase<MatrixType>::computeU()) this->m_matrixU.applyOnTheRight(p,q,j_left.transpose()); 00694 00695 m_workMatrix.applyOnTheRight(p,q,j_right); 00696 if(SVDBase<MatrixType>::computeV()) this->m_matrixV.applyOnTheRight(p,q,j_right); 00697 } 00698 } 00699 } 00700 } 00701 00702 /*** step 3. The work matrix is now diagonal, so ensure it's positive so its diagonal entries are the singular values ***/ 00703 00704 for(Index i = 0; i < this->m_diagSize; ++i) 00705 { 00706 RealScalar a = abs(m_workMatrix.coeff(i,i)); 00707 this->m_singularValues.coeffRef(i) = a; 00708 if(SVDBase<MatrixType>::computeU() && (a!=RealScalar(0))) this->m_matrixU.col(i) *= this->m_workMatrix.coeff(i,i)/a; 00709 } 00710 00711 /*** step 4. Sort singular values in descending order and compute the number of nonzero singular values ***/ 00712 00713 this->m_nonzeroSingularValues = this->m_diagSize; 00714 for(Index i = 0; i < this->m_diagSize; i++) 00715 { 00716 Index pos; 00717 RealScalar maxRemainingSingularValue = this->m_singularValues.tail(this->m_diagSize-i).maxCoeff(&pos); 00718 if(maxRemainingSingularValue == RealScalar(0)) 00719 { 00720 this->m_nonzeroSingularValues = i; 00721 break; 00722 } 00723 if(pos) 00724 { 00725 pos += i; 00726 std::swap(this->m_singularValues.coeffRef(i), this->m_singularValues.coeffRef(pos)); 00727 if(SVDBase<MatrixType>::computeU()) this->m_matrixU.col(pos).swap(this->m_matrixU.col(i)); 00728 if(SVDBase<MatrixType>::computeV()) this->m_matrixV.col(pos).swap(this->m_matrixV.col(i)); 00729 } 00730 } 00731 00732 this->m_isInitialized = true; 00733 return *this; 00734 } 00735 00736 namespace internal { 00737 template<typename _MatrixType, int QRPreconditioner, typename Rhs> 00738 struct solve_retval<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs> 00739 : solve_retval_base<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs> 00740 { 00741 typedef JacobiSVD<_MatrixType, QRPreconditioner> JacobiSVDType; 00742 EIGEN_MAKE_SOLVE_HELPERS(JacobiSVDType,Rhs) 00743 00744 template<typename Dest> void evalTo(Dest& dst) const 00745 { 00746 eigen_assert(rhs().rows() == dec().rows()); 00747 00748 // A = U S V^* 00749 // So A^{-1} = V S^{-1} U^* 00750 00751 Index diagSize = (std::min)(dec().rows(), dec().cols()); 00752 typename JacobiSVDType::SingularValuesType invertedSingVals(diagSize); 00753 00754 Index nonzeroSingVals = dec().nonzeroSingularValues(); 00755 invertedSingVals.head(nonzeroSingVals) = dec().singularValues().head(nonzeroSingVals).array().inverse(); 00756 invertedSingVals.tail(diagSize - nonzeroSingVals).setZero(); 00757 00758 dst = dec().matrixV().leftCols(diagSize) 00759 * invertedSingVals.asDiagonal() 00760 * dec().matrixU().leftCols(diagSize).adjoint() 00761 * rhs(); 00762 } 00763 }; 00764 } // end namespace internal 00765 00773 template<typename Derived> 00774 JacobiSVD<typename MatrixBase<Derived>::PlainObject> 00775 MatrixBase<Derived>::jacobiSvd(unsigned int computationOptions) const 00776 { 00777 return JacobiSVD<PlainObject>(*this, computationOptions); 00778 } 00779 00780 } // end namespace Eigen 00781 00782 #endif // EIGEN_JACOBISVD_H