Eigne Matrix Class Library
Dependents: Eigen_test Odometry_test AttitudeEstimation_usingTicker MPU9250_Quaternion_Binary_Serial ... more
JacobiSVD.h
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 00379 if(n==0) 00380 { 00381 z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q); 00382 work_matrix.row(p) *= z; 00383 if(svd.computeU()) svd.m_matrixU.col(p) *= conj(z); 00384 if(work_matrix.coeff(q,q)!=Scalar(0)) 00385 { 00386 z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q); 00387 work_matrix.row(q) *= z; 00388 if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z); 00389 } 00390 // otherwise the second row is already zero, so we have nothing to do. 00391 } 00392 else 00393 { 00394 rot.c() = conj(work_matrix.coeff(p,p)) / n; 00395 rot.s() = work_matrix.coeff(q,p) / n; 00396 work_matrix.applyOnTheLeft(p,q,rot); 00397 if(svd.computeU()) svd.m_matrixU.applyOnTheRight(p,q,rot.adjoint()); 00398 if(work_matrix.coeff(p,q) != Scalar(0)) 00399 { 00400 Scalar z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q); 00401 work_matrix.col(q) *= z; 00402 if(svd.computeV()) svd.m_matrixV.col(q) *= z; 00403 } 00404 if(work_matrix.coeff(q,q) != Scalar(0)) 00405 { 00406 z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q); 00407 work_matrix.row(q) *= z; 00408 if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z); 00409 } 00410 } 00411 } 00412 }; 00413 00414 template<typename MatrixType, typename RealScalar, typename Index> 00415 void real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q, 00416 JacobiRotation<RealScalar> *j_left, 00417 JacobiRotation<RealScalar> *j_right) 00418 { 00419 using std::sqrt; 00420 using std::abs; 00421 Matrix<RealScalar,2,2> m; 00422 m << numext::real(matrix.coeff(p,p)), numext::real(matrix.coeff(p,q)), 00423 numext::real(matrix.coeff(q,p)), numext::real(matrix.coeff(q,q)); 00424 JacobiRotation<RealScalar> rot1; 00425 RealScalar t = m.coeff(0,0) + m.coeff(1,1); 00426 RealScalar d = m.coeff(1,0) - m.coeff(0,1); 00427 if(t == RealScalar(0)) 00428 { 00429 rot1.c() = RealScalar(0); 00430 rot1.s() = d > RealScalar(0) ? RealScalar(1) : RealScalar(-1); 00431 } 00432 else 00433 { 00434 RealScalar t2d2 = numext::hypot(t,d); 00435 rot1.c() = abs(t)/t2d2; 00436 rot1.s() = d/t2d2; 00437 if(t<RealScalar(0)) 00438 rot1.s() = -rot1.s(); 00439 } 00440 m.applyOnTheLeft(0,1,rot1); 00441 j_right->makeJacobi(m,0,1); 00442 *j_left = rot1 * j_right->transpose(); 00443 } 00444 00445 } // end namespace internal 00446 00447 /** \ingroup SVD_Module 00448 * 00449 * 00450 * \class JacobiSVD 00451 * 00452 * \brief Two-sided Jacobi SVD decomposition of a rectangular matrix 00453 * 00454 * \param MatrixType the type of the matrix of which we are computing the SVD decomposition 00455 * \param QRPreconditioner this optional parameter allows to specify the type of QR decomposition that will be used internally 00456 * for the R-SVD step for non-square matrices. See discussion of possible values below. 00457 * 00458 * SVD decomposition consists in decomposing any n-by-p matrix \a A as a product 00459 * \f[ A = U S V^* \f] 00460 * 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; 00461 * 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 00462 * and right \em singular \em vectors of \a A respectively. 00463 * 00464 * Singular values are always sorted in decreasing order. 00465 * 00466 * This JacobiSVD decomposition computes only the singular values by default. If you want \a U or \a V, you need to ask for them explicitly. 00467 * 00468 * 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 00469 * 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 00470 * 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, 00471 * 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. 00472 * 00473 * Here's an example demonstrating basic usage: 00474 * \include JacobiSVD_basic.cpp 00475 * Output: \verbinclude JacobiSVD_basic.out 00476 * 00477 * This JacobiSVD class is a two-sided Jacobi R-SVD decomposition, ensuring optimal reliability and accuracy. The downside is that it's slower than 00478 * bidiagonalizing SVD algorithms for large square matrices; however its complexity is still \f$ O(n^2p) \f$ where \a n is the smaller dimension and 00479 * \a p is the greater dimension, meaning that it is still of the same order of complexity as the faster bidiagonalizing R-SVD algorithms. 00480 * In particular, like any R-SVD, it takes advantage of non-squareness in that its complexity is only linear in the greater dimension. 00481 * 00482 * If the input matrix has inf or nan coefficients, the result of the computation is undefined, but the computation is guaranteed to 00483 * terminate in finite (and reasonable) time. 00484 * 00485 * The possible values for QRPreconditioner are: 00486 * \li ColPivHouseholderQRPreconditioner is the default. In practice it's very safe. It uses column-pivoting QR. 00487 * \li FullPivHouseholderQRPreconditioner, is the safest and slowest. It uses full-pivoting QR. 00488 * Contrary to other QRs, it doesn't allow computing thin unitaries. 00489 * \li HouseholderQRPreconditioner is the fastest, and less safe and accurate than the pivoting variants. It uses non-pivoting QR. 00490 * This is very similar in safety and accuracy to the bidiagonalization process used by bidiagonalizing SVD algorithms (since bidiagonalization 00491 * is inherently non-pivoting). However the resulting SVD is still more reliable than bidiagonalizing SVDs because the Jacobi-based iterarive 00492 * process is more reliable than the optimized bidiagonal SVD iterations. 00493 * \li NoQRPreconditioner allows not to use a QR preconditioner at all. This is useful if you know that you will only be computing 00494 * JacobiSVD decompositions of square matrices. Non-square matrices require a QR preconditioner. Using this option will result in 00495 * faster compilation and smaller executable code. It won't significantly speed up computation, since JacobiSVD is always checking 00496 * if QR preconditioning is needed before applying it anyway. 00497 * 00498 * \sa MatrixBase::jacobiSvd() 00499 */ 00500 template<typename _MatrixType, int QRPreconditioner> class JacobiSVD 00501 { 00502 public: 00503 00504 typedef _MatrixType MatrixType; 00505 typedef typename MatrixType::Scalar Scalar; 00506 typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; 00507 typedef typename MatrixType::Index Index; 00508 enum { 00509 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 00510 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 00511 DiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_DYNAMIC(RowsAtCompileTime,ColsAtCompileTime), 00512 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 00513 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, 00514 MaxDiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(MaxRowsAtCompileTime,MaxColsAtCompileTime), 00515 MatrixOptions = MatrixType::Options 00516 }; 00517 00518 typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, 00519 MatrixOptions, MaxRowsAtCompileTime, MaxRowsAtCompileTime> 00520 MatrixUType ; 00521 typedef Matrix<Scalar, ColsAtCompileTime, ColsAtCompileTime, 00522 MatrixOptions, MaxColsAtCompileTime, MaxColsAtCompileTime> 00523 MatrixVType ; 00524 typedef typename internal::plain_diag_type<MatrixType, RealScalar>::type SingularValuesType; 00525 typedef typename internal::plain_row_type<MatrixType>::type RowType; 00526 typedef typename internal::plain_col_type<MatrixType>::type ColType; 00527 typedef Matrix<Scalar, DiagSizeAtCompileTime, DiagSizeAtCompileTime, 00528 MatrixOptions, MaxDiagSizeAtCompileTime, MaxDiagSizeAtCompileTime> 00529 WorkMatrixType ; 00530 00531 /** \brief Default Constructor. 00532 * 00533 * The default constructor is useful in cases in which the user intends to 00534 * perform decompositions via JacobiSVD::compute(const MatrixType&). 00535 */ 00536 JacobiSVD() 00537 : m_isInitialized(false), 00538 m_isAllocated(false), 00539 m_usePrescribedThreshold(false), 00540 m_computationOptions(0), 00541 m_rows(-1), m_cols(-1), m_diagSize(0) 00542 {} 00543 00544 00545 /** \brief Default Constructor with memory preallocation 00546 * 00547 * Like the default constructor but with preallocation of the internal data 00548 * according to the specified problem size. 00549 * \sa JacobiSVD() 00550 */ 00551 JacobiSVD(Index rows, Index cols, unsigned int computationOptions = 0) 00552 : m_isInitialized(false), 00553 m_isAllocated(false), 00554 m_usePrescribedThreshold(false), 00555 m_computationOptions(0), 00556 m_rows(-1), m_cols(-1) 00557 { 00558 allocate(rows, cols, computationOptions); 00559 } 00560 00561 /** \brief Constructor performing the decomposition of given matrix. 00562 * 00563 * \param matrix the matrix to decompose 00564 * \param computationOptions optional parameter allowing to specify if you want full or thin U or V unitaries to be computed. 00565 * By default, none is computed. This is a bit-field, the possible bits are #ComputeFullU, #ComputeThinU, 00566 * #ComputeFullV, #ComputeThinV. 00567 * 00568 * Thin unitaries are only available if your matrix type has a Dynamic number of columns (for example MatrixXf). They also are not 00569 * available with the (non-default) FullPivHouseholderQR preconditioner. 00570 */ 00571 JacobiSVD(const MatrixType& matrix, unsigned int computationOptions = 0) 00572 : m_isInitialized(false), 00573 m_isAllocated(false), 00574 m_usePrescribedThreshold(false), 00575 m_computationOptions(0), 00576 m_rows(-1), m_cols(-1) 00577 { 00578 compute(matrix, computationOptions); 00579 } 00580 00581 /** \brief Method performing the decomposition of given matrix using custom options. 00582 * 00583 * \param matrix the matrix to decompose 00584 * \param computationOptions optional parameter allowing to specify if you want full or thin U or V unitaries to be computed. 00585 * By default, none is computed. This is a bit-field, the possible bits are #ComputeFullU, #ComputeThinU, 00586 * #ComputeFullV, #ComputeThinV. 00587 * 00588 * Thin unitaries are only available if your matrix type has a Dynamic number of columns (for example MatrixXf). They also are not 00589 * available with the (non-default) FullPivHouseholderQR preconditioner. 00590 */ 00591 JacobiSVD& compute(const MatrixType& matrix, unsigned int computationOptions); 00592 00593 /** \brief Method performing the decomposition of given matrix using current options. 00594 * 00595 * \param matrix the matrix to decompose 00596 * 00597 * This method uses the current \a computationOptions, as already passed to the constructor or to compute(const MatrixType&, unsigned int). 00598 */ 00599 JacobiSVD& compute(const MatrixType& matrix) 00600 { 00601 return compute(matrix, m_computationOptions); 00602 } 00603 00604 /** \returns the \a U matrix. 00605 * 00606 * For the SVD decomposition of a n-by-p matrix, letting \a m be the minimum of \a n and \a p, 00607 * the U matrix is n-by-n if you asked for #ComputeFullU, and is n-by-m if you asked for #ComputeThinU. 00608 * 00609 * The \a m first columns of \a U are the left singular vectors of the matrix being decomposed. 00610 * 00611 * This method asserts that you asked for \a U to be computed. 00612 */ 00613 const MatrixUType & matrixU () const 00614 { 00615 eigen_assert(m_isInitialized && "JacobiSVD is not initialized."); 00616 eigen_assert(computeU () && "This JacobiSVD decomposition didn't compute U. Did you ask for it?"); 00617 return m_matrixU; 00618 } 00619 00620 /** \returns the \a V matrix. 00621 * 00622 * For the SVD decomposition of a n-by-p matrix, letting \a m be the minimum of \a n and \a p, 00623 * the V matrix is p-by-p if you asked for #ComputeFullV, and is p-by-m if you asked for ComputeThinV. 00624 * 00625 * The \a m first columns of \a V are the right singular vectors of the matrix being decomposed. 00626 * 00627 * This method asserts that you asked for \a V to be computed. 00628 */ 00629 const MatrixVType & matrixV () const 00630 { 00631 eigen_assert(m_isInitialized && "JacobiSVD is not initialized."); 00632 eigen_assert(computeV () && "This JacobiSVD decomposition didn't compute V. Did you ask for it?"); 00633 return m_matrixV; 00634 } 00635 00636 /** \returns the vector of singular values. 00637 * 00638 * For the SVD decomposition of a n-by-p matrix, letting \a m be the minimum of \a n and \a p, the 00639 * returned vector has size \a m. Singular values are always sorted in decreasing order. 00640 */ 00641 const SingularValuesType& singularValues () const 00642 { 00643 eigen_assert(m_isInitialized && "JacobiSVD is not initialized."); 00644 return m_singularValues; 00645 } 00646 00647 /** \returns true if \a U (full or thin) is asked for in this SVD decomposition */ 00648 inline bool computeU () const { return m_computeFullU || m_computeThinU; } 00649 /** \returns true if \a V (full or thin) is asked for in this SVD decomposition */ 00650 inline bool computeV () const { return m_computeFullV || m_computeThinV; } 00651 00652 /** \returns a (least squares) solution of \f$ A x = b \f$ using the current SVD decomposition of A. 00653 * 00654 * \param b the right-hand-side of the equation to solve. 00655 * 00656 * \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. 00657 * 00658 * \note SVD solving is implicitly least-squares. Thus, this method serves both purposes of exact solving and least-squares solving. 00659 * In other words, the returned solution is guaranteed to minimize the Euclidean norm \f$ \Vert A x - b \Vert \f$. 00660 */ 00661 template<typename Rhs> 00662 inline const internal::solve_retval<JacobiSVD, Rhs> 00663 solve (const MatrixBase<Rhs>& b) const 00664 { 00665 eigen_assert(m_isInitialized && "JacobiSVD is not initialized."); 00666 eigen_assert(computeU () && computeV () && "JacobiSVD::solve() requires both unitaries U and V to be computed (thin unitaries suffice)."); 00667 return internal::solve_retval<JacobiSVD, Rhs>(*this, b.derived()); 00668 } 00669 00670 /** \returns the number of singular values that are not exactly 0 */ 00671 Index nonzeroSingularValues () const 00672 { 00673 eigen_assert(m_isInitialized && "JacobiSVD is not initialized."); 00674 return m_nonzeroSingularValues; 00675 } 00676 00677 /** \returns the rank of the matrix of which \c *this is the SVD. 00678 * 00679 * \note This method has to determine which singular values should be considered nonzero. 00680 * For that, it uses the threshold value that you can control by calling 00681 * setThreshold(const RealScalar&). 00682 */ 00683 inline Index rank () const 00684 { 00685 using std::abs; 00686 eigen_assert(m_isInitialized && "JacobiSVD is not initialized."); 00687 if(m_singularValues.size()==0) return 0; 00688 RealScalar premultiplied_threshold = m_singularValues.coeff(0) * threshold(); 00689 Index i = m_nonzeroSingularValues-1; 00690 while(i>=0 && m_singularValues.coeff(i) < premultiplied_threshold) --i; 00691 return i+1; 00692 } 00693 00694 /** Allows to prescribe a threshold to be used by certain methods, such as rank() and solve(), 00695 * which need to determine when singular values are to be considered nonzero. 00696 * This is not used for the SVD decomposition itself. 00697 * 00698 * When it needs to get the threshold value, Eigen calls threshold(). 00699 * The default is \c NumTraits<Scalar>::epsilon() 00700 * 00701 * \param threshold The new value to use as the threshold. 00702 * 00703 * A singular value will be considered nonzero if its value is strictly greater than 00704 * \f$ \vert singular value \vert \leqslant threshold \times \vert max singular value \vert \f$. 00705 * 00706 * If you want to come back to the default behavior, call setThreshold(Default_t) 00707 */ 00708 JacobiSVD& setThreshold(const RealScalar& threshold) 00709 { 00710 m_usePrescribedThreshold = true; 00711 m_prescribedThreshold = threshold; 00712 return *this; 00713 } 00714 00715 /** Allows to come back to the default behavior, letting Eigen use its default formula for 00716 * determining the threshold. 00717 * 00718 * You should pass the special object Eigen::Default as parameter here. 00719 * \code svd.setThreshold(Eigen::Default); \endcode 00720 * 00721 * See the documentation of setThreshold(const RealScalar&). 00722 */ 00723 JacobiSVD& setThreshold(Default_t) 00724 { 00725 m_usePrescribedThreshold = false; 00726 return *this; 00727 } 00728 00729 /** Returns the threshold that will be used by certain methods such as rank(). 00730 * 00731 * See the documentation of setThreshold(const RealScalar&). 00732 */ 00733 RealScalar threshold() const 00734 { 00735 eigen_assert(m_isInitialized || m_usePrescribedThreshold); 00736 return m_usePrescribedThreshold ? m_prescribedThreshold 00737 : (std::max<Index>)(1,m_diagSize)*NumTraits<Scalar>::epsilon(); 00738 } 00739 00740 inline Index rows() const { return m_rows; } 00741 inline Index cols() const { return m_cols; } 00742 00743 private: 00744 void allocate(Index rows, Index cols, unsigned int computationOptions); 00745 00746 static void check_template_parameters() 00747 { 00748 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar); 00749 } 00750 00751 protected: 00752 MatrixUType m_matrixU; 00753 MatrixVType m_matrixV; 00754 SingularValuesType m_singularValues; 00755 WorkMatrixType m_workMatrix; 00756 bool m_isInitialized, m_isAllocated, m_usePrescribedThreshold; 00757 bool m_computeFullU, m_computeThinU; 00758 bool m_computeFullV, m_computeThinV; 00759 unsigned int m_computationOptions; 00760 Index m_nonzeroSingularValues, m_rows, m_cols, m_diagSize; 00761 RealScalar m_prescribedThreshold; 00762 00763 template<typename __MatrixType, int _QRPreconditioner, bool _IsComplex> 00764 friend struct internal::svd_precondition_2x2_block_to_be_real; 00765 template<typename __MatrixType, int _QRPreconditioner, int _Case, bool _DoAnything> 00766 friend struct internal::qr_preconditioner_impl; 00767 00768 internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreColsThanRows> m_qr_precond_morecols; 00769 internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreRowsThanCols> m_qr_precond_morerows; 00770 MatrixType m_scaledMatrix; 00771 }; 00772 00773 template<typename MatrixType, int QRPreconditioner> 00774 void JacobiSVD<MatrixType, QRPreconditioner>::allocate(Index rows, Index cols, unsigned int computationOptions) 00775 { 00776 eigen_assert(rows >= 0 && cols >= 0); 00777 00778 if (m_isAllocated && 00779 rows == m_rows && 00780 cols == m_cols && 00781 computationOptions == m_computationOptions) 00782 { 00783 return; 00784 } 00785 00786 m_rows = rows; 00787 m_cols = cols; 00788 m_isInitialized = false; 00789 m_isAllocated = true; 00790 m_computationOptions = computationOptions; 00791 m_computeFullU = (computationOptions & ComputeFullU) != 0; 00792 m_computeThinU = (computationOptions & ComputeThinU) != 0; 00793 m_computeFullV = (computationOptions & ComputeFullV) != 0; 00794 m_computeThinV = (computationOptions & ComputeThinV) != 0; 00795 eigen_assert(!(m_computeFullU && m_computeThinU) && "JacobiSVD: you can't ask for both full and thin U"); 00796 eigen_assert(!(m_computeFullV && m_computeThinV) && "JacobiSVD: you can't ask for both full and thin V"); 00797 eigen_assert(EIGEN_IMPLIES(m_computeThinU || m_computeThinV, MatrixType::ColsAtCompileTime==Dynamic) && 00798 "JacobiSVD: thin U and V are only available when your matrix has a dynamic number of columns."); 00799 if (QRPreconditioner == FullPivHouseholderQRPreconditioner) 00800 { 00801 eigen_assert(!(m_computeThinU || m_computeThinV) && 00802 "JacobiSVD: can't compute thin U or thin V with the FullPivHouseholderQR preconditioner. " 00803 "Use the ColPivHouseholderQR preconditioner instead."); 00804 } 00805 m_diagSize = (std::min)(m_rows, m_cols); 00806 m_singularValues.resize(m_diagSize); 00807 if(RowsAtCompileTime==Dynamic) 00808 m_matrixU.resize(m_rows, m_computeFullU ? m_rows 00809 : m_computeThinU ? m_diagSize 00810 : 0); 00811 if(ColsAtCompileTime==Dynamic) 00812 m_matrixV.resize(m_cols, m_computeFullV ? m_cols 00813 : m_computeThinV ? m_diagSize 00814 : 0); 00815 m_workMatrix.resize(m_diagSize, m_diagSize); 00816 00817 if(m_cols>m_rows) m_qr_precond_morecols.allocate(*this); 00818 if(m_rows>m_cols) m_qr_precond_morerows.allocate(*this); 00819 if(m_rows!=m_cols) m_scaledMatrix.resize(rows,cols); 00820 } 00821 00822 template<typename MatrixType, int QRPreconditioner> 00823 JacobiSVD<MatrixType, QRPreconditioner>& 00824 JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsigned int computationOptions) 00825 { 00826 check_template_parameters(); 00827 00828 using std::abs; 00829 allocate(matrix.rows(), matrix.cols(), computationOptions); 00830 00831 // currently we stop when we reach precision 2*epsilon as the last bit of precision can require an unreasonable number of iterations, 00832 // only worsening the precision of U and V as we accumulate more rotations 00833 const RealScalar precision = RealScalar(2) * NumTraits<Scalar>::epsilon(); 00834 00835 // limit for very small denormal numbers to be considered zero in order to avoid infinite loops (see bug 286) 00836 const RealScalar considerAsZero = RealScalar(2) * std::numeric_limits<RealScalar>::denorm_min(); 00837 00838 // Scaling factor to reduce over/under-flows 00839 RealScalar scale = matrix.cwiseAbs().maxCoeff(); 00840 if(scale==RealScalar(0)) scale = RealScalar(1); 00841 00842 /*** step 1. The R-SVD step: we use a QR decomposition to reduce to the case of a square matrix */ 00843 00844 if(m_rows!=m_cols) 00845 { 00846 m_scaledMatrix = matrix / scale; 00847 m_qr_precond_morecols.run(*this, m_scaledMatrix); 00848 m_qr_precond_morerows.run(*this, m_scaledMatrix); 00849 } 00850 else 00851 { 00852 m_workMatrix = matrix.block(0,0,m_diagSize,m_diagSize) / scale; 00853 if(m_computeFullU) m_matrixU.setIdentity(m_rows,m_rows); 00854 if(m_computeThinU) m_matrixU.setIdentity(m_rows,m_diagSize); 00855 if(m_computeFullV) m_matrixV.setIdentity(m_cols,m_cols); 00856 if(m_computeThinV) m_matrixV.setIdentity(m_cols, m_diagSize); 00857 } 00858 00859 /*** step 2. The main Jacobi SVD iteration. ***/ 00860 00861 bool finished = false; 00862 while(!finished) 00863 { 00864 finished = true; 00865 00866 // do a sweep: for all index pairs (p,q), perform SVD of the corresponding 2x2 sub-matrix 00867 00868 for(Index p = 1; p < m_diagSize; ++p) 00869 { 00870 for(Index q = 0; q < p; ++q) 00871 { 00872 // if this 2x2 sub-matrix is not diagonal already... 00873 // notice that this comparison will evaluate to false if any NaN is involved, ensuring that NaN's don't 00874 // keep us iterating forever. Similarly, small denormal numbers are considered zero. 00875 using std::max; 00876 RealScalar threshold = (max)(considerAsZero, precision * (max)(abs(m_workMatrix.coeff(p,p)), 00877 abs(m_workMatrix.coeff(q,q)))); 00878 // We compare both values to threshold instead of calling max to be robust to NaN (See bug 791) 00879 if(abs(m_workMatrix.coeff(p,q))>threshold || abs(m_workMatrix.coeff(q,p)) > threshold) 00880 { 00881 finished = false; 00882 00883 // perform SVD decomposition of 2x2 sub-matrix corresponding to indices p,q to make it diagonal 00884 internal::svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner>::run(m_workMatrix, *this, p, q); 00885 JacobiRotation<RealScalar> j_left, j_right; 00886 internal::real_2x2_jacobi_svd(m_workMatrix, p, q, &j_left, &j_right); 00887 00888 // accumulate resulting Jacobi rotations 00889 m_workMatrix.applyOnTheLeft(p,q,j_left); 00890 if(computeU()) m_matrixU.applyOnTheRight(p,q,j_left.transpose()); 00891 00892 m_workMatrix.applyOnTheRight(p,q,j_right); 00893 if(computeV()) m_matrixV.applyOnTheRight(p,q,j_right); 00894 } 00895 } 00896 } 00897 } 00898 00899 /*** step 3. The work matrix is now diagonal, so ensure it's positive so its diagonal entries are the singular values ***/ 00900 00901 for(Index i = 0; i < m_diagSize; ++i) 00902 { 00903 RealScalar a = abs(m_workMatrix.coeff(i,i)); 00904 m_singularValues.coeffRef(i) = a; 00905 if(computeU() && (a!=RealScalar(0))) m_matrixU.col(i) *= m_workMatrix.coeff(i,i)/a; 00906 } 00907 00908 /*** step 4. Sort singular values in descending order and compute the number of nonzero singular values ***/ 00909 00910 m_nonzeroSingularValues = m_diagSize; 00911 for(Index i = 0; i < m_diagSize; i++) 00912 { 00913 Index pos; 00914 RealScalar maxRemainingSingularValue = m_singularValues.tail(m_diagSize-i).maxCoeff(&pos); 00915 if(maxRemainingSingularValue == RealScalar(0)) 00916 { 00917 m_nonzeroSingularValues = i; 00918 break; 00919 } 00920 if(pos) 00921 { 00922 pos += i; 00923 std::swap(m_singularValues.coeffRef(i), m_singularValues.coeffRef(pos)); 00924 if(computeU()) m_matrixU.col(pos).swap(m_matrixU.col(i)); 00925 if(computeV()) m_matrixV.col(pos).swap(m_matrixV.col(i)); 00926 } 00927 } 00928 00929 m_singularValues *= scale; 00930 00931 m_isInitialized = true; 00932 return *this; 00933 } 00934 00935 namespace internal { 00936 template<typename _MatrixType, int QRPreconditioner, typename Rhs> 00937 struct solve_retval<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs> 00938 : solve_retval_base<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs> 00939 { 00940 typedef JacobiSVD<_MatrixType, QRPreconditioner> JacobiSVDType; 00941 EIGEN_MAKE_SOLVE_HELPERS(JacobiSVDType,Rhs) 00942 00943 template<typename Dest> void evalTo(Dest& dst) const 00944 { 00945 eigen_assert(rhs().rows() == dec().rows()); 00946 00947 // A = U S V^* 00948 // So A^{-1} = V S^{-1} U^* 00949 00950 Matrix<Scalar, Dynamic, Rhs::ColsAtCompileTime, 0, _MatrixType::MaxRowsAtCompileTime, Rhs::MaxColsAtCompileTime> tmp; 00951 Index rank = dec().rank(); 00952 00953 tmp.noalias() = dec().matrixU().leftCols(rank).adjoint() * rhs(); 00954 tmp = dec().singularValues().head(rank).asDiagonal().inverse() * tmp; 00955 dst = dec().matrixV().leftCols(rank) * tmp; 00956 } 00957 }; 00958 } // end namespace internal 00959 00960 /** \svd_module 00961 * 00962 * \return the singular value decomposition of \c *this computed by two-sided 00963 * Jacobi transformations. 00964 * 00965 * \sa class JacobiSVD 00966 */ 00967 template<typename Derived> 00968 JacobiSVD<typename MatrixBase<Derived>::PlainObject> 00969 MatrixBase<Derived>::jacobiSvd (unsigned int computationOptions) const 00970 { 00971 return JacobiSVD<PlainObject>(*this, computationOptions); 00972 } 00973 00974 } // end namespace Eigen 00975 00976 #endif // EIGEN_JACOBISVD_H
Generated on Tue Jul 12 2022 17:46:56 by 1.7.2