Eigen  3.4.90 (git rev a4098ac676528a83cfb73d4d26ce1b42ec05f47c)
JacobiSVD.h
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
15
16namespace Eigen {
17
18namespace internal {
19// forward declaration (needed by ICC)
20// the empty body is required by MSVC
21template<typename MatrixType, int QRPreconditioner,
22 bool IsComplex = NumTraits<typename MatrixType::Scalar>::IsComplex>
23struct svd_precondition_2x2_block_to_be_real {};
24
25/*** QR preconditioners (R-SVD)
26 ***
27 *** Their role is to reduce the problem of computing the SVD to the case of a square matrix.
28 *** This approach, known as R-SVD, is an optimization for rectangular-enough matrices, and is a requirement for
29 *** JacobiSVD which by itself is only able to work on square matrices.
30 ***/
31
32enum { PreconditionIfMoreColsThanRows, PreconditionIfMoreRowsThanCols };
33
34template<typename MatrixType, int QRPreconditioner, int Case>
35struct qr_preconditioner_should_do_anything
36{
37 enum { a = MatrixType::RowsAtCompileTime != Dynamic &&
38 MatrixType::ColsAtCompileTime != Dynamic &&
39 MatrixType::ColsAtCompileTime <= MatrixType::RowsAtCompileTime,
40 b = MatrixType::RowsAtCompileTime != Dynamic &&
41 MatrixType::ColsAtCompileTime != Dynamic &&
42 MatrixType::RowsAtCompileTime <= MatrixType::ColsAtCompileTime,
43 ret = !( (QRPreconditioner == NoQRPreconditioner) ||
44 (Case == PreconditionIfMoreColsThanRows && bool(a)) ||
45 (Case == PreconditionIfMoreRowsThanCols && bool(b)) )
46 };
47};
48
49template<typename MatrixType, int QRPreconditioner, int Case,
50 bool DoAnything = qr_preconditioner_should_do_anything<MatrixType, QRPreconditioner, Case>::ret
51> struct qr_preconditioner_impl {};
52
53template<typename MatrixType, int QRPreconditioner, int Case>
54class qr_preconditioner_impl<MatrixType, QRPreconditioner, Case, false>
55{
56public:
57 void allocate(const JacobiSVD<MatrixType, QRPreconditioner>&) {}
58 bool run(JacobiSVD<MatrixType, QRPreconditioner>&, const MatrixType&)
59 {
60 return false;
61 }
62};
63
64/*** preconditioner using FullPivHouseholderQR ***/
65
66template<typename MatrixType>
67class qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
68{
69public:
70 typedef typename MatrixType::Scalar Scalar;
71 enum
72 {
73 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
74 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime
75 };
76 typedef Matrix<Scalar, 1, RowsAtCompileTime, RowMajor, 1, MaxRowsAtCompileTime> WorkspaceType;
77
78 void allocate(const JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd)
79 {
80 if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols())
81 {
82 m_qr.~QRType();
83 ::new (&m_qr) QRType(svd.rows(), svd.cols());
84 }
85 if (svd.m_computeFullU) m_workspace.resize(svd.rows());
86 }
87
88 bool run(JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
89 {
90 if(matrix.rows() > matrix.cols())
91 {
92 m_qr.compute(matrix);
93 svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>();
94 if(svd.m_computeFullU) m_qr.matrixQ().evalTo(svd.m_matrixU, m_workspace);
95 if(svd.computeV()) svd.m_matrixV = m_qr.colsPermutation();
96 return true;
97 }
98 return false;
99 }
100private:
101 typedef FullPivHouseholderQR<MatrixType> QRType;
102 QRType m_qr;
103 WorkspaceType m_workspace;
104};
105
106template<typename MatrixType>
107class qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true>
108{
109public:
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
120 typedef typename internal::make_proper_matrix_type<
121 Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime
122 >::type TransposeTypeWithSameStorageOrder;
123
124 void allocate(const JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd)
125 {
126 if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols())
127 {
128 m_qr.~QRType();
129 ::new (&m_qr) QRType(svd.cols(), svd.rows());
130 }
132 if (svd.m_computeFullV) m_workspace.resize(svd.cols());
133 }
134
135 bool run(JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
136 {
137 if(matrix.cols() > matrix.rows())
138 {
142 if(svd.m_computeFullV) m_qr.matrixQ().evalTo(svd.m_matrixV, m_workspace);
143 if(svd.computeU()) svd.m_matrixU = m_qr.colsPermutation();
144 return true;
145 }
146 else return false;
147 }
148private:
149 typedef FullPivHouseholderQR<TransposeTypeWithSameStorageOrder> QRType;
150 QRType m_qr;
152 typename internal::plain_row_type<MatrixType>::type m_workspace;
153};
154
155/*** preconditioner using ColPivHouseholderQR ***/
156
157template<typename MatrixType>
158class qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
159{
160public:
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
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
190private:
191 typedef ColPivHouseholderQR<MatrixType> QRType;
192 QRType m_qr;
193 typename internal::plain_col_type<MatrixType>::type m_workspace;
194};
195
196template<typename MatrixType>
197class qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true>
198{
199public:
200 typedef typename MatrixType::Scalar Scalar;
201 enum
202 {
203 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
204 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
205 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
206 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
207 Options = MatrixType::Options
208 };
209
210 typedef typename internal::make_proper_matrix_type<
211 Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime
212 >::type TransposeTypeWithSameStorageOrder;
213
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());
224 }
225
226 bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
227 {
228 if(matrix.cols() > matrix.rows())
229 {
232
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
246private:
247 typedef ColPivHouseholderQR<TransposeTypeWithSameStorageOrder> QRType;
248 QRType m_qr;
250 typename internal::plain_row_type<MatrixType>::type m_workspace;
251};
252
253/*** preconditioner using HouseholderQR ***/
254
255template<typename MatrixType>
256class qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
257{
258public:
259 void allocate(const JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd)
260 {
261 if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols())
262 {
263 m_qr.~QRType();
264 ::new (&m_qr) QRType(svd.rows(), svd.cols());
265 }
266 if (svd.m_computeFullU) m_workspace.resize(svd.rows());
267 else if (svd.m_computeThinU) m_workspace.resize(svd.cols());
268 }
269
270 bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix)
271 {
272 if(matrix.rows() > matrix.cols())
273 {
274 m_qr.compute(matrix);
275 svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>();
276 if(svd.m_computeFullU) m_qr.householderQ().evalTo(svd.m_matrixU, m_workspace);
277 else if(svd.m_computeThinU)
278 {
279 svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols());
280 m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_workspace);
281 }
282 if(svd.computeV()) svd.m_matrixV.setIdentity(matrix.cols(), matrix.cols());
283 return true;
284 }
285 return false;
286 }
287private:
288 typedef HouseholderQR<MatrixType> QRType;
289 QRType m_qr;
290 typename internal::plain_col_type<MatrixType>::type m_workspace;
291};
292
293template<typename MatrixType>
294class qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true>
295{
296public:
297 typedef typename MatrixType::Scalar Scalar;
298 enum
299 {
300 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
301 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
302 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
303 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
304 Options = MatrixType::Options
305 };
306
307 typedef typename internal::make_proper_matrix_type<
308 Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime
309 >::type TransposeTypeWithSameStorageOrder;
310
311 void allocate(const JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd)
312 {
313 if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols())
314 {
315 m_qr.~QRType();
316 ::new (&m_qr) QRType(svd.cols(), svd.rows());
317 }
318 if (svd.m_computeFullV) m_workspace.resize(svd.cols());
319 else if (svd.m_computeThinV) m_workspace.resize(svd.rows());
321 }
322
323 bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix)
324 {
325 if(matrix.cols() > matrix.rows())
326 {
329
331 if(svd.m_computeFullV) m_qr.householderQ().evalTo(svd.m_matrixV, m_workspace);
332 else if(svd.m_computeThinV)
333 {
334 svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows());
335 m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_workspace);
336 }
337 if(svd.computeU()) svd.m_matrixU.setIdentity(matrix.rows(), matrix.rows());
338 return true;
339 }
340 else return false;
341 }
342
343private:
344 typedef HouseholderQR<TransposeTypeWithSameStorageOrder> QRType;
345 QRType m_qr;
347 typename internal::plain_row_type<MatrixType>::type m_workspace;
348};
349
350/*** 2x2 SVD implementation
351 ***
352 *** JacobiSVD consists in performing a series of 2x2 SVD subproblems
353 ***/
354
355template<typename MatrixType, int QRPreconditioner>
356struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, false>
357{
358 typedef JacobiSVD<MatrixType, QRPreconditioner> SVD;
359 typedef typename MatrixType::RealScalar RealScalar;
360 static bool run(typename SVD::WorkMatrixType&, SVD&, Index, Index, RealScalar&) { return true; }
361};
362
363template<typename MatrixType, int QRPreconditioner>
364struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, true>
365{
366 typedef JacobiSVD<MatrixType, QRPreconditioner> SVD;
367 typedef typename MatrixType::Scalar Scalar;
368 typedef typename MatrixType::RealScalar RealScalar;
369 static bool run(typename SVD::WorkMatrixType& work_matrix, SVD& svd, Index p, Index q, RealScalar& maxDiagEntry)
370 {
371 using std::sqrt;
372 using std::abs;
375 RealScalar n = sqrt(numext::abs2(work_matrix.coeff(p,p)) + numext::abs2(work_matrix.coeff(q,p)));
376
377 const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)();
378 const RealScalar precision = NumTraits<Scalar>::epsilon();
379
380 if(n==0)
381 {
382 // make sure first column is zero
383 work_matrix.coeffRef(p,p) = work_matrix.coeffRef(q,p) = Scalar(0);
384
385 if(abs(numext::imag(work_matrix.coeff(p,q)))>considerAsZero)
386 {
387 // 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
388 z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q);
389 work_matrix.row(p) *= z;
390 if(svd.computeU()) svd.m_matrixU.col(p) *= conj(z);
391 }
392 if(abs(numext::imag(work_matrix.coeff(q,q)))>considerAsZero)
393 {
394 z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q);
395 work_matrix.row(q) *= z;
396 if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z);
397 }
398 // otherwise the second row is already zero, so we have nothing to do.
399 }
400 else
401 {
402 rot.c() = conj(work_matrix.coeff(p,p)) / n;
403 rot.s() = work_matrix.coeff(q,p) / n;
404 work_matrix.applyOnTheLeft(p,q,rot);
406 if(abs(numext::imag(work_matrix.coeff(p,q)))>considerAsZero)
407 {
408 z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q);
409 work_matrix.col(q) *= z;
410 if(svd.computeV()) svd.m_matrixV.col(q) *= z;
411 }
412 if(abs(numext::imag(work_matrix.coeff(q,q)))>considerAsZero)
413 {
414 z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q);
415 work_matrix.row(q) *= z;
416 if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z);
417 }
418 }
419
420 // update largest diagonal entry
421 maxDiagEntry = numext::maxi<RealScalar>(maxDiagEntry,numext::maxi<RealScalar>(abs(work_matrix.coeff(p,p)), abs(work_matrix.coeff(q,q))));
422 // and check whether the 2x2 block is already diagonal
423 RealScalar threshold = numext::maxi<RealScalar>(considerAsZero, precision * maxDiagEntry);
424 return abs(work_matrix.coeff(p,q))>threshold || abs(work_matrix.coeff(q,p)) > threshold;
425 }
426};
427
428template<typename MatrixType_, int QRPreconditioner>
429struct traits<JacobiSVD<MatrixType_,QRPreconditioner> >
430 : traits<MatrixType_>
431{
432 typedef MatrixType_ MatrixType;
433};
434
435} // end namespace internal
436
490template<typename MatrixType_, int QRPreconditioner> class JacobiSVD
491 : public SVDBase<JacobiSVD<MatrixType_,QRPreconditioner> >
492{
493 typedef SVDBase<JacobiSVD> Base;
494 public:
495
496 typedef MatrixType_ MatrixType;
497 typedef typename MatrixType::Scalar Scalar;
498 typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
499 enum {
500 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
501 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
502 DiagSizeAtCompileTime = internal::min_size_prefer_dynamic(RowsAtCompileTime,ColsAtCompileTime),
503 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
504 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
505 MaxDiagSizeAtCompileTime = internal::min_size_prefer_fixed(MaxRowsAtCompileTime,MaxColsAtCompileTime),
506 MatrixOptions = MatrixType::Options
507 };
508
509 typedef typename Base::MatrixUType MatrixUType;
510 typedef typename Base::MatrixVType MatrixVType;
511 typedef typename Base::SingularValuesType SingularValuesType;
512
513 typedef typename internal::plain_row_type<MatrixType>::type RowType;
514 typedef typename internal::plain_col_type<MatrixType>::type ColType;
515 typedef Matrix<Scalar, DiagSizeAtCompileTime, DiagSizeAtCompileTime,
516 MatrixOptions, MaxDiagSizeAtCompileTime, MaxDiagSizeAtCompileTime>
518
525 {}
526
527
534 JacobiSVD(Index rows, Index cols, unsigned int computationOptions = 0)
535 {
536 allocate(rows, cols, computationOptions);
537 }
538
549 explicit JacobiSVD(const MatrixType& matrix, unsigned int computationOptions = 0)
550 {
551 compute(matrix, computationOptions);
552 }
553
564 JacobiSVD& compute(const MatrixType& matrix, unsigned int computationOptions);
565
572 JacobiSVD& compute(const MatrixType& matrix)
573 {
574 return compute(matrix, m_computationOptions);
575 }
576
577 using Base::computeU;
578 using Base::computeV;
579 using Base::rows;
580 using Base::cols;
581 using Base::rank;
582
583 private:
584 void allocate(Index rows, Index cols, unsigned int computationOptions);
585
586 protected:
587 using Base::m_matrixU;
588 using Base::m_matrixV;
589 using Base::m_singularValues;
590 using Base::m_info;
591 using Base::m_isInitialized;
592 using Base::m_isAllocated;
593 using Base::m_usePrescribedThreshold;
594 using Base::m_computeFullU;
595 using Base::m_computeThinU;
596 using Base::m_computeFullV;
597 using Base::m_computeThinV;
598 using Base::m_computationOptions;
599 using Base::m_nonzeroSingularValues;
600 using Base::m_rows;
601 using Base::m_cols;
602 using Base::m_diagSize;
603 using Base::m_prescribedThreshold;
604 WorkMatrixType m_workMatrix;
605
606 template<typename MatrixType__, int QRPreconditioner_, bool IsComplex_>
607 friend struct internal::svd_precondition_2x2_block_to_be_real;
608 template<typename MatrixType__, int QRPreconditioner_, int Case_, bool DoAnything_>
609 friend struct internal::qr_preconditioner_impl;
610
611 internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreColsThanRows> m_qr_precond_morecols;
612 internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreRowsThanCols> m_qr_precond_morerows;
613 MatrixType m_scaledMatrix;
614};
615
616template<typename MatrixType, int QRPreconditioner>
617void JacobiSVD<MatrixType, QRPreconditioner>::allocate(Eigen::Index rows, Eigen::Index cols, unsigned int computationOptions)
618{
619 eigen_assert(rows >= 0 && cols >= 0);
620
621 if (m_isAllocated &&
622 rows == m_rows &&
623 cols == m_cols &&
624 computationOptions == m_computationOptions)
625 {
626 return;
627 }
628
629 m_rows = rows;
630 m_cols = cols;
631 m_info = Success;
632 m_isInitialized = false;
633 m_isAllocated = true;
634 m_computationOptions = computationOptions;
635 m_computeFullU = (computationOptions & ComputeFullU) != 0;
636 m_computeThinU = (computationOptions & ComputeThinU) != 0;
637 m_computeFullV = (computationOptions & ComputeFullV) != 0;
638 m_computeThinV = (computationOptions & ComputeThinV) != 0;
639 eigen_assert(!(m_computeFullU && m_computeThinU) && "JacobiSVD: you can't ask for both full and thin U");
640 eigen_assert(!(m_computeFullV && m_computeThinV) && "JacobiSVD: you can't ask for both full and thin V");
641 eigen_assert(internal::check_implication(m_computeThinU || m_computeThinV, MatrixType::ColsAtCompileTime==Dynamic) &&
642 "JacobiSVD: thin U and V are only available when your matrix has a dynamic number of columns.");
643 if (QRPreconditioner == FullPivHouseholderQRPreconditioner)
644 {
645 eigen_assert(!(m_computeThinU || m_computeThinV) &&
646 "JacobiSVD: can't compute thin U or thin V with the FullPivHouseholderQR preconditioner. "
647 "Use the ColPivHouseholderQR preconditioner instead.");
648 }
649 m_diagSize = (std::min)(m_rows, m_cols);
650 m_singularValues.resize(m_diagSize);
651 if(RowsAtCompileTime==Dynamic)
652 m_matrixU.resize(m_rows, m_computeFullU ? m_rows
653 : m_computeThinU ? m_diagSize
654 : 0);
655 if(ColsAtCompileTime==Dynamic)
656 m_matrixV.resize(m_cols, m_computeFullV ? m_cols
657 : m_computeThinV ? m_diagSize
658 : 0);
659 m_workMatrix.resize(m_diagSize, m_diagSize);
660
661 if(m_cols>m_rows) m_qr_precond_morecols.allocate(*this);
662 if(m_rows>m_cols) m_qr_precond_morerows.allocate(*this);
663 if(m_rows!=m_cols) m_scaledMatrix.resize(rows,cols);
664}
665
666template<typename MatrixType, int QRPreconditioner>
667JacobiSVD<MatrixType, QRPreconditioner>&
668JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsigned int computationOptions)
669{
670 using std::abs;
671 allocate(matrix.rows(), matrix.cols(), computationOptions);
672
673 // currently we stop when we reach precision 2*epsilon as the last bit of precision can require an unreasonable number of iterations,
674 // only worsening the precision of U and V as we accumulate more rotations
675 const RealScalar precision = RealScalar(2) * NumTraits<Scalar>::epsilon();
676
677 // limit for denormal numbers to be considered zero in order to avoid infinite loops (see bug 286)
678 const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)();
679
680 // Scaling factor to reduce over/under-flows
681 RealScalar scale = matrix.cwiseAbs().template maxCoeff<PropagateNaN>();
682 if (!(numext::isfinite)(scale)) {
683 m_isInitialized = true;
684 m_info = InvalidInput;
685 return *this;
686 }
687 if(scale==RealScalar(0)) scale = RealScalar(1);
688
689 /*** step 1. The R-SVD step: we use a QR decomposition to reduce to the case of a square matrix */
690
691 if(m_rows!=m_cols)
692 {
693 m_scaledMatrix = matrix / scale;
694 m_qr_precond_morecols.run(*this, m_scaledMatrix);
695 m_qr_precond_morerows.run(*this, m_scaledMatrix);
696 }
697 else
698 {
699 m_workMatrix = matrix.block(0,0,m_diagSize,m_diagSize) / scale;
700 if(m_computeFullU) m_matrixU.setIdentity(m_rows,m_rows);
701 if(m_computeThinU) m_matrixU.setIdentity(m_rows,m_diagSize);
702 if(m_computeFullV) m_matrixV.setIdentity(m_cols,m_cols);
703 if(m_computeThinV) m_matrixV.setIdentity(m_cols, m_diagSize);
704 }
705
706 /*** step 2. The main Jacobi SVD iteration. ***/
707 RealScalar maxDiagEntry = m_workMatrix.cwiseAbs().diagonal().maxCoeff();
708
709 bool finished = false;
710 while(!finished)
711 {
712 finished = true;
713
714 // do a sweep: for all index pairs (p,q), perform SVD of the corresponding 2x2 sub-matrix
715
716 for(Index p = 1; p < m_diagSize; ++p)
717 {
718 for(Index q = 0; q < p; ++q)
719 {
720 // if this 2x2 sub-matrix is not diagonal already...
721 // notice that this comparison will evaluate to false if any NaN is involved, ensuring that NaN's don't
722 // keep us iterating forever. Similarly, small denormal numbers are considered zero.
723 RealScalar threshold = numext::maxi<RealScalar>(considerAsZero, precision * maxDiagEntry);
724 if(abs(m_workMatrix.coeff(p,q))>threshold || abs(m_workMatrix.coeff(q,p)) > threshold)
725 {
726 finished = false;
727 // perform SVD decomposition of 2x2 sub-matrix corresponding to indices p,q to make it diagonal
728 // the complex to real operation returns true if the updated 2x2 block is not already diagonal
729 if(internal::svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner>::run(m_workMatrix, *this, p, q, maxDiagEntry))
730 {
731 JacobiRotation<RealScalar> j_left, j_right;
732 internal::real_2x2_jacobi_svd(m_workMatrix, p, q, &j_left, &j_right);
733
734 // accumulate resulting Jacobi rotations
735 m_workMatrix.applyOnTheLeft(p,q,j_left);
736 if(computeU()) m_matrixU.applyOnTheRight(p,q,j_left.transpose());
737
738 m_workMatrix.applyOnTheRight(p,q,j_right);
739 if(computeV()) m_matrixV.applyOnTheRight(p,q,j_right);
740
741 // keep track of the largest diagonal coefficient
742 maxDiagEntry = numext::maxi<RealScalar>(maxDiagEntry,numext::maxi<RealScalar>(abs(m_workMatrix.coeff(p,p)), abs(m_workMatrix.coeff(q,q))));
743 }
744 }
745 }
746 }
747 }
748
749 /*** step 3. The work matrix is now diagonal, so ensure it's positive so its diagonal entries are the singular values ***/
750
751 for(Index i = 0; i < m_diagSize; ++i)
752 {
753 // For a complex matrix, some diagonal coefficients might note have been
754 // treated by svd_precondition_2x2_block_to_be_real, and the imaginary part
755 // of some diagonal entry might not be null.
756 if(NumTraits<Scalar>::IsComplex && abs(numext::imag(m_workMatrix.coeff(i,i)))>considerAsZero)
757 {
758 RealScalar a = abs(m_workMatrix.coeff(i,i));
759 m_singularValues.coeffRef(i) = abs(a);
760 if(computeU()) m_matrixU.col(i) *= m_workMatrix.coeff(i,i)/a;
761 }
762 else
763 {
764 // m_workMatrix.coeff(i,i) is already real, no difficulty:
765 RealScalar a = numext::real(m_workMatrix.coeff(i,i));
766 m_singularValues.coeffRef(i) = abs(a);
767 if(computeU() && (a<RealScalar(0))) m_matrixU.col(i) = -m_matrixU.col(i);
768 }
769 }
770
771 m_singularValues *= scale;
772
773 /*** step 4. Sort singular values in descending order and compute the number of nonzero singular values ***/
774
775 m_nonzeroSingularValues = m_diagSize;
776 for(Index i = 0; i < m_diagSize; i++)
777 {
778 Index pos;
779 RealScalar maxRemainingSingularValue = m_singularValues.tail(m_diagSize-i).maxCoeff(&pos);
780 if(maxRemainingSingularValue == RealScalar(0))
781 {
782 m_nonzeroSingularValues = i;
783 break;
784 }
785 if(pos)
786 {
787 pos += i;
788 std::swap(m_singularValues.coeffRef(i), m_singularValues.coeffRef(pos));
789 if(computeU()) m_matrixU.col(pos).swap(m_matrixU.col(i));
790 if(computeV()) m_matrixV.col(pos).swap(m_matrixV.col(i));
791 }
792 }
793
794 m_isInitialized = true;
795 return *this;
796}
797
805template<typename Derived>
807MatrixBase<Derived>::jacobiSvd(unsigned int computationOptions) const
808{
809 return JacobiSVD<PlainObject>(*this, computationOptions);
810}
811
812} // end namespace Eigen
813
814#endif // EIGEN_JACOBISVD_H
internal::traits< Homogeneous< MatrixType, Direction_ > >::Scalar Scalar
Definition: DenseBase.h:61
Rotation given by a cosine-sine pair.
Definition: Jacobi.h:37
Definition: Jacobi.h:69
JacobiRotation transpose() const
Definition: Jacobi.h:65
Two-sided Jacobi SVD decomposition of a rectangular matrix.
Definition: JacobiSVD.h:492
EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT
Definition: EigenBase.h:65
JacobiSVD & compute(const MatrixType &matrix, unsigned int computationOptions)
Method performing the decomposition of given matrix using custom options.
Definition: JacobiSVD.h:668
JacobiSVD(Index rows, Index cols, unsigned int computationOptions=0)
Default Constructor with memory preallocation.
Definition: JacobiSVD.h:534
EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: EigenBase.h:62
JacobiSVD & compute(const MatrixType &matrix)
Method performing the decomposition of given matrix using current options.
Definition: JacobiSVD.h:572
JacobiSVD()
Default Constructor.
Definition: JacobiSVD.h:524
JacobiSVD(const MatrixType &matrix, unsigned int computationOptions=0)
Constructor performing the decomposition of given matrix.
Definition: JacobiSVD.h:549
JacobiSVD< PlainObject > jacobiSvd(unsigned int computationOptions=0) const
Definition: JacobiSVD.h:807
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:182
Base class of SVD algorithms.
Definition: SVDBase.h:66
Index rank() const
Definition: SVDBase.h:150
bool computeV() const
Definition: SVDBase.h:212
bool computeU() const
Definition: SVDBase.h:210
@ NoQRPreconditioner
Definition: Constants.h:427
@ HouseholderQRPreconditioner
Definition: Constants.h:429
@ ColPivHouseholderQRPreconditioner
Definition: Constants.h:431
@ FullPivHouseholderQRPreconditioner
Definition: Constants.h:433
@ InvalidInput
Definition: Constants.h:451
@ Success
Definition: Constants.h:444
@ ComputeFullV
Definition: Constants.h:399
@ ComputeThinV
Definition: Constants.h:401
@ ComputeFullU
Definition: Constants.h:395
@ ComputeThinU
Definition: Constants.h:397
Namespace containing all symbols from the Eigen library.
Definition: B01_Experimental.dox:1
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_sqrt_op< typename Derived::Scalar >, const Derived > sqrt(const Eigen::ArrayBase< Derived > &x)
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_conjugate_op< typename Derived::Scalar >, const Derived > conj(const Eigen::ArrayBase< Derived > &x)
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:59
const int Dynamic
Definition: Constants.h:24
Definition: EigenBase.h:32
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:41
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:235