Lines Matching refs:SVD
22 /*** QR preconditioners (R-SVD)
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
77 void allocate(const JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd)
79 if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols())
82 ::new (&m_qr) QRType(svd.rows(), svd.cols());
84 if (svd.m_computeFullU) m_workspace.resize(svd.rows());
87 bool run(JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd, const MatrixType& 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();
122 void allocate(const JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd)
124 if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols())
127 ::new (&m_qr) QRType(svd.cols(), svd.rows());
129 m_adjoint.resize(svd.cols(), svd.rows());
130 if (svd.m_computeFullV) m_workspace.resize(svd.cols());
133 bool run(JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
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();
161 void allocate(const JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd)
163 if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols())
166 ::new (&m_qr) QRType(svd.rows(), svd.cols());
168 if (svd.m_computeFullU) m_workspace.resize(svd.rows());
169 else if (svd.m_computeThinU) m_workspace.resize(svd.cols());
172 bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& 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)
181 svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols());
182 m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_workspace);
184 if(svd.computeV()) svd.m_matrixV = m_qr.colsPermutation();
214 void allocate(const JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd)
216 if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols())
219 ::new (&m_qr) QRType(svd.cols(), svd.rows());
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());
226 bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
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)
237 svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows());
238 m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_workspace);
240 if(svd.computeU()) svd.m_matrixU = m_qr.colsPermutation();
261 void allocate(const JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd)
263 if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols())
266 ::new (&m_qr) QRType(svd.rows(), svd.cols());
268 if (svd.m_computeFullU) m_workspace.resize(svd.rows());
269 else if (svd.m_computeThinU) m_workspace.resize(svd.cols());
272 bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& 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)
281 svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols());
282 m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_workspace);
284 if(svd.computeV()) svd.m_matrixV.setIdentity(matrix.cols(), matrix.cols());
313 void allocate(const JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd)
315 if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols())
318 ::new (&m_qr) QRType(svd.cols(), svd.rows());
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());
325 bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix)
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)
336 svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows());
337 m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_workspace);
339 if(svd.computeU()) svd.m_matrixU.setIdentity(matrix.rows(), matrix.rows());
352 /*** 2x2 SVD implementation
354 *** JacobiSVD consists in performing a series of 2x2 SVD subproblems
360 typedef JacobiSVD<MatrixType, QRPreconditioner> SVD;
361 typedef typename SVD::Index Index;
362 static void run(typename SVD::WorkMatrixType&, SVD&, Index, Index) {}
368 typedef JacobiSVD<MatrixType, QRPreconditioner> SVD;
371 typedef typename SVD::Index Index;
372 static void run(typename SVD::WorkMatrixType& work_matrix, SVD& svd, Index p, Index q)
382 if(svd.computeU()) svd.m_matrixU.col(p) *= conj(z);
385 if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z);
392 if(svd.computeU()) svd.m_matrixU.applyOnTheRight(p,q,rot.adjoint());
397 if(svd.computeV()) svd.m_matrixV.col(q) *= z;
403 if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z);
444 * \brief Two-sided Jacobi SVD decomposition of a rectangular matrix
446 * \param MatrixType the type of the matrix of which we are computing the SVD decomposition
448 * for the R-SVD step for non-square matrices. See discussion of possible values below.
450 * SVD decomposition consists in decomposing any n-by-p matrix \a A as a product
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.
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.
585 /** \returns a (least squares) solution of \f$ A x = b \f$ using the current SVD decomposition of A.
591 * \note SVD solving is implicitly least-squares. Thus, this method serves both purposes of exact solving and least-squares solving.
652 /*** step 1. The R-SVD step: we use a QR decomposition to reduce to the case of a square matrix */
663 /*** step 2. The main Jacobi SVD iteration. ***/
670 // do a sweep: for all index pairs (p,q), perform SVD of the corresponding 2x2 sub-matrix
686 // perform SVD decomposition of 2x2 sub-matrix corresponding to indices p,q to make it diagonal
