Eigen  5.0.1-dev+60122df6
 
Loading...
Searching...
No Matches
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
14// IWYU pragma: private
15#include "./InternalHeaderCheck.h"
16
17namespace Eigen {
18
19namespace internal {
20
21// forward declaration (needed by ICC)
22// the empty body is required by MSVC
23template <typename MatrixType, int Options, bool IsComplex = NumTraits<typename MatrixType::Scalar>::IsComplex>
24struct svd_precondition_2x2_block_to_be_real {};
25
26/*** QR preconditioners (R-SVD)
27 ***
28 *** Their role is to reduce the problem of computing the SVD to the case of a square matrix.
29 *** This approach, known as R-SVD, is an optimization for rectangular-enough matrices, and is a requirement for
30 *** JacobiSVD which by itself is only able to work on square matrices.
31 ***/
32
33enum { PreconditionIfMoreColsThanRows, PreconditionIfMoreRowsThanCols };
34
35template <typename MatrixType, int QRPreconditioner, int Case>
36struct qr_preconditioner_should_do_anything {
37 enum {
38 a = MatrixType::RowsAtCompileTime != Dynamic && MatrixType::ColsAtCompileTime != Dynamic &&
39 MatrixType::ColsAtCompileTime <= MatrixType::RowsAtCompileTime,
40 b = MatrixType::RowsAtCompileTime != Dynamic && MatrixType::ColsAtCompileTime != Dynamic &&
41 MatrixType::RowsAtCompileTime <= MatrixType::ColsAtCompileTime,
42 ret = !((QRPreconditioner == NoQRPreconditioner) || (Case == PreconditionIfMoreColsThanRows && bool(a)) ||
43 (Case == PreconditionIfMoreRowsThanCols && bool(b)))
44 };
45};
46
47template <typename MatrixType, int Options, int QRPreconditioner, int Case,
48 bool DoAnything = qr_preconditioner_should_do_anything<MatrixType, QRPreconditioner, Case>::ret>
49struct qr_preconditioner_impl {};
50
51template <typename MatrixType, int Options, int QRPreconditioner, int Case>
52class qr_preconditioner_impl<MatrixType, Options, QRPreconditioner, Case, false> {
53 public:
54 void allocate(const JacobiSVD<MatrixType, Options>&) {}
55 template <typename Xpr>
56 bool run(JacobiSVD<MatrixType, Options>&, const Xpr&) {
57 return false;
58 }
59};
60
61/*** preconditioner using FullPivHouseholderQR ***/
62
63template <typename MatrixType, int Options>
64class qr_preconditioner_impl<MatrixType, Options, FullPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols,
65 true> {
66 public:
67 typedef typename MatrixType::Scalar Scalar;
68 typedef JacobiSVD<MatrixType, Options> SVDType;
69
70 enum { WorkspaceSize = MatrixType::RowsAtCompileTime, MaxWorkspaceSize = MatrixType::MaxRowsAtCompileTime };
71
72 typedef Matrix<Scalar, 1, WorkspaceSize, RowMajor, 1, MaxWorkspaceSize> WorkspaceType;
73
74 void allocate(const SVDType& svd) {
75 if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols()) {
76 internal::destroy_at(&m_qr);
77 internal::construct_at(&m_qr, svd.rows(), svd.cols());
78 }
79 if (svd.m_computeFullU) m_workspace.resize(svd.rows());
80 }
81 template <typename Xpr>
82 bool run(SVDType& svd, const Xpr& matrix) {
83 if (matrix.rows() > matrix.cols()) {
84 m_qr.compute(matrix);
85 svd.m_workMatrix = m_qr.matrixQR().block(0, 0, matrix.cols(), matrix.cols()).template triangularView<Upper>();
86 if (svd.m_computeFullU) m_qr.matrixQ().evalTo(svd.m_matrixU, m_workspace);
87 if (svd.computeV()) svd.m_matrixV = m_qr.colsPermutation();
88 return true;
89 }
90 return false;
91 }
92
93 private:
94 typedef FullPivHouseholderQR<MatrixType> QRType;
95 QRType m_qr;
96 WorkspaceType m_workspace;
97};
98
99template <typename MatrixType, int Options>
100class qr_preconditioner_impl<MatrixType, Options, FullPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows,
101 true> {
102 public:
103 typedef typename MatrixType::Scalar Scalar;
104 typedef JacobiSVD<MatrixType, Options> SVDType;
105
106 enum {
107 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
108 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
109 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
110 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
111 MatrixOptions = traits<MatrixType>::Options
112 };
113
114 typedef typename internal::make_proper_matrix_type<Scalar, ColsAtCompileTime, RowsAtCompileTime, MatrixOptions,
115 MaxColsAtCompileTime, MaxRowsAtCompileTime>::type
116 TransposeTypeWithSameStorageOrder;
117
118 void allocate(const SVDType& svd) {
119 if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols()) {
120 internal::destroy_at(&m_qr);
121 internal::construct_at(&m_qr, svd.cols(), svd.rows());
122 }
123 if (svd.m_computeFullV) m_workspace.resize(svd.cols());
124 }
125 template <typename Xpr>
126 bool run(SVDType& svd, const Xpr& matrix) {
127 if (matrix.cols() > matrix.rows()) {
128 m_qr.compute(matrix.adjoint());
129 svd.m_workMatrix =
130 m_qr.matrixQR().block(0, 0, matrix.rows(), matrix.rows()).template triangularView<Upper>().adjoint();
131 if (svd.m_computeFullV) m_qr.matrixQ().evalTo(svd.m_matrixV, m_workspace);
132 if (svd.computeU()) svd.m_matrixU = m_qr.colsPermutation();
133 return true;
134 } else
135 return false;
136 }
137
138 private:
139 typedef FullPivHouseholderQR<TransposeTypeWithSameStorageOrder> QRType;
140 QRType m_qr;
141 typename plain_row_type<MatrixType>::type m_workspace;
142};
143
144/*** preconditioner using ColPivHouseholderQR ***/
145
146template <typename MatrixType, int Options>
147class qr_preconditioner_impl<MatrixType, Options, ColPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols,
148 true> {
149 public:
150 typedef typename MatrixType::Scalar Scalar;
151 typedef JacobiSVD<MatrixType, Options> SVDType;
152
153 enum {
154 WorkspaceSize = internal::traits<SVDType>::MatrixUColsAtCompileTime,
155 MaxWorkspaceSize = internal::traits<SVDType>::MatrixUMaxColsAtCompileTime
156 };
157
158 typedef Matrix<Scalar, 1, WorkspaceSize, RowMajor, 1, MaxWorkspaceSize> WorkspaceType;
159
160 void allocate(const SVDType& svd) {
161 if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols()) {
162 internal::destroy_at(&m_qr);
163 internal::construct_at(&m_qr, svd.rows(), svd.cols());
164 }
165 if (svd.m_computeFullU)
166 m_workspace.resize(svd.rows());
167 else if (svd.m_computeThinU)
168 m_workspace.resize(svd.cols());
169 }
170 template <typename Xpr>
171 bool run(SVDType& svd, const Xpr& matrix) {
172 if (matrix.rows() > matrix.cols()) {
173 m_qr.compute(matrix);
174 svd.m_workMatrix = m_qr.matrixQR().block(0, 0, matrix.cols(), matrix.cols()).template triangularView<Upper>();
175 if (svd.m_computeFullU)
176 m_qr.householderQ().evalTo(svd.m_matrixU, m_workspace);
177 else if (svd.m_computeThinU) {
178 svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols());
179 m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_workspace);
180 }
181 if (svd.computeV()) svd.m_matrixV = m_qr.colsPermutation();
182 return true;
183 }
184 return false;
185 }
186
187 private:
188 typedef ColPivHouseholderQR<MatrixType> QRType;
189 QRType m_qr;
190 WorkspaceType m_workspace;
191};
192
193template <typename MatrixType, int Options>
194class qr_preconditioner_impl<MatrixType, Options, ColPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows,
195 true> {
196 public:
197 typedef typename MatrixType::Scalar Scalar;
198 typedef JacobiSVD<MatrixType, Options> SVDType;
199
200 enum {
201 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
202 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
203 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
204 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
205 MatrixOptions = internal::traits<MatrixType>::Options,
206 WorkspaceSize = internal::traits<SVDType>::MatrixVColsAtCompileTime,
207 MaxWorkspaceSize = internal::traits<SVDType>::MatrixVMaxColsAtCompileTime
208 };
209
210 typedef Matrix<Scalar, WorkspaceSize, 1, ColMajor, MaxWorkspaceSize, 1> WorkspaceType;
211
212 typedef typename internal::make_proper_matrix_type<Scalar, ColsAtCompileTime, RowsAtCompileTime, MatrixOptions,
213 MaxColsAtCompileTime, MaxRowsAtCompileTime>::type
214 TransposeTypeWithSameStorageOrder;
215
216 void allocate(const SVDType& svd) {
217 if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols()) {
218 internal::destroy_at(&m_qr);
219 internal::construct_at(&m_qr, svd.cols(), svd.rows());
220 }
221 if (svd.m_computeFullV)
222 m_workspace.resize(svd.cols());
223 else if (svd.m_computeThinV)
224 m_workspace.resize(svd.rows());
225 }
226 template <typename Xpr>
227 bool run(SVDType& svd, const Xpr& matrix) {
228 if (matrix.cols() > matrix.rows()) {
229 m_qr.compute(matrix.adjoint());
230
231 svd.m_workMatrix =
232 m_qr.matrixQR().block(0, 0, matrix.rows(), matrix.rows()).template triangularView<Upper>().adjoint();
233 if (svd.m_computeFullV)
234 m_qr.householderQ().evalTo(svd.m_matrixV, m_workspace);
235 else if (svd.m_computeThinV) {
236 svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows());
237 m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_workspace);
238 }
239 if (svd.computeU()) svd.m_matrixU = m_qr.colsPermutation();
240 return true;
241 } else
242 return false;
243 }
244
245 private:
246 typedef ColPivHouseholderQR<TransposeTypeWithSameStorageOrder> QRType;
247 QRType m_qr;
248 WorkspaceType m_workspace;
249};
250
251/*** preconditioner using HouseholderQR ***/
252
253template <typename MatrixType, int Options>
254class qr_preconditioner_impl<MatrixType, Options, HouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true> {
255 public:
256 typedef typename MatrixType::Scalar Scalar;
257 typedef JacobiSVD<MatrixType, Options> SVDType;
258
259 enum {
260 WorkspaceSize = internal::traits<SVDType>::MatrixUColsAtCompileTime,
261 MaxWorkspaceSize = internal::traits<SVDType>::MatrixUMaxColsAtCompileTime
262 };
263
264 typedef Matrix<Scalar, 1, WorkspaceSize, RowMajor, 1, MaxWorkspaceSize> WorkspaceType;
265
266 void allocate(const SVDType& svd) {
267 if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols()) {
268 internal::destroy_at(&m_qr);
269 internal::construct_at(&m_qr, svd.rows(), svd.cols());
270 }
271 if (svd.m_computeFullU)
272 m_workspace.resize(svd.rows());
273 else if (svd.m_computeThinU)
274 m_workspace.resize(svd.cols());
275 }
276 template <typename Xpr>
277 bool run(SVDType& svd, const Xpr& matrix) {
278 if (matrix.rows() > matrix.cols()) {
279 m_qr.compute(matrix);
280 svd.m_workMatrix = m_qr.matrixQR().block(0, 0, matrix.cols(), matrix.cols()).template triangularView<Upper>();
281 if (svd.m_computeFullU)
282 m_qr.householderQ().evalTo(svd.m_matrixU, m_workspace);
283 else if (svd.m_computeThinU) {
284 svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols());
285 m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_workspace);
286 }
287 if (svd.computeV()) svd.m_matrixV.setIdentity(matrix.cols(), matrix.cols());
288 return true;
289 }
290 return false;
291 }
292
293 private:
294 typedef HouseholderQR<MatrixType> QRType;
295 QRType m_qr;
296 WorkspaceType m_workspace;
297};
298
299template <typename MatrixType, int Options>
300class qr_preconditioner_impl<MatrixType, Options, HouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true> {
301 public:
302 typedef typename MatrixType::Scalar Scalar;
303 typedef JacobiSVD<MatrixType, Options> SVDType;
304
305 enum {
306 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
307 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
308 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
309 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
310 MatrixOptions = internal::traits<MatrixType>::Options,
311 WorkspaceSize = internal::traits<SVDType>::MatrixVColsAtCompileTime,
312 MaxWorkspaceSize = internal::traits<SVDType>::MatrixVMaxColsAtCompileTime
313 };
314
315 typedef Matrix<Scalar, WorkspaceSize, 1, ColMajor, MaxWorkspaceSize, 1> WorkspaceType;
316
317 typedef typename internal::make_proper_matrix_type<Scalar, ColsAtCompileTime, RowsAtCompileTime, MatrixOptions,
318 MaxColsAtCompileTime, MaxRowsAtCompileTime>::type
319 TransposeTypeWithSameStorageOrder;
320
321 void allocate(const SVDType& svd) {
322 if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols()) {
323 internal::destroy_at(&m_qr);
324 internal::construct_at(&m_qr, svd.cols(), svd.rows());
325 }
326 if (svd.m_computeFullV)
327 m_workspace.resize(svd.cols());
328 else if (svd.m_computeThinV)
329 m_workspace.resize(svd.rows());
330 }
331
332 template <typename Xpr>
333 bool run(SVDType& svd, const Xpr& matrix) {
334 if (matrix.cols() > matrix.rows()) {
335 m_qr.compute(matrix.adjoint());
336
337 svd.m_workMatrix =
338 m_qr.matrixQR().block(0, 0, matrix.rows(), matrix.rows()).template triangularView<Upper>().adjoint();
339 if (svd.m_computeFullV)
340 m_qr.householderQ().evalTo(svd.m_matrixV, m_workspace);
341 else if (svd.m_computeThinV) {
342 svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows());
343 m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_workspace);
344 }
345 if (svd.computeU()) svd.m_matrixU.setIdentity(matrix.rows(), matrix.rows());
346 return true;
347 } else
348 return false;
349 }
350
351 private:
352 typedef HouseholderQR<TransposeTypeWithSameStorageOrder> QRType;
353 QRType m_qr;
354 WorkspaceType m_workspace;
355};
356
357/*** 2x2 SVD implementation
358 ***
359 *** JacobiSVD consists in performing a series of 2x2 SVD subproblems
360 ***/
361
362template <typename MatrixType, int Options>
363struct svd_precondition_2x2_block_to_be_real<MatrixType, Options, false> {
364 typedef JacobiSVD<MatrixType, Options> SVD;
365 typedef typename MatrixType::RealScalar RealScalar;
366 static bool run(typename SVD::WorkMatrixType&, SVD&, Index, Index, RealScalar&) { return true; }
367};
368
369template <typename MatrixType, int Options>
370struct svd_precondition_2x2_block_to_be_real<MatrixType, Options, true> {
371 typedef JacobiSVD<MatrixType, Options> SVD;
372 typedef typename MatrixType::Scalar Scalar;
373 typedef typename MatrixType::RealScalar RealScalar;
374 static bool run(typename SVD::WorkMatrixType& work_matrix, SVD& svd, Index p, Index q, RealScalar& maxDiagEntry) {
375 using std::abs;
376 using std::sqrt;
377 Scalar z;
378 JacobiRotation<Scalar> rot;
379 RealScalar n = sqrt(numext::abs2(work_matrix.coeff(p, p)) + numext::abs2(work_matrix.coeff(q, p)));
380
381 const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)();
382 const RealScalar precision = NumTraits<Scalar>::epsilon();
383
384 if (numext::is_exactly_zero(n)) {
385 // make sure first column is zero
386 work_matrix.coeffRef(p, p) = work_matrix.coeffRef(q, p) = Scalar(0);
387
388 if (abs(numext::imag(work_matrix.coeff(p, q))) > considerAsZero) {
389 // work_matrix.coeff(p,q) can be zero if work_matrix.coeff(q,p) is not zero but small enough to underflow when
390 // computing n
391 z = abs(work_matrix.coeff(p, q)) / work_matrix.coeff(p, q);
392 work_matrix.row(p) *= z;
393 if (svd.computeU()) svd.m_matrixU.col(p) *= conj(z);
394 }
395 if (abs(numext::imag(work_matrix.coeff(q, q))) > considerAsZero) {
396 z = abs(work_matrix.coeff(q, q)) / work_matrix.coeff(q, q);
397 work_matrix.row(q) *= z;
398 if (svd.computeU()) svd.m_matrixU.col(q) *= conj(z);
399 }
400 // otherwise the second row is already zero, so we have nothing to do.
401 } else {
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);
405 if (svd.computeU()) svd.m_matrixU.applyOnTheRight(p, q, rot.adjoint());
406 if (abs(numext::imag(work_matrix.coeff(p, q))) > considerAsZero) {
407 z = abs(work_matrix.coeff(p, q)) / work_matrix.coeff(p, q);
408 work_matrix.col(q) *= z;
409 if (svd.computeV()) svd.m_matrixV.col(q) *= z;
410 }
411 if (abs(numext::imag(work_matrix.coeff(q, q))) > considerAsZero) {
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>(
420 maxDiagEntry, numext::maxi<RealScalar>(abs(work_matrix.coeff(p, p)), abs(work_matrix.coeff(q, q))));
421 // and check whether the 2x2 block is already diagonal
422 RealScalar threshold = numext::maxi<RealScalar>(considerAsZero, precision * maxDiagEntry);
423 return abs(work_matrix.coeff(p, q)) > threshold || abs(work_matrix.coeff(q, p)) > threshold;
424 }
425};
426
427template <typename MatrixType_, int Options>
428struct traits<JacobiSVD<MatrixType_, Options> > : svd_traits<MatrixType_, Options> {
429 typedef MatrixType_ MatrixType;
430};
431
432} // end namespace internal
433
499template <typename MatrixType_, int Options_>
500class JacobiSVD : public SVDBase<JacobiSVD<MatrixType_, Options_> > {
501 typedef SVDBase<JacobiSVD> Base;
502
503 public:
504 typedef MatrixType_ MatrixType;
505 typedef typename Base::Scalar Scalar;
506 typedef typename Base::RealScalar RealScalar;
507 enum : int {
508 Options = Options_,
509 QRPreconditioner = internal::get_qr_preconditioner(Options),
510 RowsAtCompileTime = Base::RowsAtCompileTime,
511 ColsAtCompileTime = Base::ColsAtCompileTime,
512 DiagSizeAtCompileTime = Base::DiagSizeAtCompileTime,
513 MaxRowsAtCompileTime = Base::MaxRowsAtCompileTime,
514 MaxColsAtCompileTime = Base::MaxColsAtCompileTime,
515 MaxDiagSizeAtCompileTime = Base::MaxDiagSizeAtCompileTime,
516 MatrixOptions = Base::MatrixOptions
517 };
518
519 typedef typename Base::MatrixUType MatrixUType;
520 typedef typename Base::MatrixVType MatrixVType;
521 typedef typename Base::SingularValuesType SingularValuesType;
522 typedef Matrix<Scalar, DiagSizeAtCompileTime, DiagSizeAtCompileTime, MatrixOptions, MaxDiagSizeAtCompileTime,
523 MaxDiagSizeAtCompileTime>
524 WorkMatrixType;
525
531 JacobiSVD() {}
532
540 JacobiSVD(Index rows, Index cols) { allocate(rows, cols, internal::get_computation_options(Options)); }
541
558 EIGEN_DEPRECATED_WITH_REASON("Options should be specified using the class template parameter.")
559 JacobiSVD(Index rows, Index cols, unsigned int computationOptions) {
560 internal::check_svd_options_assertions<MatrixType, Options>(computationOptions, rows, cols);
561 allocate(rows, cols, computationOptions);
562 }
563
569 template <typename Derived>
570 explicit JacobiSVD(const MatrixBase<Derived>& matrix) {
571 compute_impl(matrix, internal::get_computation_options(Options));
572 }
573
574 template <typename Derived>
575 explicit JacobiSVD(const TriangularBase<Derived>& matrix) {
576 compute_impl(matrix, internal::get_computation_options(Options));
577 }
578
591 // EIGEN_DEPRECATED // TODO(cantonios): re-enable after fixing a few 3p libraries that error on deprecation warnings.
592 template <typename Derived>
593 JacobiSVD(const MatrixBase<Derived>& matrix, unsigned int computationOptions) {
594 internal::check_svd_options_assertions<MatrixBase<Derived>, Options>(computationOptions, matrix.rows(),
595 matrix.cols());
596 compute_impl(matrix, computationOptions);
597 }
598
604 template <typename Derived>
605 JacobiSVD& compute(const MatrixBase<Derived>& matrix) {
606 return compute_impl(matrix, m_computationOptions);
607 }
608
609 template <typename Derived>
610 JacobiSVD& compute(const TriangularBase<Derived>& matrix) {
611 return compute_impl(matrix, m_computationOptions);
612 }
613
623 template <typename Derived>
624 EIGEN_DEPRECATED_WITH_REASON("Options should be specified using the class template parameter.")
625 JacobiSVD& compute(const MatrixBase<Derived>& matrix, unsigned int computationOptions) {
626 internal::check_svd_options_assertions<MatrixBase<Derived>, Options>(m_computationOptions, matrix.rows(),
627 matrix.cols());
628 return compute_impl(matrix, computationOptions);
629 }
630
631 using Base::cols;
632 using Base::computeU;
633 using Base::computeV;
634 using Base::diagSize;
635 using Base::rank;
636 using Base::rows;
637
638 void allocate(Index rows_, Index cols_, unsigned int computationOptions) {
639 if (Base::allocate(rows_, cols_, computationOptions)) return;
640 eigen_assert(!(ShouldComputeThinU && int(QRPreconditioner) == int(FullPivHouseholderQRPreconditioner)) &&
641 !(ShouldComputeThinU && int(QRPreconditioner) == int(FullPivHouseholderQRPreconditioner)) &&
642 "JacobiSVD: can't compute thin U or thin V with the FullPivHouseholderQR preconditioner. "
643 "Use the ColPivHouseholderQR preconditioner instead.");
644
645 m_workMatrix.resize(diagSize(), diagSize());
646 if (cols() > rows()) m_qr_precond_morecols.allocate(*this);
647 if (rows() > cols()) m_qr_precond_morerows.allocate(*this);
648 }
649
650 private:
651 template <typename Derived>
652 JacobiSVD& compute_impl(const TriangularBase<Derived>& matrix, unsigned int computationOptions);
653 template <typename Derived>
654 JacobiSVD& compute_impl(const MatrixBase<Derived>& matrix, unsigned int computationOptions);
655
656 protected:
657 using Base::m_computationOptions;
658 using Base::m_computeFullU;
659 using Base::m_computeFullV;
660 using Base::m_computeThinU;
661 using Base::m_computeThinV;
662 using Base::m_info;
663 using Base::m_isAllocated;
664 using Base::m_isInitialized;
665 using Base::m_matrixU;
666 using Base::m_matrixV;
667 using Base::m_nonzeroSingularValues;
668 using Base::m_prescribedThreshold;
669 using Base::m_singularValues;
670 using Base::m_usePrescribedThreshold;
671 using Base::ShouldComputeThinU;
672 using Base::ShouldComputeThinV;
673
674 EIGEN_STATIC_ASSERT(!(ShouldComputeThinU && int(QRPreconditioner) == int(FullPivHouseholderQRPreconditioner)) &&
675 !(ShouldComputeThinU && int(QRPreconditioner) == int(FullPivHouseholderQRPreconditioner)),
676 "JacobiSVD: can't compute thin U or thin V with the FullPivHouseholderQR preconditioner. "
677 "Use the ColPivHouseholderQR preconditioner instead.")
678
679 template <typename MatrixType__, int Options__, bool IsComplex_>
680 friend struct internal::svd_precondition_2x2_block_to_be_real;
681 template <typename MatrixType__, int Options__, int QRPreconditioner_, int Case_, bool DoAnything_>
682 friend struct internal::qr_preconditioner_impl;
683
684 internal::qr_preconditioner_impl<MatrixType, Options, QRPreconditioner, internal::PreconditionIfMoreColsThanRows>
685 m_qr_precond_morecols;
686 internal::qr_preconditioner_impl<MatrixType, Options, QRPreconditioner, internal::PreconditionIfMoreRowsThanCols>
687 m_qr_precond_morerows;
688 WorkMatrixType m_workMatrix;
689};
690
691template <typename MatrixType, int Options>
692template <typename Derived>
693JacobiSVD<MatrixType, Options>& JacobiSVD<MatrixType, Options>::compute_impl(const TriangularBase<Derived>& matrix,
694 unsigned int computationOptions) {
695 return compute_impl(matrix.toDenseMatrix(), computationOptions);
696}
697
698template <typename MatrixType, int Options>
699template <typename Derived>
700JacobiSVD<MatrixType, Options>& JacobiSVD<MatrixType, Options>::compute_impl(const MatrixBase<Derived>& matrix,
701 unsigned int computationOptions) {
702 EIGEN_STATIC_ASSERT_SAME_MATRIX_SIZE(Derived, MatrixType);
703 EIGEN_STATIC_ASSERT((std::is_same<typename Derived::Scalar, typename MatrixType::Scalar>::value),
704 Input matrix must have the same Scalar type as the BDCSVD object.);
705
706 using std::abs;
707
708 allocate(matrix.rows(), matrix.cols(), computationOptions);
709
710 // currently we stop when we reach precision 2*epsilon as the last bit of precision can require an unreasonable number
711 // of iterations, only worsening the precision of U and V as we accumulate more rotations
712 const RealScalar precision = RealScalar(2) * NumTraits<Scalar>::epsilon();
713
714 // limit for denormal numbers to be considered zero in order to avoid infinite loops (see bug 286)
715 const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)();
716
717 // Scaling factor to reduce over/under-flows
718 RealScalar scale = matrix.cwiseAbs().template maxCoeff<PropagateNaN>();
719 if (!(numext::isfinite)(scale)) {
720 m_isInitialized = true;
721 m_info = InvalidInput;
722 m_nonzeroSingularValues = 0;
723 return *this;
724 }
725 if (numext::is_exactly_zero(scale)) scale = RealScalar(1);
726
727 /*** step 1. The R-SVD step: we use a QR decomposition to reduce to the case of a square matrix */
728
729 if (rows() != cols()) {
730 m_qr_precond_morecols.run(*this, matrix / scale);
731 m_qr_precond_morerows.run(*this, matrix / scale);
732 } else {
733 m_workMatrix =
734 matrix.template topLeftCorner<DiagSizeAtCompileTime, DiagSizeAtCompileTime>(diagSize(), diagSize()) / scale;
735 if (m_computeFullU) m_matrixU.setIdentity(rows(), rows());
736 if (m_computeThinU) m_matrixU.setIdentity(rows(), diagSize());
737 if (m_computeFullV) m_matrixV.setIdentity(cols(), cols());
738 if (m_computeThinV) m_matrixV.setIdentity(cols(), diagSize());
739 }
740
741 /*** step 2. The main Jacobi SVD iteration. ***/
742 RealScalar maxDiagEntry = m_workMatrix.cwiseAbs().diagonal().maxCoeff();
743
744 bool finished = false;
745 while (!finished) {
746 finished = true;
747
748 // do a sweep: for all index pairs (p,q), perform SVD of the corresponding 2x2 sub-matrix
749
750 for (Index p = 1; p < diagSize(); ++p) {
751 for (Index q = 0; q < p; ++q) {
752 // if this 2x2 sub-matrix is not diagonal already...
753 // notice that this comparison will evaluate to false if any NaN is involved, ensuring that NaN's don't
754 // keep us iterating forever. Similarly, small denormal numbers are considered zero.
755 RealScalar threshold = numext::maxi<RealScalar>(considerAsZero, precision * maxDiagEntry);
756 if (abs(m_workMatrix.coeff(p, q)) > threshold || abs(m_workMatrix.coeff(q, p)) > threshold) {
757 finished = false;
758 // perform SVD decomposition of 2x2 sub-matrix corresponding to indices p,q to make it diagonal
759 // the complex to real operation returns true if the updated 2x2 block is not already diagonal
760 if (internal::svd_precondition_2x2_block_to_be_real<MatrixType, Options>::run(m_workMatrix, *this, p, q,
761 maxDiagEntry)) {
762 JacobiRotation<RealScalar> j_left, j_right;
763 internal::real_2x2_jacobi_svd(m_workMatrix, p, q, &j_left, &j_right);
764
765 // accumulate resulting Jacobi rotations
766 m_workMatrix.applyOnTheLeft(p, q, j_left);
767 if (computeU()) m_matrixU.applyOnTheRight(p, q, j_left.transpose());
768
769 m_workMatrix.applyOnTheRight(p, q, j_right);
770 if (computeV()) m_matrixV.applyOnTheRight(p, q, j_right);
771
772 // keep track of the largest diagonal coefficient
773 maxDiagEntry = numext::maxi<RealScalar>(
774 maxDiagEntry, numext::maxi<RealScalar>(abs(m_workMatrix.coeff(p, p)), abs(m_workMatrix.coeff(q, q))));
775 }
776 }
777 }
778 }
779 }
780
781 /*** step 3. The work matrix is now diagonal, so ensure it's positive so its diagonal entries are the singular values
782 * ***/
783
784 for (Index i = 0; i < diagSize(); ++i) {
785 // For a complex matrix, some diagonal coefficients might note have been
786 // treated by svd_precondition_2x2_block_to_be_real, and the imaginary part
787 // of some diagonal entry might not be null.
788 if (NumTraits<Scalar>::IsComplex && abs(numext::imag(m_workMatrix.coeff(i, i))) > considerAsZero) {
789 RealScalar a = abs(m_workMatrix.coeff(i, i));
790 m_singularValues.coeffRef(i) = abs(a);
791 if (computeU()) m_matrixU.col(i) *= m_workMatrix.coeff(i, i) / a;
792 } else {
793 // m_workMatrix.coeff(i,i) is already real, no difficulty:
794 RealScalar a = numext::real(m_workMatrix.coeff(i, i));
795 m_singularValues.coeffRef(i) = abs(a);
796 if (computeU() && (a < RealScalar(0))) m_matrixU.col(i) = -m_matrixU.col(i);
797 }
798 }
799
800 m_singularValues *= scale;
801
802 /*** step 4. Sort singular values in descending order and compute the number of nonzero singular values ***/
803
804 m_nonzeroSingularValues = diagSize();
805 for (Index i = 0; i < diagSize(); i++) {
806 Index pos;
807 RealScalar maxRemainingSingularValue = m_singularValues.tail(diagSize() - i).maxCoeff(&pos);
808 if (numext::is_exactly_zero(maxRemainingSingularValue)) {
809 m_nonzeroSingularValues = i;
810 break;
811 }
812 if (pos) {
813 pos += i;
814 std::swap(m_singularValues.coeffRef(i), m_singularValues.coeffRef(pos));
815 if (computeU()) m_matrixU.col(pos).swap(m_matrixU.col(i));
816 if (computeV()) m_matrixV.col(pos).swap(m_matrixV.col(i));
817 }
818 }
819
820 m_isInitialized = true;
821 return *this;
822}
823
831template <typename Derived>
832template <int Options>
833JacobiSVD<typename MatrixBase<Derived>::PlainObject, Options> MatrixBase<Derived>::jacobiSvd() const {
834 return JacobiSVD<PlainObject, Options>(*this);
835}
836
837template <typename Derived>
838template <int Options>
839JacobiSVD<typename MatrixBase<Derived>::PlainObject, Options> MatrixBase<Derived>::jacobiSvd(
840 unsigned int computationOptions) const {
841 return JacobiSVD<PlainObject, Options>(*this, computationOptions);
842}
843
844} // end namespace Eigen
845
846#endif // EIGEN_JACOBISVD_H
Eigen::JacobiSVD
Two-sided Jacobi SVD decomposition of a rectangular matrix.
Definition JacobiSVD.h:500
Eigen::JacobiSVD::JacobiSVD
JacobiSVD()
Default Constructor.
Definition JacobiSVD.h:531
Eigen::JacobiSVD::computeV
bool computeV() const
Definition SVDBase.h:275
Eigen::JacobiSVD::compute
JacobiSVD & compute(const MatrixBase< Derived > &matrix)
Method performing the decomposition of given matrix. Computes Thin/Full unitaries U/V if specified us...
Definition JacobiSVD.h:605
Eigen::JacobiSVD::computeU
bool computeU() const
Definition SVDBase.h:273
Eigen::JacobiSVD::JacobiSVD
JacobiSVD(Index rows, Index cols)
Default Constructor with memory preallocation.
Definition JacobiSVD.h:540
Eigen::JacobiSVD::JacobiSVD
JacobiSVD(const MatrixBase< Derived > &matrix)
Constructor performing the decomposition of given matrix, using the custom options specified with the...
Definition JacobiSVD.h:570
Eigen::JacobiSVD::JacobiSVD
JacobiSVD(const MatrixBase< Derived > &matrix, unsigned int computationOptions)
Constructor performing the decomposition of given matrix using specified options for computing unitar...
Definition JacobiSVD.h:593
Eigen::MatrixBase
Base class for all dense matrices, vectors, and expressions.
Definition MatrixBase.h:52
Eigen::Matrix
The matrix class, also used for vectors and row-vectors.
Definition Matrix.h:186
Eigen::PlainObjectBase::resize
constexpr void resize(Index rows, Index cols)
Definition PlainObjectBase.h:268
Eigen::SVDBase< JacobiSVD >::rank
Index rank() const
Definition SVDBase.h:217
Eigen::SVDBase< JacobiSVD >::computeV
bool computeV() const
Definition SVDBase.h:275
Eigen::SVDBase< JacobiSVD >::computeU
bool computeU() const
Definition SVDBase.h:273
Eigen::SVDBase< JacobiSVD< MatrixType_, Options_ > >::threshold
RealScalar threshold() const
Definition SVDBase.h:265
Eigen::TriangularBase
Base class for triangular part in a matrix.
Definition TriangularMatrix.h:32
Eigen::NoQRPreconditioner
@ NoQRPreconditioner
Definition Constants.h:423
Eigen::HouseholderQRPreconditioner
@ HouseholderQRPreconditioner
Definition Constants.h:425
Eigen::ColPivHouseholderQRPreconditioner
@ ColPivHouseholderQRPreconditioner
Definition Constants.h:421
Eigen::FullPivHouseholderQRPreconditioner
@ FullPivHouseholderQRPreconditioner
Definition Constants.h:427
Eigen::InvalidInput
@ InvalidInput
Definition Constants.h:447
Eigen
Namespace containing all symbols from the Eigen library.
Definition B01_Experimental.dox:1
Eigen::sqrt
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_sqrt_op< typename Derived::Scalar >, const Derived > sqrt(const Eigen::ArrayBase< Derived > &x)
Eigen::EIGEN_DEPRECATED_WITH_REASON
EIGEN_DEPRECATED_WITH_REASON("Initialization is no longer needed.") inline void initParallel()
Definition Parallelizer.h:50
Eigen::conj
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_conjugate_op< typename Derived::Scalar >, const Derived > conj(const Eigen::ArrayBase< Derived > &x)
Eigen::abs
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition Meta.h:82
Eigen::Dynamic
const int Dynamic
Definition Constants.h:25