Eigen  3.2.91
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 namespace Eigen {
15 
16 namespace internal {
17 // forward declaration (needed by ICC)
18 // the empty body is required by MSVC
19 template<typename MatrixType, int QRPreconditioner,
20  bool IsComplex = NumTraits<typename MatrixType::Scalar>::IsComplex>
21 struct svd_precondition_2x2_block_to_be_real {};
22 
23 /*** QR preconditioners (R-SVD)
24  ***
25  *** Their role is to reduce the problem of computing the SVD to the case of a square matrix.
26  *** This approach, known as R-SVD, is an optimization for rectangular-enough matrices, and is a requirement for
27  *** JacobiSVD which by itself is only able to work on square matrices.
28  ***/
29 
30 enum { PreconditionIfMoreColsThanRows, PreconditionIfMoreRowsThanCols };
31 
32 template<typename MatrixType, int QRPreconditioner, int Case>
33 struct qr_preconditioner_should_do_anything
34 {
35  enum { a = MatrixType::RowsAtCompileTime != Dynamic &&
36  MatrixType::ColsAtCompileTime != Dynamic &&
37  MatrixType::ColsAtCompileTime <= MatrixType::RowsAtCompileTime,
38  b = MatrixType::RowsAtCompileTime != Dynamic &&
39  MatrixType::ColsAtCompileTime != Dynamic &&
40  MatrixType::RowsAtCompileTime <= MatrixType::ColsAtCompileTime,
41  ret = !( (QRPreconditioner == NoQRPreconditioner) ||
42  (Case == PreconditionIfMoreColsThanRows && bool(a)) ||
43  (Case == PreconditionIfMoreRowsThanCols && bool(b)) )
44  };
45 };
46 
47 template<typename MatrixType, int QRPreconditioner, int Case,
48  bool DoAnything = qr_preconditioner_should_do_anything<MatrixType, QRPreconditioner, Case>::ret
49 > struct qr_preconditioner_impl {};
50 
51 template<typename MatrixType, int QRPreconditioner, int Case>
52 class qr_preconditioner_impl<MatrixType, QRPreconditioner, Case, false>
53 {
54 public:
55  void allocate(const JacobiSVD<MatrixType, QRPreconditioner>&) {}
56  bool run(JacobiSVD<MatrixType, QRPreconditioner>&, const MatrixType&)
57  {
58  return false;
59  }
60 };
61 
62 /*** preconditioner using FullPivHouseholderQR ***/
63 
64 template<typename MatrixType>
65 class qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
66 {
67 public:
68  typedef typename MatrixType::Scalar Scalar;
69  enum
70  {
71  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
72  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime
73  };
74  typedef Matrix<Scalar, 1, RowsAtCompileTime, RowMajor, 1, MaxRowsAtCompileTime> WorkspaceType;
75 
76  void allocate(const JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd)
77  {
78  if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols())
79  {
80  m_qr.~QRType();
81  ::new (&m_qr) QRType(svd.rows(), svd.cols());
82  }
83  if (svd.m_computeFullU) m_workspace.resize(svd.rows());
84  }
85 
86  bool run(JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
87  {
88  if(matrix.rows() > matrix.cols())
89  {
90  m_qr.compute(matrix);
91  svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>();
92  if(svd.m_computeFullU) m_qr.matrixQ().evalTo(svd.m_matrixU, m_workspace);
93  if(svd.computeV()) svd.m_matrixV = m_qr.colsPermutation();
94  return true;
95  }
96  return false;
97  }
98 private:
99  typedef FullPivHouseholderQR<MatrixType> QRType;
100  QRType m_qr;
101  WorkspaceType m_workspace;
102 };
103 
104 template<typename MatrixType>
105 class qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true>
106 {
107 public:
108  typedef typename MatrixType::Scalar Scalar;
109  enum
110  {
111  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
112  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
113  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
114  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
115  Options = MatrixType::Options
116  };
117  typedef Matrix<Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime>
118  TransposeTypeWithSameStorageOrder;
119 
120  void allocate(const JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd)
121  {
122  if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols())
123  {
124  m_qr.~QRType();
125  ::new (&m_qr) QRType(svd.cols(), svd.rows());
126  }
127  m_adjoint.resize(svd.cols(), svd.rows());
128  if (svd.m_computeFullV) m_workspace.resize(svd.cols());
129  }
130 
131  bool run(JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
132  {
133  if(matrix.cols() > matrix.rows())
134  {
135  m_adjoint = matrix.adjoint();
136  m_qr.compute(m_adjoint);
137  svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint();
138  if(svd.m_computeFullV) m_qr.matrixQ().evalTo(svd.m_matrixV, m_workspace);
139  if(svd.computeU()) svd.m_matrixU = m_qr.colsPermutation();
140  return true;
141  }
142  else return false;
143  }
144 private:
145  typedef FullPivHouseholderQR<TransposeTypeWithSameStorageOrder> QRType;
146  QRType m_qr;
147  TransposeTypeWithSameStorageOrder m_adjoint;
148  typename internal::plain_row_type<MatrixType>::type m_workspace;
149 };
150 
151 /*** preconditioner using ColPivHouseholderQR ***/
152 
153 template<typename MatrixType>
154 class qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
155 {
156 public:
157  void allocate(const JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd)
158  {
159  if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols())
160  {
161  m_qr.~QRType();
162  ::new (&m_qr) QRType(svd.rows(), svd.cols());
163  }
164  if (svd.m_computeFullU) m_workspace.resize(svd.rows());
165  else if (svd.m_computeThinU) m_workspace.resize(svd.cols());
166  }
167 
168  bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
169  {
170  if(matrix.rows() > matrix.cols())
171  {
172  m_qr.compute(matrix);
173  svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>();
174  if(svd.m_computeFullU) m_qr.householderQ().evalTo(svd.m_matrixU, m_workspace);
175  else if(svd.m_computeThinU)
176  {
177  svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols());
178  m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_workspace);
179  }
180  if(svd.computeV()) svd.m_matrixV = m_qr.colsPermutation();
181  return true;
182  }
183  return false;
184  }
185 
186 private:
187  typedef ColPivHouseholderQR<MatrixType> QRType;
188  QRType m_qr;
189  typename internal::plain_col_type<MatrixType>::type m_workspace;
190 };
191 
192 template<typename MatrixType>
193 class qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true>
194 {
195 public:
196  typedef typename MatrixType::Scalar Scalar;
197  enum
198  {
199  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
200  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
201  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
202  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
203  Options = MatrixType::Options
204  };
205 
206  typedef Matrix<Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime>
207  TransposeTypeWithSameStorageOrder;
208 
209  void allocate(const JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd)
210  {
211  if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols())
212  {
213  m_qr.~QRType();
214  ::new (&m_qr) QRType(svd.cols(), svd.rows());
215  }
216  if (svd.m_computeFullV) m_workspace.resize(svd.cols());
217  else if (svd.m_computeThinV) m_workspace.resize(svd.rows());
218  m_adjoint.resize(svd.cols(), svd.rows());
219  }
220 
221  bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
222  {
223  if(matrix.cols() > matrix.rows())
224  {
225  m_adjoint = matrix.adjoint();
226  m_qr.compute(m_adjoint);
227 
228  svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint();
229  if(svd.m_computeFullV) m_qr.householderQ().evalTo(svd.m_matrixV, m_workspace);
230  else if(svd.m_computeThinV)
231  {
232  svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows());
233  m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_workspace);
234  }
235  if(svd.computeU()) svd.m_matrixU = m_qr.colsPermutation();
236  return true;
237  }
238  else return false;
239  }
240 
241 private:
242  typedef ColPivHouseholderQR<TransposeTypeWithSameStorageOrder> QRType;
243  QRType m_qr;
244  TransposeTypeWithSameStorageOrder m_adjoint;
245  typename internal::plain_row_type<MatrixType>::type m_workspace;
246 };
247 
248 /*** preconditioner using HouseholderQR ***/
249 
250 template<typename MatrixType>
251 class qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
252 {
253 public:
254  void allocate(const JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd)
255  {
256  if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols())
257  {
258  m_qr.~QRType();
259  ::new (&m_qr) QRType(svd.rows(), svd.cols());
260  }
261  if (svd.m_computeFullU) m_workspace.resize(svd.rows());
262  else if (svd.m_computeThinU) m_workspace.resize(svd.cols());
263  }
264 
265  bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix)
266  {
267  if(matrix.rows() > matrix.cols())
268  {
269  m_qr.compute(matrix);
270  svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>();
271  if(svd.m_computeFullU) m_qr.householderQ().evalTo(svd.m_matrixU, m_workspace);
272  else if(svd.m_computeThinU)
273  {
274  svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols());
275  m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_workspace);
276  }
277  if(svd.computeV()) svd.m_matrixV.setIdentity(matrix.cols(), matrix.cols());
278  return true;
279  }
280  return false;
281  }
282 private:
283  typedef HouseholderQR<MatrixType> QRType;
284  QRType m_qr;
285  typename internal::plain_col_type<MatrixType>::type m_workspace;
286 };
287 
288 template<typename MatrixType>
289 class qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true>
290 {
291 public:
292  typedef typename MatrixType::Scalar Scalar;
293  enum
294  {
295  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
296  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
297  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
298  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
299  Options = MatrixType::Options
300  };
301 
302  typedef Matrix<Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime>
303  TransposeTypeWithSameStorageOrder;
304 
305  void allocate(const JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd)
306  {
307  if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols())
308  {
309  m_qr.~QRType();
310  ::new (&m_qr) QRType(svd.cols(), svd.rows());
311  }
312  if (svd.m_computeFullV) m_workspace.resize(svd.cols());
313  else if (svd.m_computeThinV) m_workspace.resize(svd.rows());
314  m_adjoint.resize(svd.cols(), svd.rows());
315  }
316 
317  bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix)
318  {
319  if(matrix.cols() > matrix.rows())
320  {
321  m_adjoint = matrix.adjoint();
322  m_qr.compute(m_adjoint);
323 
324  svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint();
325  if(svd.m_computeFullV) m_qr.householderQ().evalTo(svd.m_matrixV, m_workspace);
326  else if(svd.m_computeThinV)
327  {
328  svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows());
329  m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_workspace);
330  }
331  if(svd.computeU()) svd.m_matrixU.setIdentity(matrix.rows(), matrix.rows());
332  return true;
333  }
334  else return false;
335  }
336 
337 private:
338  typedef HouseholderQR<TransposeTypeWithSameStorageOrder> QRType;
339  QRType m_qr;
340  TransposeTypeWithSameStorageOrder m_adjoint;
341  typename internal::plain_row_type<MatrixType>::type m_workspace;
342 };
343 
344 /*** 2x2 SVD implementation
345  ***
346  *** JacobiSVD consists in performing a series of 2x2 SVD subproblems
347  ***/
348 
349 template<typename MatrixType, int QRPreconditioner>
350 struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, false>
351 {
352  typedef JacobiSVD<MatrixType, QRPreconditioner> SVD;
353  static void run(typename SVD::WorkMatrixType&, SVD&, Index, Index) {}
354 };
355 
356 template<typename MatrixType, int QRPreconditioner>
357 struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, true>
358 {
359  typedef JacobiSVD<MatrixType, QRPreconditioner> SVD;
360  typedef typename MatrixType::Scalar Scalar;
361  typedef typename MatrixType::RealScalar RealScalar;
362  static void run(typename SVD::WorkMatrixType& work_matrix, SVD& svd, Index p, Index q)
363  {
364  using std::sqrt;
365  Scalar z;
366  JacobiRotation<Scalar> rot;
367  RealScalar n = sqrt(numext::abs2(work_matrix.coeff(p,p)) + numext::abs2(work_matrix.coeff(q,p)));
368 
369  if(n==0)
370  {
371  z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q);
372  work_matrix.row(p) *= z;
373  if(svd.computeU()) svd.m_matrixU.col(p) *= conj(z);
374  if(work_matrix.coeff(q,q)!=Scalar(0))
375  {
376  z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q);
377  work_matrix.row(q) *= z;
378  if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z);
379  }
380  // otherwise the second row is already zero, so we have nothing to do.
381  }
382  else
383  {
384  rot.c() = conj(work_matrix.coeff(p,p)) / n;
385  rot.s() = work_matrix.coeff(q,p) / n;
386  work_matrix.applyOnTheLeft(p,q,rot);
387  if(svd.computeU()) svd.m_matrixU.applyOnTheRight(p,q,rot.adjoint());
388  if(work_matrix.coeff(p,q) != Scalar(0))
389  {
390  z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q);
391  work_matrix.col(q) *= z;
392  if(svd.computeV()) svd.m_matrixV.col(q) *= z;
393  }
394  if(work_matrix.coeff(q,q) != Scalar(0))
395  {
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  }
401  }
402 };
403 
404 template<typename MatrixType, typename RealScalar, typename Index>
405 void real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q,
406  JacobiRotation<RealScalar> *j_left,
407  JacobiRotation<RealScalar> *j_right)
408 {
409  using std::sqrt;
410  using std::abs;
411  Matrix<RealScalar,2,2> m;
412  m << numext::real(matrix.coeff(p,p)), numext::real(matrix.coeff(p,q)),
413  numext::real(matrix.coeff(q,p)), numext::real(matrix.coeff(q,q));
414  JacobiRotation<RealScalar> rot1;
415  RealScalar t = m.coeff(0,0) + m.coeff(1,1);
416  RealScalar d = m.coeff(1,0) - m.coeff(0,1);
417 
418  if(d == RealScalar(0))
419  {
420  rot1.s() = RealScalar(0);
421  rot1.c() = RealScalar(1);
422  }
423  else
424  {
425  // If d!=0, then t/d cannot overflow because the magnitude of the
426  // entries forming d are not too small compared to the ones forming t.
427  RealScalar u = t / d;
428  RealScalar tmp = sqrt(RealScalar(1) + numext::abs2(u));
429  rot1.s() = RealScalar(1) / tmp;
430  rot1.c() = u / tmp;
431  }
432  m.applyOnTheLeft(0,1,rot1);
433  j_right->makeJacobi(m,0,1);
434  *j_left = rot1 * j_right->transpose();
435 }
436 
437 template<typename _MatrixType, int QRPreconditioner>
438 struct traits<JacobiSVD<_MatrixType,QRPreconditioner> >
439 {
440  typedef _MatrixType MatrixType;
441 };
442 
443 } // end namespace internal
444 
498 template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
499  : public SVDBase<JacobiSVD<_MatrixType,QRPreconditioner> >
500 {
501  typedef SVDBase<JacobiSVD> Base;
502  public:
503 
504  typedef _MatrixType MatrixType;
505  typedef typename MatrixType::Scalar Scalar;
506  typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
507  enum {
508  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
509  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
510  DiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_DYNAMIC(RowsAtCompileTime,ColsAtCompileTime),
511  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
512  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
513  MaxDiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(MaxRowsAtCompileTime,MaxColsAtCompileTime),
514  MatrixOptions = MatrixType::Options
515  };
516 
517  typedef typename Base::MatrixUType MatrixUType;
518  typedef typename Base::MatrixVType MatrixVType;
519  typedef typename Base::SingularValuesType SingularValuesType;
520 
521  typedef typename internal::plain_row_type<MatrixType>::type RowType;
522  typedef typename internal::plain_col_type<MatrixType>::type ColType;
523  typedef Matrix<Scalar, DiagSizeAtCompileTime, DiagSizeAtCompileTime,
524  MatrixOptions, MaxDiagSizeAtCompileTime, MaxDiagSizeAtCompileTime>
525  WorkMatrixType;
526 
533  {}
534 
535 
542  explicit JacobiSVD(Index rows, Index cols, unsigned int computationOptions = 0)
543  {
544  allocate(rows, cols, computationOptions);
545  }
546 
557  explicit JacobiSVD(const MatrixType& matrix, unsigned int computationOptions = 0)
558  {
559  compute(matrix, computationOptions);
560  }
561 
572  JacobiSVD& compute(const MatrixType& matrix, unsigned int computationOptions);
573 
580  JacobiSVD& compute(const MatrixType& matrix)
581  {
582  return compute(matrix, m_computationOptions);
583  }
584 
585  using Base::computeU;
586  using Base::computeV;
587  using Base::rows;
588  using Base::cols;
589  using Base::rank;
590 
591  private:
592  void allocate(Index rows, Index cols, unsigned int computationOptions);
593 
594  protected:
595  using Base::m_matrixU;
596  using Base::m_matrixV;
597  using Base::m_singularValues;
598  using Base::m_isInitialized;
599  using Base::m_isAllocated;
600  using Base::m_usePrescribedThreshold;
601  using Base::m_computeFullU;
602  using Base::m_computeThinU;
603  using Base::m_computeFullV;
604  using Base::m_computeThinV;
605  using Base::m_computationOptions;
606  using Base::m_nonzeroSingularValues;
607  using Base::m_rows;
608  using Base::m_cols;
609  using Base::m_diagSize;
610  using Base::m_prescribedThreshold;
611  WorkMatrixType m_workMatrix;
612 
613  template<typename __MatrixType, int _QRPreconditioner, bool _IsComplex>
614  friend struct internal::svd_precondition_2x2_block_to_be_real;
615  template<typename __MatrixType, int _QRPreconditioner, int _Case, bool _DoAnything>
616  friend struct internal::qr_preconditioner_impl;
617 
618  internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreColsThanRows> m_qr_precond_morecols;
619  internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreRowsThanCols> m_qr_precond_morerows;
620  MatrixType m_scaledMatrix;
621 };
622 
623 template<typename MatrixType, int QRPreconditioner>
624 void JacobiSVD<MatrixType, QRPreconditioner>::allocate(Index rows, Index cols, unsigned int computationOptions)
625 {
626  eigen_assert(rows >= 0 && cols >= 0);
627 
628  if (m_isAllocated &&
629  rows == m_rows &&
630  cols == m_cols &&
631  computationOptions == m_computationOptions)
632  {
633  return;
634  }
635 
636  m_rows = rows;
637  m_cols = cols;
638  m_isInitialized = false;
639  m_isAllocated = true;
640  m_computationOptions = computationOptions;
641  m_computeFullU = (computationOptions & ComputeFullU) != 0;
642  m_computeThinU = (computationOptions & ComputeThinU) != 0;
643  m_computeFullV = (computationOptions & ComputeFullV) != 0;
644  m_computeThinV = (computationOptions & ComputeThinV) != 0;
645  eigen_assert(!(m_computeFullU && m_computeThinU) && "JacobiSVD: you can't ask for both full and thin U");
646  eigen_assert(!(m_computeFullV && m_computeThinV) && "JacobiSVD: you can't ask for both full and thin V");
647  eigen_assert(EIGEN_IMPLIES(m_computeThinU || m_computeThinV, MatrixType::ColsAtCompileTime==Dynamic) &&
648  "JacobiSVD: thin U and V are only available when your matrix has a dynamic number of columns.");
649  if (QRPreconditioner == FullPivHouseholderQRPreconditioner)
650  {
651  eigen_assert(!(m_computeThinU || m_computeThinV) &&
652  "JacobiSVD: can't compute thin U or thin V with the FullPivHouseholderQR preconditioner. "
653  "Use the ColPivHouseholderQR preconditioner instead.");
654  }
655  m_diagSize = (std::min)(m_rows, m_cols);
656  m_singularValues.resize(m_diagSize);
657  if(RowsAtCompileTime==Dynamic)
658  m_matrixU.resize(m_rows, m_computeFullU ? m_rows
659  : m_computeThinU ? m_diagSize
660  : 0);
661  if(ColsAtCompileTime==Dynamic)
662  m_matrixV.resize(m_cols, m_computeFullV ? m_cols
663  : m_computeThinV ? m_diagSize
664  : 0);
665  m_workMatrix.resize(m_diagSize, m_diagSize);
666 
667  if(m_cols>m_rows) m_qr_precond_morecols.allocate(*this);
668  if(m_rows>m_cols) m_qr_precond_morerows.allocate(*this);
669  if(m_cols!=m_cols) m_scaledMatrix.resize(rows,cols);
670 }
671 
672 template<typename MatrixType, int QRPreconditioner>
673 JacobiSVD<MatrixType, QRPreconditioner>&
674 JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsigned int computationOptions)
675 {
676  using std::abs;
677  allocate(matrix.rows(), matrix.cols(), computationOptions);
678 
679  // currently we stop when we reach precision 2*epsilon as the last bit of precision can require an unreasonable number of iterations,
680  // only worsening the precision of U and V as we accumulate more rotations
681  const RealScalar precision = RealScalar(2) * NumTraits<Scalar>::epsilon();
682 
683  // limit for very small denormal numbers to be considered zero in order to avoid infinite loops (see bug 286)
684  // FIXME What about considerering any denormal numbers as zero, using:
685  // const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)();
686  const RealScalar considerAsZero = RealScalar(2) * std::numeric_limits<RealScalar>::denorm_min();
687 
688  // Scaling factor to reduce over/under-flows
689  RealScalar scale = matrix.cwiseAbs().maxCoeff();
690  if(scale==RealScalar(0)) scale = RealScalar(1);
691 
692  /*** step 1. The R-SVD step: we use a QR decomposition to reduce to the case of a square matrix */
693 
694  if(m_rows!=m_cols)
695  {
696  m_scaledMatrix = matrix / scale;
697  m_qr_precond_morecols.run(*this, m_scaledMatrix);
698  m_qr_precond_morerows.run(*this, m_scaledMatrix);
699  }
700  else
701  {
702  m_workMatrix = matrix.block(0,0,m_diagSize,m_diagSize) / scale;
703  if(m_computeFullU) m_matrixU.setIdentity(m_rows,m_rows);
704  if(m_computeThinU) m_matrixU.setIdentity(m_rows,m_diagSize);
705  if(m_computeFullV) m_matrixV.setIdentity(m_cols,m_cols);
706  if(m_computeThinV) m_matrixV.setIdentity(m_cols, m_diagSize);
707  }
708 
709  /*** step 2. The main Jacobi SVD iteration. ***/
710 
711  bool finished = false;
712  while(!finished)
713  {
714  finished = true;
715 
716  // do a sweep: for all index pairs (p,q), perform SVD of the corresponding 2x2 sub-matrix
717 
718  for(Index p = 1; p < m_diagSize; ++p)
719  {
720  for(Index q = 0; q < p; ++q)
721  {
722  // if this 2x2 sub-matrix is not diagonal already...
723  // notice that this comparison will evaluate to false if any NaN is involved, ensuring that NaN's don't
724  // keep us iterating forever. Similarly, small denormal numbers are considered zero.
725  RealScalar threshold = numext::maxi<RealScalar>(considerAsZero,
726  precision * numext::maxi<RealScalar>(abs(m_workMatrix.coeff(p,p)),
727  abs(m_workMatrix.coeff(q,q))));
728  // We compare both values to threshold instead of calling max to be robust to NaN (See bug 791)
729  if(abs(m_workMatrix.coeff(p,q))>threshold || abs(m_workMatrix.coeff(q,p)) > threshold)
730  {
731  finished = false;
732 
733  // perform SVD decomposition of 2x2 sub-matrix corresponding to indices p,q to make it diagonal
734  internal::svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner>::run(m_workMatrix, *this, p, q);
735  JacobiRotation<RealScalar> j_left, j_right;
736  internal::real_2x2_jacobi_svd(m_workMatrix, p, q, &j_left, &j_right);
737 
738  // accumulate resulting Jacobi rotations
739  m_workMatrix.applyOnTheLeft(p,q,j_left);
740  if(computeU()) m_matrixU.applyOnTheRight(p,q,j_left.transpose());
741 
742  m_workMatrix.applyOnTheRight(p,q,j_right);
743  if(computeV()) m_matrixV.applyOnTheRight(p,q,j_right);
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  RealScalar a = abs(m_workMatrix.coeff(i,i));
754  m_singularValues.coeffRef(i) = a;
755  if(computeU() && (a!=RealScalar(0))) m_matrixU.col(i) *= m_workMatrix.coeff(i,i)/a;
756  }
757 
758  m_singularValues *= scale;
759 
760  /*** step 4. Sort singular values in descending order and compute the number of nonzero singular values ***/
761 
762  m_nonzeroSingularValues = m_diagSize;
763  for(Index i = 0; i < m_diagSize; i++)
764  {
765  Index pos;
766  RealScalar maxRemainingSingularValue = m_singularValues.tail(m_diagSize-i).maxCoeff(&pos);
767  if(maxRemainingSingularValue == RealScalar(0))
768  {
769  m_nonzeroSingularValues = i;
770  break;
771  }
772  if(pos)
773  {
774  pos += i;
775  std::swap(m_singularValues.coeffRef(i), m_singularValues.coeffRef(pos));
776  if(computeU()) m_matrixU.col(pos).swap(m_matrixU.col(i));
777  if(computeV()) m_matrixV.col(pos).swap(m_matrixV.col(i));
778  }
779  }
780 
781  m_isInitialized = true;
782  return *this;
783 }
784 
785 #ifndef __CUDACC__
786 
793 template<typename Derived>
795 MatrixBase<Derived>::jacobiSvd(unsigned int computationOptions) const
796 {
797  return JacobiSVD<PlainObject>(*this, computationOptions);
798 }
799 #endif // __CUDACC__
800 
801 } // end namespace Eigen
802 
803 #endif // EIGEN_JACOBISVD_H
Definition: Constants.h:375
Definition: Constants.h:407
Definition: Constants.h:381
Definition: LDLT.h:16
Rotation given by a cosine-sine pair.
Definition: ForwardDeclarations.h:259
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:107
Definition: Constants.h:409
JacobiRotation transpose() const
Definition: Jacobi.h:59
Definition: Constants.h:411
JacobiSVD()
Default Constructor.
Definition: JacobiSVD.h:532
Definition: Eigen_Colamd.h:54
JacobiSVD & compute(const MatrixType &matrix)
Method performing the decomposition of given matrix using current options.
Definition: JacobiSVD.h:580
Two-sided Jacobi SVD decomposition of a rectangular matrix.
Definition: ForwardDeclarations.h:254
JacobiSVD & compute(const MatrixType &matrix, unsigned int computationOptions)
Method performing the decomposition of given matrix using custom options.
Definition: JacobiSVD.h:674
Definition: Constants.h:379
Definition: Constants.h:377
JacobiSVD(Index rows, Index cols, unsigned int computationOptions=0)
Default Constructor with memory preallocation.
Definition: JacobiSVD.h:542
JacobiSVD(const MatrixType &matrix, unsigned int computationOptions=0)
Constructor performing the decomposition of given matrix.
Definition: JacobiSVD.h:557
JacobiSVD< PlainObject > jacobiSvd(unsigned int computationOptions=0) const
Definition: JacobiSVD.h:795