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

source code of qtmultimedia/src/3rdparty/eigen/Eigen/src/SVD/JacobiSVD.h