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