1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2009-2010 Benoit Jacob <jacob.benoit.1@gmail.com> 5 // 6 // This Source Code Form is subject to the terms of the Mozilla 7 // Public License v. 2.0. If a copy of the MPL was not distributed 8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9 10 #ifndef EIGEN_JACOBISVD_H 11 #define EIGEN_JACOBISVD_H 12 13 namespace Eigen { 14 15 namespace internal { 16 // forward declaration (needed by ICC) 17 // the empty body is required by MSVC 18 template<typename MatrixType, int QRPreconditioner, 19 bool IsComplex = NumTraits<typename MatrixType::Scalar>::IsComplex> 20 struct svd_precondition_2x2_block_to_be_real {}; 21 22 /*** QR preconditioners (R-SVD) 23 *** 24 *** Their role is to reduce the problem of computing the SVD to the case of a square matrix. 25 *** This approach, known as R-SVD, is an optimization for rectangular-enough matrices, and is a requirement for 26 *** JacobiSVD which by itself is only able to work on square matrices. 27 ***/ 28 29 enum { PreconditionIfMoreColsThanRows, PreconditionIfMoreRowsThanCols }; 30 31 template<typename MatrixType, int QRPreconditioner, int Case> 32 struct qr_preconditioner_should_do_anything 33 { 34 enum { a = MatrixType::RowsAtCompileTime != Dynamic && 35 MatrixType::ColsAtCompileTime != Dynamic && 36 MatrixType::ColsAtCompileTime <= MatrixType::RowsAtCompileTime, 37 b = MatrixType::RowsAtCompileTime != Dynamic && 38 MatrixType::ColsAtCompileTime != Dynamic && 39 MatrixType::RowsAtCompileTime <= MatrixType::ColsAtCompileTime, 40 ret = !( (QRPreconditioner == NoQRPreconditioner) || 41 (Case == PreconditionIfMoreColsThanRows && bool(a)) || 42 (Case == PreconditionIfMoreRowsThanCols && bool(b)) ) 43 }; 44 }; 45 46 template<typename MatrixType, int QRPreconditioner, int Case, 47 bool DoAnything = qr_preconditioner_should_do_anything<MatrixType, QRPreconditioner, Case>::ret 48 > struct qr_preconditioner_impl {}; 49 50 template<typename MatrixType, int QRPreconditioner, int Case> 51 class qr_preconditioner_impl<MatrixType, QRPreconditioner, Case, false> 52 { 53 public: 54 typedef typename MatrixType::Index Index; allocate(const JacobiSVD<MatrixType,QRPreconditioner> &)55 void allocate(const JacobiSVD<MatrixType, QRPreconditioner>&) {} run(JacobiSVD<MatrixType,QRPreconditioner> &,const MatrixType &)56 bool run(JacobiSVD<MatrixType, QRPreconditioner>&, const MatrixType&) 57 { 58 return false; 59 } 60 }; 61 62 /*** preconditioner using FullPivHouseholderQR ***/ 63 64 template<typename MatrixType> 65 class qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true> 66 { 67 public: 68 typedef typename MatrixType::Index Index; 69 typedef typename MatrixType::Scalar Scalar; 70 enum 71 { 72 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 73 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime 74 }; 75 typedef Matrix<Scalar, 1, RowsAtCompileTime, RowMajor, 1, MaxRowsAtCompileTime> WorkspaceType; 76 allocate(const JacobiSVD<MatrixType,FullPivHouseholderQRPreconditioner> & svd)77 void allocate(const JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd) 78 { 79 if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols()) 80 { 81 m_qr.~QRType(); 82 ::new (&m_qr) QRType(svd.rows(), svd.cols()); 83 } 84 if (svd.m_computeFullU) m_workspace.resize(svd.rows()); 85 } 86 run(JacobiSVD<MatrixType,FullPivHouseholderQRPreconditioner> & svd,const MatrixType & matrix)87 bool run(JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix) 88 { 89 if(matrix.rows() > matrix.cols()) 90 { 91 m_qr.compute(matrix); 92 svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>(); 93 if(svd.m_computeFullU) m_qr.matrixQ().evalTo(svd.m_matrixU, m_workspace); 94 if(svd.computeV()) svd.m_matrixV = m_qr.colsPermutation(); 95 return true; 96 } 97 return false; 98 } 99 private: 100 typedef FullPivHouseholderQR<MatrixType> QRType; 101 QRType m_qr; 102 WorkspaceType m_workspace; 103 }; 104 105 template<typename MatrixType> 106 class qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true> 107 { 108 public: 109 typedef typename MatrixType::Index Index; 110 typedef typename MatrixType::Scalar Scalar; 111 enum 112 { 113 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 114 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 115 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 116 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, 117 Options = MatrixType::Options 118 }; 119 typedef Matrix<Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime> 120 TransposeTypeWithSameStorageOrder; 121 allocate(const JacobiSVD<MatrixType,FullPivHouseholderQRPreconditioner> & svd)122 void allocate(const JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd) 123 { 124 if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols()) 125 { 126 m_qr.~QRType(); 127 ::new (&m_qr) QRType(svd.cols(), svd.rows()); 128 } 129 m_adjoint.resize(svd.cols(), svd.rows()); 130 if (svd.m_computeFullV) m_workspace.resize(svd.cols()); 131 } 132 run(JacobiSVD<MatrixType,FullPivHouseholderQRPreconditioner> & svd,const MatrixType & matrix)133 bool run(JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix) 134 { 135 if(matrix.cols() > matrix.rows()) 136 { 137 m_adjoint = matrix.adjoint(); 138 m_qr.compute(m_adjoint); 139 svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint(); 140 if(svd.m_computeFullV) m_qr.matrixQ().evalTo(svd.m_matrixV, m_workspace); 141 if(svd.computeU()) svd.m_matrixU = m_qr.colsPermutation(); 142 return true; 143 } 144 else return false; 145 } 146 private: 147 typedef FullPivHouseholderQR<TransposeTypeWithSameStorageOrder> QRType; 148 QRType m_qr; 149 TransposeTypeWithSameStorageOrder m_adjoint; 150 typename internal::plain_row_type<MatrixType>::type m_workspace; 151 }; 152 153 /*** preconditioner using ColPivHouseholderQR ***/ 154 155 template<typename MatrixType> 156 class qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true> 157 { 158 public: 159 typedef typename MatrixType::Index Index; 160 allocate(const JacobiSVD<MatrixType,ColPivHouseholderQRPreconditioner> & svd)161 void allocate(const JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd) 162 { 163 if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols()) 164 { 165 m_qr.~QRType(); 166 ::new (&m_qr) QRType(svd.rows(), svd.cols()); 167 } 168 if (svd.m_computeFullU) m_workspace.resize(svd.rows()); 169 else if (svd.m_computeThinU) m_workspace.resize(svd.cols()); 170 } 171 run(JacobiSVD<MatrixType,ColPivHouseholderQRPreconditioner> & svd,const MatrixType & matrix)172 bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix) 173 { 174 if(matrix.rows() > matrix.cols()) 175 { 176 m_qr.compute(matrix); 177 svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>(); 178 if(svd.m_computeFullU) m_qr.householderQ().evalTo(svd.m_matrixU, m_workspace); 179 else if(svd.m_computeThinU) 180 { 181 svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols()); 182 m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_workspace); 183 } 184 if(svd.computeV()) svd.m_matrixV = m_qr.colsPermutation(); 185 return true; 186 } 187 return false; 188 } 189 190 private: 191 typedef ColPivHouseholderQR<MatrixType> QRType; 192 QRType m_qr; 193 typename internal::plain_col_type<MatrixType>::type m_workspace; 194 }; 195 196 template<typename MatrixType> 197 class qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true> 198 { 199 public: 200 typedef typename MatrixType::Index Index; 201 typedef typename MatrixType::Scalar Scalar; 202 enum 203 { 204 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 205 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 206 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 207 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, 208 Options = MatrixType::Options 209 }; 210 211 typedef Matrix<Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime> 212 TransposeTypeWithSameStorageOrder; 213 allocate(const JacobiSVD<MatrixType,ColPivHouseholderQRPreconditioner> & svd)214 void allocate(const JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd) 215 { 216 if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols()) 217 { 218 m_qr.~QRType(); 219 ::new (&m_qr) QRType(svd.cols(), svd.rows()); 220 } 221 if (svd.m_computeFullV) m_workspace.resize(svd.cols()); 222 else if (svd.m_computeThinV) m_workspace.resize(svd.rows()); 223 m_adjoint.resize(svd.cols(), svd.rows()); 224 } 225 run(JacobiSVD<MatrixType,ColPivHouseholderQRPreconditioner> & svd,const MatrixType & matrix)226 bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix) 227 { 228 if(matrix.cols() > matrix.rows()) 229 { 230 m_adjoint = matrix.adjoint(); 231 m_qr.compute(m_adjoint); 232 233 svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint(); 234 if(svd.m_computeFullV) m_qr.householderQ().evalTo(svd.m_matrixV, m_workspace); 235 else if(svd.m_computeThinV) 236 { 237 svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows()); 238 m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_workspace); 239 } 240 if(svd.computeU()) svd.m_matrixU = m_qr.colsPermutation(); 241 return true; 242 } 243 else return false; 244 } 245 246 private: 247 typedef ColPivHouseholderQR<TransposeTypeWithSameStorageOrder> QRType; 248 QRType m_qr; 249 TransposeTypeWithSameStorageOrder m_adjoint; 250 typename internal::plain_row_type<MatrixType>::type m_workspace; 251 }; 252 253 /*** preconditioner using HouseholderQR ***/ 254 255 template<typename MatrixType> 256 class qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true> 257 { 258 public: 259 typedef typename MatrixType::Index Index; 260 allocate(const JacobiSVD<MatrixType,HouseholderQRPreconditioner> & svd)261 void allocate(const JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd) 262 { 263 if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols()) 264 { 265 m_qr.~QRType(); 266 ::new (&m_qr) QRType(svd.rows(), svd.cols()); 267 } 268 if (svd.m_computeFullU) m_workspace.resize(svd.rows()); 269 else if (svd.m_computeThinU) m_workspace.resize(svd.cols()); 270 } 271 run(JacobiSVD<MatrixType,HouseholderQRPreconditioner> & svd,const MatrixType & matrix)272 bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix) 273 { 274 if(matrix.rows() > matrix.cols()) 275 { 276 m_qr.compute(matrix); 277 svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>(); 278 if(svd.m_computeFullU) m_qr.householderQ().evalTo(svd.m_matrixU, m_workspace); 279 else if(svd.m_computeThinU) 280 { 281 svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols()); 282 m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_workspace); 283 } 284 if(svd.computeV()) svd.m_matrixV.setIdentity(matrix.cols(), matrix.cols()); 285 return true; 286 } 287 return false; 288 } 289 private: 290 typedef HouseholderQR<MatrixType> QRType; 291 QRType m_qr; 292 typename internal::plain_col_type<MatrixType>::type m_workspace; 293 }; 294 295 template<typename MatrixType> 296 class qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true> 297 { 298 public: 299 typedef typename MatrixType::Index Index; 300 typedef typename MatrixType::Scalar Scalar; 301 enum 302 { 303 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 304 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 305 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 306 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, 307 Options = MatrixType::Options 308 }; 309 310 typedef Matrix<Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime> 311 TransposeTypeWithSameStorageOrder; 312 allocate(const JacobiSVD<MatrixType,HouseholderQRPreconditioner> & svd)313 void allocate(const JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd) 314 { 315 if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols()) 316 { 317 m_qr.~QRType(); 318 ::new (&m_qr) QRType(svd.cols(), svd.rows()); 319 } 320 if (svd.m_computeFullV) m_workspace.resize(svd.cols()); 321 else if (svd.m_computeThinV) m_workspace.resize(svd.rows()); 322 m_adjoint.resize(svd.cols(), svd.rows()); 323 } 324 run(JacobiSVD<MatrixType,HouseholderQRPreconditioner> & svd,const MatrixType & matrix)325 bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix) 326 { 327 if(matrix.cols() > matrix.rows()) 328 { 329 m_adjoint = matrix.adjoint(); 330 m_qr.compute(m_adjoint); 331 332 svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint(); 333 if(svd.m_computeFullV) m_qr.householderQ().evalTo(svd.m_matrixV, m_workspace); 334 else if(svd.m_computeThinV) 335 { 336 svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows()); 337 m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_workspace); 338 } 339 if(svd.computeU()) svd.m_matrixU.setIdentity(matrix.rows(), matrix.rows()); 340 return true; 341 } 342 else return false; 343 } 344 345 private: 346 typedef HouseholderQR<TransposeTypeWithSameStorageOrder> QRType; 347 QRType m_qr; 348 TransposeTypeWithSameStorageOrder m_adjoint; 349 typename internal::plain_row_type<MatrixType>::type m_workspace; 350 }; 351 352 /*** 2x2 SVD implementation 353 *** 354 *** JacobiSVD consists in performing a series of 2x2 SVD subproblems 355 ***/ 356 357 template<typename MatrixType, int QRPreconditioner> 358 struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, false> 359 { 360 typedef JacobiSVD<MatrixType, QRPreconditioner> SVD; 361 typedef typename SVD::Index Index; 362 static void run(typename SVD::WorkMatrixType&, SVD&, Index, Index) {} 363 }; 364 365 template<typename MatrixType, int QRPreconditioner> 366 struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, true> 367 { 368 typedef JacobiSVD<MatrixType, QRPreconditioner> SVD; 369 typedef typename MatrixType::Scalar Scalar; 370 typedef typename MatrixType::RealScalar RealScalar; 371 typedef typename SVD::Index Index; 372 static void run(typename SVD::WorkMatrixType& work_matrix, SVD& svd, Index p, Index q) 373 { 374 using std::sqrt; 375 Scalar z; 376 JacobiRotation<Scalar> rot; 377 RealScalar n = sqrt(numext::abs2(work_matrix.coeff(p,p)) + numext::abs2(work_matrix.coeff(q,p))); 378 if(n==0) 379 { 380 z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q); 381 work_matrix.row(p) *= z; 382 if(svd.computeU()) svd.m_matrixU.col(p) *= conj(z); 383 z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q); 384 work_matrix.row(q) *= z; 385 if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z); 386 } 387 else 388 { 389 rot.c() = conj(work_matrix.coeff(p,p)) / n; 390 rot.s() = work_matrix.coeff(q,p) / n; 391 work_matrix.applyOnTheLeft(p,q,rot); 392 if(svd.computeU()) svd.m_matrixU.applyOnTheRight(p,q,rot.adjoint()); 393 if(work_matrix.coeff(p,q) != Scalar(0)) 394 { 395 Scalar z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q); 396 work_matrix.col(q) *= z; 397 if(svd.computeV()) svd.m_matrixV.col(q) *= z; 398 } 399 if(work_matrix.coeff(q,q) != Scalar(0)) 400 { 401 z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q); 402 work_matrix.row(q) *= z; 403 if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z); 404 } 405 } 406 } 407 }; 408 409 template<typename MatrixType, typename RealScalar, typename Index> 410 void real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q, 411 JacobiRotation<RealScalar> *j_left, 412 JacobiRotation<RealScalar> *j_right) 413 { 414 using std::sqrt; 415 Matrix<RealScalar,2,2> m; 416 m << numext::real(matrix.coeff(p,p)), numext::real(matrix.coeff(p,q)), 417 numext::real(matrix.coeff(q,p)), numext::real(matrix.coeff(q,q)); 418 JacobiRotation<RealScalar> rot1; 419 RealScalar t = m.coeff(0,0) + m.coeff(1,1); 420 RealScalar d = m.coeff(1,0) - m.coeff(0,1); 421 if(t == RealScalar(0)) 422 { 423 rot1.c() = RealScalar(0); 424 rot1.s() = d > RealScalar(0) ? RealScalar(1) : RealScalar(-1); 425 } 426 else 427 { 428 RealScalar u = d / t; 429 rot1.c() = RealScalar(1) / sqrt(RealScalar(1) + numext::abs2(u)); 430 rot1.s() = rot1.c() * u; 431 } 432 m.applyOnTheLeft(0,1,rot1); 433 j_right->makeJacobi(m,0,1); 434 *j_left = rot1 * j_right->transpose(); 435 } 436 437 } // end namespace internal 438 439 /** \ingroup SVD_Module 440 * 441 * 442 * \class JacobiSVD 443 * 444 * \brief Two-sided Jacobi SVD decomposition of a rectangular matrix 445 * 446 * \param MatrixType the type of the matrix of which we are computing the SVD decomposition 447 * \param QRPreconditioner this optional parameter allows to specify the type of QR decomposition that will be used internally 448 * for the R-SVD step for non-square matrices. See discussion of possible values below. 449 * 450 * SVD decomposition consists in decomposing any n-by-p matrix \a A as a product 451 * \f[ A = U S V^* \f] 452 * 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; 453 * 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 454 * and right \em singular \em vectors of \a A respectively. 455 * 456 * Singular values are always sorted in decreasing order. 457 * 458 * 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. 459 * 460 * 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 461 * 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 462 * 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, 463 * 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. 464 * 465 * Here's an example demonstrating basic usage: 466 * \include JacobiSVD_basic.cpp 467 * Output: \verbinclude JacobiSVD_basic.out 468 * 469 * This JacobiSVD class is a two-sided Jacobi R-SVD decomposition, ensuring optimal reliability and accuracy. The downside is that it's slower than 470 * 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 471 * \a p is the greater dimension, meaning that it is still of the same order of complexity as the faster bidiagonalizing R-SVD algorithms. 472 * In particular, like any R-SVD, it takes advantage of non-squareness in that its complexity is only linear in the greater dimension. 473 * 474 * If the input matrix has inf or nan coefficients, the result of the computation is undefined, but the computation is guaranteed to 475 * terminate in finite (and reasonable) time. 476 * 477 * The possible values for QRPreconditioner are: 478 * \li ColPivHouseholderQRPreconditioner is the default. In practice it's very safe. It uses column-pivoting QR. 479 * \li FullPivHouseholderQRPreconditioner, is the safest and slowest. It uses full-pivoting QR. 480 * Contrary to other QRs, it doesn't allow computing thin unitaries. 481 * \li HouseholderQRPreconditioner is the fastest, and less safe and accurate than the pivoting variants. It uses non-pivoting QR. 482 * This is very similar in safety and accuracy to the bidiagonalization process used by bidiagonalizing SVD algorithms (since bidiagonalization 483 * is inherently non-pivoting). However the resulting SVD is still more reliable than bidiagonalizing SVDs because the Jacobi-based iterarive 484 * process is more reliable than the optimized bidiagonal SVD iterations. 485 * \li NoQRPreconditioner allows not to use a QR preconditioner at all. This is useful if you know that you will only be computing 486 * JacobiSVD decompositions of square matrices. Non-square matrices require a QR preconditioner. Using this option will result in 487 * faster compilation and smaller executable code. It won't significantly speed up computation, since JacobiSVD is always checking 488 * if QR preconditioning is needed before applying it anyway. 489 * 490 * \sa MatrixBase::jacobiSvd() 491 */ 492 template<typename _MatrixType, int QRPreconditioner> 493 class JacobiSVD : public SVDBase<_MatrixType> 494 { 495 public: 496 497 typedef _MatrixType MatrixType; 498 typedef typename MatrixType::Scalar Scalar; 499 typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; 500 typedef typename MatrixType::Index Index; 501 enum { 502 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 503 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 504 DiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_DYNAMIC(RowsAtCompileTime,ColsAtCompileTime), 505 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 506 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, 507 MaxDiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(MaxRowsAtCompileTime,MaxColsAtCompileTime), 508 MatrixOptions = MatrixType::Options 509 }; 510 511 typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, 512 MatrixOptions, MaxRowsAtCompileTime, MaxRowsAtCompileTime> 513 MatrixUType; 514 typedef Matrix<Scalar, ColsAtCompileTime, ColsAtCompileTime, 515 MatrixOptions, MaxColsAtCompileTime, MaxColsAtCompileTime> 516 MatrixVType; 517 typedef typename internal::plain_diag_type<MatrixType, RealScalar>::type SingularValuesType; 518 typedef typename internal::plain_row_type<MatrixType>::type RowType; 519 typedef typename internal::plain_col_type<MatrixType>::type ColType; 520 typedef Matrix<Scalar, DiagSizeAtCompileTime, DiagSizeAtCompileTime, 521 MatrixOptions, MaxDiagSizeAtCompileTime, MaxDiagSizeAtCompileTime> 522 WorkMatrixType; 523 524 /** \brief Default Constructor. 525 * 526 * The default constructor is useful in cases in which the user intends to 527 * perform decompositions via JacobiSVD::compute(const MatrixType&). 528 */ 529 JacobiSVD() 530 : SVDBase<_MatrixType>::SVDBase() 531 {} 532 533 534 /** \brief Default Constructor with memory preallocation 535 * 536 * Like the default constructor but with preallocation of the internal data 537 * according to the specified problem size. 538 * \sa JacobiSVD() 539 */ 540 JacobiSVD(Index rows, Index cols, unsigned int computationOptions = 0) 541 : SVDBase<_MatrixType>::SVDBase() 542 { 543 allocate(rows, cols, computationOptions); 544 } 545 546 /** \brief Constructor performing the decomposition of given matrix. 547 * 548 * \param matrix the matrix to decompose 549 * \param computationOptions optional parameter allowing to specify if you want full or thin U or V unitaries to be computed. 550 * By default, none is computed. This is a bit-field, the possible bits are #ComputeFullU, #ComputeThinU, 551 * #ComputeFullV, #ComputeThinV. 552 * 553 * Thin unitaries are only available if your matrix type has a Dynamic number of columns (for example MatrixXf). They also are not 554 * available with the (non-default) FullPivHouseholderQR preconditioner. 555 */ 556 JacobiSVD(const MatrixType& matrix, unsigned int computationOptions = 0) 557 : SVDBase<_MatrixType>::SVDBase() 558 { 559 compute(matrix, computationOptions); 560 } 561 562 /** \brief Method performing the decomposition of given matrix using custom options. 563 * 564 * \param matrix the matrix to decompose 565 * \param computationOptions optional parameter allowing to specify if you want full or thin U or V unitaries to be computed. 566 * By default, none is computed. This is a bit-field, the possible bits are #ComputeFullU, #ComputeThinU, 567 * #ComputeFullV, #ComputeThinV. 568 * 569 * Thin unitaries are only available if your matrix type has a Dynamic number of columns (for example MatrixXf). They also are not 570 * available with the (non-default) FullPivHouseholderQR preconditioner. 571 */ 572 SVDBase<MatrixType>& compute(const MatrixType& matrix, unsigned int computationOptions); 573 574 /** \brief Method performing the decomposition of given matrix using current options. 575 * 576 * \param matrix the matrix to decompose 577 * 578 * This method uses the current \a computationOptions, as already passed to the constructor or to compute(const MatrixType&, unsigned int). 579 */ 580 SVDBase<MatrixType>& compute(const MatrixType& matrix) 581 { 582 return compute(matrix, this->m_computationOptions); 583 } 584 585 /** \returns a (least squares) solution of \f$ A x = b \f$ using the current SVD decomposition of A. 586 * 587 * \param b the right-hand-side of the equation to solve. 588 * 589 * \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. 590 * 591 * \note SVD solving is implicitly least-squares. Thus, this method serves both purposes of exact solving and least-squares solving. 592 * In other words, the returned solution is guaranteed to minimize the Euclidean norm \f$ \Vert A x - b \Vert \f$. 593 */ 594 template<typename Rhs> 595 inline const internal::solve_retval<JacobiSVD, Rhs> 596 solve(const MatrixBase<Rhs>& b) const 597 { 598 eigen_assert(this->m_isInitialized && "JacobiSVD is not initialized."); 599 eigen_assert(SVDBase<MatrixType>::computeU() && SVDBase<MatrixType>::computeV() && "JacobiSVD::solve() requires both unitaries U and V to be computed (thin unitaries suffice)."); 600 return internal::solve_retval<JacobiSVD, Rhs>(*this, b.derived()); 601 } 602 603 604 605 private: 606 void allocate(Index rows, Index cols, unsigned int computationOptions); 607 608 protected: 609 WorkMatrixType m_workMatrix; 610 611 template<typename __MatrixType, int _QRPreconditioner, bool _IsComplex> 612 friend struct internal::svd_precondition_2x2_block_to_be_real; 613 template<typename __MatrixType, int _QRPreconditioner, int _Case, bool _DoAnything> 614 friend struct internal::qr_preconditioner_impl; 615 616 internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreColsThanRows> m_qr_precond_morecols; 617 internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreRowsThanCols> m_qr_precond_morerows; 618 }; 619 620 template<typename MatrixType, int QRPreconditioner> 621 void JacobiSVD<MatrixType, QRPreconditioner>::allocate(Index rows, Index cols, unsigned int computationOptions) 622 { 623 if (SVDBase<MatrixType>::allocate(rows, cols, computationOptions)) return; 624 625 if (QRPreconditioner == FullPivHouseholderQRPreconditioner) 626 { 627 eigen_assert(!(this->m_computeThinU || this->m_computeThinV) && 628 "JacobiSVD: can't compute thin U or thin V with the FullPivHouseholderQR preconditioner. " 629 "Use the ColPivHouseholderQR preconditioner instead."); 630 } 631 632 m_workMatrix.resize(this->m_diagSize, this->m_diagSize); 633 634 if(this->m_cols>this->m_rows) m_qr_precond_morecols.allocate(*this); 635 if(this->m_rows>this->m_cols) m_qr_precond_morerows.allocate(*this); 636 } 637 638 template<typename MatrixType, int QRPreconditioner> 639 SVDBase<MatrixType>& 640 JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsigned int computationOptions) 641 { 642 using std::abs; 643 allocate(matrix.rows(), matrix.cols(), computationOptions); 644 645 // currently we stop when we reach precision 2*epsilon as the last bit of precision can require an unreasonable number of iterations, 646 // only worsening the precision of U and V as we accumulate more rotations 647 const RealScalar precision = RealScalar(2) * NumTraits<Scalar>::epsilon(); 648 649 // limit for very small denormal numbers to be considered zero in order to avoid infinite loops (see bug 286) 650 const RealScalar considerAsZero = RealScalar(2) * std::numeric_limits<RealScalar>::denorm_min(); 651 652 /*** step 1. The R-SVD step: we use a QR decomposition to reduce to the case of a square matrix */ 653 654 if(!m_qr_precond_morecols.run(*this, matrix) && !m_qr_precond_morerows.run(*this, matrix)) 655 { 656 m_workMatrix = matrix.block(0,0,this->m_diagSize,this->m_diagSize); 657 if(this->m_computeFullU) this->m_matrixU.setIdentity(this->m_rows,this->m_rows); 658 if(this->m_computeThinU) this->m_matrixU.setIdentity(this->m_rows,this->m_diagSize); 659 if(this->m_computeFullV) this->m_matrixV.setIdentity(this->m_cols,this->m_cols); 660 if(this->m_computeThinV) this->m_matrixV.setIdentity(this->m_cols, this->m_diagSize); 661 } 662 663 /*** step 2. The main Jacobi SVD iteration. ***/ 664 665 bool finished = false; 666 while(!finished) 667 { 668 finished = true; 669 670 // do a sweep: for all index pairs (p,q), perform SVD of the corresponding 2x2 sub-matrix 671 672 for(Index p = 1; p < this->m_diagSize; ++p) 673 { 674 for(Index q = 0; q < p; ++q) 675 { 676 // if this 2x2 sub-matrix is not diagonal already... 677 // notice that this comparison will evaluate to false if any NaN is involved, ensuring that NaN's don't 678 // keep us iterating forever. Similarly, small denormal numbers are considered zero. 679 using std::max; 680 RealScalar threshold = (max)(considerAsZero, precision * (max)(abs(m_workMatrix.coeff(p,p)), 681 abs(m_workMatrix.coeff(q,q)))); 682 if((max)(abs(m_workMatrix.coeff(p,q)),abs(m_workMatrix.coeff(q,p))) > threshold) 683 { 684 finished = false; 685 686 // perform SVD decomposition of 2x2 sub-matrix corresponding to indices p,q to make it diagonal 687 internal::svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner>::run(m_workMatrix, *this, p, q); 688 JacobiRotation<RealScalar> j_left, j_right; 689 internal::real_2x2_jacobi_svd(m_workMatrix, p, q, &j_left, &j_right); 690 691 // accumulate resulting Jacobi rotations 692 m_workMatrix.applyOnTheLeft(p,q,j_left); 693 if(SVDBase<MatrixType>::computeU()) this->m_matrixU.applyOnTheRight(p,q,j_left.transpose()); 694 695 m_workMatrix.applyOnTheRight(p,q,j_right); 696 if(SVDBase<MatrixType>::computeV()) this->m_matrixV.applyOnTheRight(p,q,j_right); 697 } 698 } 699 } 700 } 701 702 /*** step 3. The work matrix is now diagonal, so ensure it's positive so its diagonal entries are the singular values ***/ 703 704 for(Index i = 0; i < this->m_diagSize; ++i) 705 { 706 RealScalar a = abs(m_workMatrix.coeff(i,i)); 707 this->m_singularValues.coeffRef(i) = a; 708 if(SVDBase<MatrixType>::computeU() && (a!=RealScalar(0))) this->m_matrixU.col(i) *= this->m_workMatrix.coeff(i,i)/a; 709 } 710 711 /*** step 4. Sort singular values in descending order and compute the number of nonzero singular values ***/ 712 713 this->m_nonzeroSingularValues = this->m_diagSize; 714 for(Index i = 0; i < this->m_diagSize; i++) 715 { 716 Index pos; 717 RealScalar maxRemainingSingularValue = this->m_singularValues.tail(this->m_diagSize-i).maxCoeff(&pos); 718 if(maxRemainingSingularValue == RealScalar(0)) 719 { 720 this->m_nonzeroSingularValues = i; 721 break; 722 } 723 if(pos) 724 { 725 pos += i; 726 std::swap(this->m_singularValues.coeffRef(i), this->m_singularValues.coeffRef(pos)); 727 if(SVDBase<MatrixType>::computeU()) this->m_matrixU.col(pos).swap(this->m_matrixU.col(i)); 728 if(SVDBase<MatrixType>::computeV()) this->m_matrixV.col(pos).swap(this->m_matrixV.col(i)); 729 } 730 } 731 732 this->m_isInitialized = true; 733 return *this; 734 } 735 736 namespace internal { 737 template<typename _MatrixType, int QRPreconditioner, typename Rhs> 738 struct solve_retval<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs> 739 : solve_retval_base<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs> 740 { 741 typedef JacobiSVD<_MatrixType, QRPreconditioner> JacobiSVDType; 742 EIGEN_MAKE_SOLVE_HELPERS(JacobiSVDType,Rhs) 743 744 template<typename Dest> void evalTo(Dest& dst) const 745 { 746 eigen_assert(rhs().rows() == dec().rows()); 747 748 // A = U S V^* 749 // So A^{-1} = V S^{-1} U^* 750 751 Index diagSize = (std::min)(dec().rows(), dec().cols()); 752 typename JacobiSVDType::SingularValuesType invertedSingVals(diagSize); 753 754 Index nonzeroSingVals = dec().nonzeroSingularValues(); 755 invertedSingVals.head(nonzeroSingVals) = dec().singularValues().head(nonzeroSingVals).array().inverse(); 756 invertedSingVals.tail(diagSize - nonzeroSingVals).setZero(); 757 758 dst = dec().matrixV().leftCols(diagSize) 759 * invertedSingVals.asDiagonal() 760 * dec().matrixU().leftCols(diagSize).adjoint() 761 * rhs(); 762 } 763 }; 764 } // end namespace internal 765 766 /** \svd_module 767 * 768 * \return the singular value decomposition of \c *this computed by two-sided 769 * Jacobi transformations. 770 * 771 * \sa class JacobiSVD 772 */ 773 template<typename Derived> 774 JacobiSVD<typename MatrixBase<Derived>::PlainObject> 775 MatrixBase<Derived>::jacobiSvd(unsigned int computationOptions) const 776 { 777 return JacobiSVD<PlainObject>(*this, computationOptions); 778 } 779 780 } // end namespace Eigen 781 782 #endif // EIGEN_JACOBISVD_H 783