JacobiSVD.h
Go to the documentation of this file.
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 //
6 // This Source Code Form is subject to the terms of the Mozilla
7 // Public License v. 2.0. If a copy of the MPL was not distributed
8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9 
10 #ifndef EIGEN_JACOBISVD_H
11 #define EIGEN_JACOBISVD_H
12 
13 namespace Eigen {
14 
15 namespace internal {
16 // forward declaration (needed by ICC)
17 // the empty body is required by MSVC
18 template<typename MatrixType, int QRPreconditioner,
21 
22 /*** QR preconditioners (R-SVD)
23  ***
24  *** Their role is to reduce the problem of computing the SVD to the case of a square matrix.
25  *** This approach, known as R-SVD, is an optimization for rectangular-enough matrices, and is a requirement for
26  *** JacobiSVD which by itself is only able to work on square matrices.
27  ***/
28 
30 
31 template<typename MatrixType, int QRPreconditioner, int Case>
33 {
34  enum { a = MatrixType::RowsAtCompileTime != Dynamic &&
35  MatrixType::ColsAtCompileTime != Dynamic &&
36  MatrixType::ColsAtCompileTime <= MatrixType::RowsAtCompileTime,
37  b = MatrixType::RowsAtCompileTime != Dynamic &&
38  MatrixType::ColsAtCompileTime != Dynamic &&
39  MatrixType::RowsAtCompileTime <= MatrixType::ColsAtCompileTime,
40  ret = !( (QRPreconditioner == NoQRPreconditioner) ||
41  (Case == PreconditionIfMoreColsThanRows && bool(a)) ||
42  (Case == PreconditionIfMoreRowsThanCols && bool(b)) )
43  };
44 };
45 
46 template<typename MatrixType, int QRPreconditioner, int Case,
49 
50 template<typename MatrixType, int QRPreconditioner, int Case>
51 class qr_preconditioner_impl<MatrixType, QRPreconditioner, Case, false>
52 {
53 public:
54  typedef typename MatrixType::Index Index;
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::Index Index;
69  typedef typename MatrixType::Scalar Scalar;
70  enum
71  {
72  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
73  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime
74  };
76 
78  {
79  if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols())
80  {
81  m_qr.~QRType();
82  ::new (&m_qr) QRType(svd.rows(), svd.cols());
83  }
84  if (svd.m_computeFullU) m_workspace.resize(svd.rows());
85  }
86 
88  {
89  if(matrix.rows() > matrix.cols())
90  {
91  m_qr.compute(matrix);
92  svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>();
93  if(svd.m_computeFullU) m_qr.matrixQ().evalTo(svd.m_matrixU, m_workspace);
94  if(svd.computeV()) svd.m_matrixV = m_qr.colsPermutation();
95  return true;
96  }
97  return false;
98  }
99 private:
101  QRType m_qr;
102  WorkspaceType m_workspace;
103 };
104 
105 template<typename MatrixType>
107 {
108 public:
109  typedef typename MatrixType::Index Index;
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  };
121 
123  {
124  if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols())
125  {
126  m_qr.~QRType();
127  ::new (&m_qr) QRType(svd.cols(), svd.rows());
128  }
129  m_adjoint.resize(svd.cols(), svd.rows());
130  if (svd.m_computeFullV) m_workspace.resize(svd.cols());
131  }
132 
134  {
135  if(matrix.cols() > matrix.rows())
136  {
137  m_adjoint = matrix.adjoint();
138  m_qr.compute(m_adjoint);
139  svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint();
140  if(svd.m_computeFullV) m_qr.matrixQ().evalTo(svd.m_matrixV, m_workspace);
141  if(svd.computeU()) svd.m_matrixU = m_qr.colsPermutation();
142  return true;
143  }
144  else return false;
145  }
146 private:
148  QRType m_qr;
151 };
152 
153 /*** preconditioner using ColPivHouseholderQR ***/
154 
155 template<typename MatrixType>
156 class qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
157 {
158 public:
159  typedef typename MatrixType::Index Index;
160 
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 
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 
190 private:
192  QRType m_qr;
194 };
195 
196 template<typename MatrixType>
198 {
199 public:
200  typedef typename MatrixType::Index Index;
201  typedef typename MatrixType::Scalar Scalar;
202  enum
203  {
204  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
205  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
206  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
207  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
208  Options = MatrixType::Options
209  };
210 
213 
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());
223  m_adjoint.resize(svd.cols(), svd.rows());
224  }
225 
227  {
228  if(matrix.cols() > matrix.rows())
229  {
230  m_adjoint = matrix.adjoint();
231  m_qr.compute(m_adjoint);
232 
233  svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint();
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 
246 private:
248  QRType m_qr;
251 };
252 
253 /*** preconditioner using HouseholderQR ***/
254 
255 template<typename MatrixType>
256 class qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
257 {
258 public:
259  typedef typename MatrixType::Index Index;
260 
262  {
263  if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols())
264  {
265  m_qr.~QRType();
266  ::new (&m_qr) QRType(svd.rows(), svd.cols());
267  }
268  if (svd.m_computeFullU) m_workspace.resize(svd.rows());
269  else if (svd.m_computeThinU) m_workspace.resize(svd.cols());
270  }
271 
272  bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix)
273  {
274  if(matrix.rows() > matrix.cols())
275  {
276  m_qr.compute(matrix);
277  svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>();
278  if(svd.m_computeFullU) m_qr.householderQ().evalTo(svd.m_matrixU, m_workspace);
279  else if(svd.m_computeThinU)
280  {
281  svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols());
282  m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_workspace);
283  }
284  if(svd.computeV()) svd.m_matrixV.setIdentity(matrix.cols(), matrix.cols());
285  return true;
286  }
287  return false;
288  }
289 private:
291  QRType m_qr;
293 };
294 
295 template<typename MatrixType>
297 {
298 public:
299  typedef typename MatrixType::Index Index;
300  typedef typename MatrixType::Scalar Scalar;
301  enum
302  {
303  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
304  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
305  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
306  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
307  Options = MatrixType::Options
308  };
309 
312 
314  {
315  if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols())
316  {
317  m_qr.~QRType();
318  ::new (&m_qr) QRType(svd.cols(), svd.rows());
319  }
320  if (svd.m_computeFullV) m_workspace.resize(svd.cols());
321  else if (svd.m_computeThinV) m_workspace.resize(svd.rows());
322  m_adjoint.resize(svd.cols(), svd.rows());
323  }
324 
325  bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix)
326  {
327  if(matrix.cols() > matrix.rows())
328  {
329  m_adjoint = matrix.adjoint();
330  m_qr.compute(m_adjoint);
331 
332  svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint();
333  if(svd.m_computeFullV) m_qr.householderQ().evalTo(svd.m_matrixV, m_workspace);
334  else if(svd.m_computeThinV)
335  {
336  svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows());
337  m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_workspace);
338  }
339  if(svd.computeU()) svd.m_matrixU.setIdentity(matrix.rows(), matrix.rows());
340  return true;
341  }
342  else return false;
343  }
344 
345 private:
347  QRType m_qr;
350 };
351 
352 /*** 2x2 SVD implementation
353  ***
354  *** JacobiSVD consists in performing a series of 2x2 SVD subproblems
355  ***/
356 
357 template<typename MatrixType, int QRPreconditioner>
358 struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, false>
359 {
361  typedef typename SVD::Index Index;
362  static void run(typename SVD::WorkMatrixType&, SVD&, Index, Index) {}
363 };
364 
365 template<typename MatrixType, int QRPreconditioner>
366 struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, true>
367 {
369  typedef typename MatrixType::Scalar Scalar;
370  typedef typename MatrixType::RealScalar RealScalar;
371  typedef typename SVD::Index Index;
372  static void run(typename SVD::WorkMatrixType& work_matrix, SVD& svd, Index p, Index q)
373  {
374  using std::sqrt;
375  Scalar z;
377  RealScalar n = sqrt(numext::abs2(work_matrix.coeff(p,p)) + numext::abs2(work_matrix.coeff(q,p)));
378  if(n==0)
379  {
380  z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q);
381  work_matrix.row(p) *= z;
382  if(svd.computeU()) svd.m_matrixU.col(p) *= conj(z);
383  z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q);
384  work_matrix.row(q) *= z;
385  if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z);
386  }
387  else
388  {
389  rot.c() = conj(work_matrix.coeff(p,p)) / n;
390  rot.s() = work_matrix.coeff(q,p) / n;
391  work_matrix.applyOnTheLeft(p,q,rot);
392  if(svd.computeU()) svd.m_matrixU.applyOnTheRight(p,q,rot.adjoint());
393  if(work_matrix.coeff(p,q) != Scalar(0))
394  {
395  Scalar z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q);
396  work_matrix.col(q) *= z;
397  if(svd.computeV()) svd.m_matrixV.col(q) *= z;
398  }
399  if(work_matrix.coeff(q,q) != Scalar(0))
400  {
401  z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q);
402  work_matrix.row(q) *= z;
403  if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z);
404  }
405  }
406  }
407 };
408 
409 template<typename MatrixType, typename RealScalar, typename Index>
410 void real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q,
413 {
414  using std::sqrt;
416  m << numext::real(matrix.coeff(p,p)), numext::real(matrix.coeff(p,q)),
417  numext::real(matrix.coeff(q,p)), numext::real(matrix.coeff(q,q));
419  RealScalar t = m.coeff(0,0) + m.coeff(1,1);
420  RealScalar d = m.coeff(1,0) - m.coeff(0,1);
421  if(t == RealScalar(0))
422  {
423  rot1.c() = RealScalar(0);
424  rot1.s() = d > RealScalar(0) ? RealScalar(1) : RealScalar(-1);
425  }
426  else
427  {
428  RealScalar u = d / t;
429  rot1.c() = RealScalar(1) / sqrt(RealScalar(1) + numext::abs2(u));
430  rot1.s() = rot1.c() * u;
431  }
432  m.applyOnTheLeft(0,1,rot1);
433  j_right->makeJacobi(m,0,1);
434  *j_left = rot1 * j_right->transpose();
435 }
436 
437 } // end namespace internal
438 
492 template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
493 {
494  public:
495 
496  typedef _MatrixType MatrixType;
497  typedef typename MatrixType::Scalar Scalar;
499  typedef typename MatrixType::Index Index;
500  enum {
501  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
502  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
503  DiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_DYNAMIC(RowsAtCompileTime,ColsAtCompileTime),
504  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
505  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
506  MaxDiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(MaxRowsAtCompileTime,MaxColsAtCompileTime),
507  MatrixOptions = MatrixType::Options
508  };
509 
510  typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime,
511  MatrixOptions, MaxRowsAtCompileTime, MaxRowsAtCompileTime>
513  typedef Matrix<Scalar, ColsAtCompileTime, ColsAtCompileTime,
514  MatrixOptions, MaxColsAtCompileTime, MaxColsAtCompileTime>
519  typedef Matrix<Scalar, DiagSizeAtCompileTime, DiagSizeAtCompileTime,
520  MatrixOptions, MaxDiagSizeAtCompileTime, MaxDiagSizeAtCompileTime>
522 
529  : m_isInitialized(false),
530  m_isAllocated(false),
531  m_computationOptions(0),
532  m_rows(-1), m_cols(-1)
533  {}
534 
535 
542  JacobiSVD(Index rows, Index cols, unsigned int computationOptions = 0)
543  : m_isInitialized(false),
544  m_isAllocated(false),
545  m_computationOptions(0),
546  m_rows(-1), m_cols(-1)
547  {
548  allocate(rows, cols, computationOptions);
549  }
550 
561  JacobiSVD(const MatrixType& matrix, unsigned int computationOptions = 0)
562  : m_isInitialized(false),
563  m_isAllocated(false),
564  m_computationOptions(0),
565  m_rows(-1), m_cols(-1)
566  {
567  compute(matrix, computationOptions);
568  }
569 
580  JacobiSVD& compute(const MatrixType& matrix, unsigned int computationOptions);
581 
588  JacobiSVD& compute(const MatrixType& matrix)
589  {
590  return compute(matrix, m_computationOptions);
591  }
592 
602  const MatrixUType& matrixU() const
603  {
604  eigen_assert(m_isInitialized && "JacobiSVD is not initialized.");
605  eigen_assert(computeU() && "This JacobiSVD decomposition didn't compute U. Did you ask for it?");
606  return m_matrixU;
607  }
608 
618  const MatrixVType& matrixV() const
619  {
620  eigen_assert(m_isInitialized && "JacobiSVD is not initialized.");
621  eigen_assert(computeV() && "This JacobiSVD decomposition didn't compute V. Did you ask for it?");
622  return m_matrixV;
623  }
624 
630  const SingularValuesType& singularValues() const
631  {
632  eigen_assert(m_isInitialized && "JacobiSVD is not initialized.");
633  return m_singularValues;
634  }
635 
637  inline bool computeU() const { return m_computeFullU || m_computeThinU; }
639  inline bool computeV() const { return m_computeFullV || m_computeThinV; }
640 
650  template<typename Rhs>
652  solve(const MatrixBase<Rhs>& b) const
653  {
654  eigen_assert(m_isInitialized && "JacobiSVD is not initialized.");
655  eigen_assert(computeU() && computeV() && "JacobiSVD::solve() requires both unitaries U and V to be computed (thin unitaries suffice).");
656  return internal::solve_retval<JacobiSVD, Rhs>(*this, b.derived());
657  }
658 
660  Index nonzeroSingularValues() const
661  {
662  eigen_assert(m_isInitialized && "JacobiSVD is not initialized.");
663  return m_nonzeroSingularValues;
664  }
665 
666  inline Index rows() const { return m_rows; }
667  inline Index cols() const { return m_cols; }
668 
669  private:
670  void allocate(Index rows, Index cols, unsigned int computationOptions);
671 
672  protected:
675  SingularValuesType m_singularValues;
677  bool m_isInitialized, m_isAllocated;
678  bool m_computeFullU, m_computeThinU;
679  bool m_computeFullV, m_computeThinV;
680  unsigned int m_computationOptions;
681  Index m_nonzeroSingularValues, m_rows, m_cols, m_diagSize;
682 
683  template<typename __MatrixType, int _QRPreconditioner, bool _IsComplex>
685  template<typename __MatrixType, int _QRPreconditioner, int _Case, bool _DoAnything>
687 
690 };
691 
692 template<typename MatrixType, int QRPreconditioner>
693 void JacobiSVD<MatrixType, QRPreconditioner>::allocate(Index rows, Index cols, unsigned int computationOptions)
694 {
695  eigen_assert(rows >= 0 && cols >= 0);
696 
697  if (m_isAllocated &&
698  rows == m_rows &&
699  cols == m_cols &&
700  computationOptions == m_computationOptions)
701  {
702  return;
703  }
704 
705  m_rows = rows;
706  m_cols = cols;
707  m_isInitialized = false;
708  m_isAllocated = true;
709  m_computationOptions = computationOptions;
710  m_computeFullU = (computationOptions & ComputeFullU) != 0;
711  m_computeThinU = (computationOptions & ComputeThinU) != 0;
712  m_computeFullV = (computationOptions & ComputeFullV) != 0;
713  m_computeThinV = (computationOptions & ComputeThinV) != 0;
714  eigen_assert(!(m_computeFullU && m_computeThinU) && "JacobiSVD: you can't ask for both full and thin U");
715  eigen_assert(!(m_computeFullV && m_computeThinV) && "JacobiSVD: you can't ask for both full and thin V");
716  eigen_assert(EIGEN_IMPLIES(m_computeThinU || m_computeThinV, MatrixType::ColsAtCompileTime==Dynamic) &&
717  "JacobiSVD: thin U and V are only available when your matrix has a dynamic number of columns.");
718  if (QRPreconditioner == FullPivHouseholderQRPreconditioner)
719  {
720  eigen_assert(!(m_computeThinU || m_computeThinV) &&
721  "JacobiSVD: can't compute thin U or thin V with the FullPivHouseholderQR preconditioner. "
722  "Use the ColPivHouseholderQR preconditioner instead.");
723  }
724  m_diagSize = (std::min)(m_rows, m_cols);
725  m_singularValues.resize(m_diagSize);
726  if(RowsAtCompileTime==Dynamic)
727  m_matrixU.resize(m_rows, m_computeFullU ? m_rows
728  : m_computeThinU ? m_diagSize
729  : 0);
730  if(ColsAtCompileTime==Dynamic)
731  m_matrixV.resize(m_cols, m_computeFullV ? m_cols
732  : m_computeThinV ? m_diagSize
733  : 0);
734  m_workMatrix.resize(m_diagSize, m_diagSize);
735 
736  if(m_cols>m_rows) m_qr_precond_morecols.allocate(*this);
737  if(m_rows>m_cols) m_qr_precond_morerows.allocate(*this);
738 }
739 
740 template<typename MatrixType, int QRPreconditioner>
742 JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsigned int computationOptions)
743 {
744  using std::abs;
745  allocate(matrix.rows(), matrix.cols(), computationOptions);
746 
747  // currently we stop when we reach precision 2*epsilon as the last bit of precision can require an unreasonable number of iterations,
748  // only worsening the precision of U and V as we accumulate more rotations
750 
751  // limit for very small denormal numbers to be considered zero in order to avoid infinite loops (see bug 286)
752  const RealScalar considerAsZero = RealScalar(2) * std::numeric_limits<RealScalar>::denorm_min();
753 
754  /*** step 1. The R-SVD step: we use a QR decomposition to reduce to the case of a square matrix */
755 
756  if(!m_qr_precond_morecols.run(*this, matrix) && !m_qr_precond_morerows.run(*this, matrix))
757  {
758  m_workMatrix = matrix.block(0,0,m_diagSize,m_diagSize);
759  if(m_computeFullU) m_matrixU.setIdentity(m_rows,m_rows);
760  if(m_computeThinU) m_matrixU.setIdentity(m_rows,m_diagSize);
761  if(m_computeFullV) m_matrixV.setIdentity(m_cols,m_cols);
762  if(m_computeThinV) m_matrixV.setIdentity(m_cols, m_diagSize);
763  }
764 
765  /*** step 2. The main Jacobi SVD iteration. ***/
766 
767  bool finished = false;
768  while(!finished)
769  {
770  finished = true;
771 
772  // do a sweep: for all index pairs (p,q), perform SVD of the corresponding 2x2 sub-matrix
773 
774  for(Index p = 1; p < m_diagSize; ++p)
775  {
776  for(Index q = 0; q < p; ++q)
777  {
778  // if this 2x2 sub-matrix is not diagonal already...
779  // notice that this comparison will evaluate to false if any NaN is involved, ensuring that NaN's don't
780  // keep us iterating forever. Similarly, small denormal numbers are considered zero.
781  using std::max;
782  RealScalar threshold = (max)(considerAsZero, precision * (max)(abs(m_workMatrix.coeff(p,p)),
783  abs(m_workMatrix.coeff(q,q))));
784  if((max)(abs(m_workMatrix.coeff(p,q)),abs(m_workMatrix.coeff(q,p))) > threshold)
785  {
786  finished = false;
787 
788  // perform SVD decomposition of 2x2 sub-matrix corresponding to indices p,q to make it diagonal
790  JacobiRotation<RealScalar> j_left, j_right;
791  internal::real_2x2_jacobi_svd(m_workMatrix, p, q, &j_left, &j_right);
792 
793  // accumulate resulting Jacobi rotations
794  m_workMatrix.applyOnTheLeft(p,q,j_left);
795  if(computeU()) m_matrixU.applyOnTheRight(p,q,j_left.transpose());
796 
797  m_workMatrix.applyOnTheRight(p,q,j_right);
798  if(computeV()) m_matrixV.applyOnTheRight(p,q,j_right);
799  }
800  }
801  }
802  }
803 
804  /*** step 3. The work matrix is now diagonal, so ensure it's positive so its diagonal entries are the singular values ***/
805 
806  for(Index i = 0; i < m_diagSize; ++i)
807  {
808  RealScalar a = abs(m_workMatrix.coeff(i,i));
809  m_singularValues.coeffRef(i) = a;
810  if(computeU() && (a!=RealScalar(0))) m_matrixU.col(i) *= m_workMatrix.coeff(i,i)/a;
811  }
812 
813  /*** step 4. Sort singular values in descending order and compute the number of nonzero singular values ***/
814 
815  m_nonzeroSingularValues = m_diagSize;
816  for(Index i = 0; i < m_diagSize; i++)
817  {
818  Index pos;
819  RealScalar maxRemainingSingularValue = m_singularValues.tail(m_diagSize-i).maxCoeff(&pos);
820  if(maxRemainingSingularValue == RealScalar(0))
821  {
822  m_nonzeroSingularValues = i;
823  break;
824  }
825  if(pos)
826  {
827  pos += i;
828  std::swap(m_singularValues.coeffRef(i), m_singularValues.coeffRef(pos));
829  if(computeU()) m_matrixU.col(pos).swap(m_matrixU.col(i));
830  if(computeV()) m_matrixV.col(pos).swap(m_matrixV.col(i));
831  }
832  }
833 
834  m_isInitialized = true;
835  return *this;
836 }
837 
838 namespace internal {
839 template<typename _MatrixType, int QRPreconditioner, typename Rhs>
840 struct solve_retval<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs>
841  : solve_retval_base<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs>
842 {
844  EIGEN_MAKE_SOLVE_HELPERS(JacobiSVDType,Rhs)
845 
846  template<typename Dest> void evalTo(Dest& dst) const
847  {
848  eigen_assert(rhs().rows() == dec().rows());
849 
850  // A = U S V^*
851  // So A^{-1} = V S^{-1} U^*
852 
854  Index nonzeroSingVals = dec().nonzeroSingularValues();
855 
856  tmp.noalias() = dec().matrixU().leftCols(nonzeroSingVals).adjoint() * rhs();
857  tmp = dec().singularValues().head(nonzeroSingVals).asDiagonal().inverse() * tmp;
858  dst = dec().matrixV().leftCols(nonzeroSingVals) * tmp;
859  }
860 };
861 } // end namespace internal
862 
870 template<typename Derived>
872 MatrixBase<Derived>::jacobiSvd(unsigned int computationOptions) const
873 {
874  return JacobiSVD<PlainObject>(*this, computationOptions);
875 }
876 
877 } // end namespace Eigen
878 
879 #endif // EIGEN_JACOBISVD_H
Matrix< Scalar, DiagSizeAtCompileTime, DiagSizeAtCompileTime, MatrixOptions, MaxDiagSizeAtCompileTime, MaxDiagSizeAtCompileTime > WorkMatrixType
Definition: JacobiSVD.h:521
Matrix< Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime > TransposeTypeWithSameStorageOrder
Definition: JacobiSVD.h:120
bool run(JacobiSVD< MatrixType, QRPreconditioner > &, const MatrixType &)
Definition: JacobiSVD.h:56
Scalar & c()
Definition: Jacobi.h:45
const AutoDiffScalar< DerType > & conj(const AutoDiffScalar< DerType > &x)
internal::plain_col_type< MatrixType >::type ColType
Definition: JacobiSVD.h:518
IntermediateState sqrt(const Expression &arg)
Matrix< Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime > TransposeTypeWithSameStorageOrder
Definition: JacobiSVD.h:311
void allocate(const JacobiSVD< MatrixType, QRPreconditioner > &)
Definition: JacobiSVD.h:55
bool m_isInitialized
Definition: JacobiSVD.h:677
Index cols() const
Definition: JacobiSVD.h:667
unsigned int m_computationOptions
Definition: JacobiSVD.h:680
static void run(typename SVD::WorkMatrixType &work_matrix, SVD &svd, Index p, Index q)
Definition: JacobiSVD.h:372
const internal::solve_retval< JacobiSVD, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: JacobiSVD.h:652
iterative scaling algorithm to equilibrate rows and column norms in matrices
Definition: matrix.hpp:471
Rotation given by a cosine-sine pair.
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:88
const SingularValuesType & singularValues() const
Definition: JacobiSVD.h:630
bool computeV() const
Definition: JacobiSVD.h:639
#define EIGEN_SIZE_MIN_PREFER_FIXED(a, b)
JacobiRotation transpose() const
Definition: Jacobi.h:59
EIGEN_STRONG_INLINE const CwiseUnaryOp< internal::scalar_abs2_op< Scalar >, const Derived > abs2() const
internal::plain_diag_type< MatrixType, RealScalar >::type SingularValuesType
Definition: JacobiSVD.h:516
MatrixType::Scalar Scalar
Definition: JacobiSVD.h:497
EIGEN_STRONG_INLINE const CwiseUnaryOp< internal::scalar_abs_op< Scalar >, const Derived > abs() const
RealReturnType real() const
bool run(JacobiSVD< MatrixType, ColPivHouseholderQRPreconditioner > &svd, const MatrixType &matrix)
Definition: JacobiSVD.h:172
void real_2x2_jacobi_svd(const MatrixType &matrix, Index p, Index q, JacobiRotation< RealScalar > *j_left, JacobiRotation< RealScalar > *j_right)
Definition: JacobiSVD.h:410
EIGEN_STRONG_INLINE const Scalar & coeff(Index rowId, Index colId) const
MatrixVType m_matrixV
Definition: JacobiSVD.h:674
Index nonzeroSingularValues() const
Definition: JacobiSVD.h:660
static void run(typename SVD::WorkMatrixType &, SVD &, Index, Index)
Definition: JacobiSVD.h:362
JacobiSVD()
Default Constructor.
Definition: JacobiSVD.h:528
MatrixType::Index Index
Definition: JacobiSVD.h:499
Provides a generic way to set and pass user-specified options.
Definition: options.hpp:65
bool makeJacobi(const MatrixBase< Derived > &, typename Derived::Index p, typename Derived::Index q)
Definition: Jacobi.h:126
internal::qr_preconditioner_impl< MatrixType, QRPreconditioner, internal::PreconditionIfMoreColsThanRows > m_qr_precond_morecols
Definition: JacobiSVD.h:688
Index rows() const
Definition: JacobiSVD.h:666
_MatrixType MatrixType
Definition: JacobiSVD.h:496
JacobiRotation adjoint() const
Definition: Jacobi.h:62
Matrix< Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime > TransposeTypeWithSameStorageOrder
Definition: JacobiSVD.h:212
Matrix< Scalar, RowsAtCompileTime, RowsAtCompileTime, MatrixOptions, MaxRowsAtCompileTime, MaxRowsAtCompileTime > MatrixUType
Definition: JacobiSVD.h:512
SingularValuesType m_singularValues
Definition: JacobiSVD.h:675
void rhs(const real_t *x, real_t *f)
bool run(JacobiSVD< MatrixType, HouseholderQRPreconditioner > &svd, const MatrixType &matrix)
Definition: JacobiSVD.h:272
bool run(JacobiSVD< MatrixType, FullPivHouseholderQRPreconditioner > &svd, const MatrixType &matrix)
Definition: JacobiSVD.h:87
bool run(JacobiSVD< MatrixType, HouseholderQRPreconditioner > &svd, const MatrixType &matrix)
Definition: JacobiSVD.h:325
#define EIGEN_SIZE_MIN_PREFER_DYNAMIC(a, b)
JacobiSVD & compute(const MatrixType &matrix)
Method performing the decomposition of given matrix using current options.
Definition: JacobiSVD.h:588
Matrix< Scalar, ColsAtCompileTime, ColsAtCompileTime, MatrixOptions, MaxColsAtCompileTime, MaxColsAtCompileTime > MatrixVType
Definition: JacobiSVD.h:515
void allocate(Index rows, Index cols, unsigned int computationOptions)
Definition: JacobiSVD.h:693
Two-sided Jacobi SVD decomposition of a rectangular matrix.
Scalar & s()
Definition: Jacobi.h:47
internal::plain_row_type< MatrixType >::type RowType
Definition: JacobiSVD.h:517
NumTraits< typename MatrixType::Scalar >::Real RealScalar
Definition: JacobiSVD.h:498
JacobiSVD & compute(const MatrixType &matrix, unsigned int computationOptions)
Method performing the decomposition of given matrix using custom options.
Definition: JacobiSVD.h:742
bool run(JacobiSVD< MatrixType, FullPivHouseholderQRPreconditioner > &svd, const MatrixType &matrix)
Definition: JacobiSVD.h:133
const MatrixUType & matrixU() const
Definition: JacobiSVD.h:602
#define EIGEN_MAKE_SOLVE_HELPERS(DecompositionType, Rhs)
Definition: Solve.h:61
#define EIGEN_IMPLIES(a, b)
MatrixUType m_matrixU
Definition: JacobiSVD.h:673
internal::qr_preconditioner_impl< MatrixType, QRPreconditioner, internal::PreconditionIfMoreRowsThanCols > m_qr_precond_morerows
Definition: JacobiSVD.h:689
#define eigen_assert(x)
JacobiSVD(Index rows, Index cols, unsigned int computationOptions=0)
Default Constructor with memory preallocation.
Definition: JacobiSVD.h:542
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
bool run(JacobiSVD< MatrixType, ColPivHouseholderQRPreconditioner > &svd, const MatrixType &matrix)
Definition: JacobiSVD.h:226
WorkMatrixType m_workMatrix
Definition: JacobiSVD.h:676
const MatrixVType & matrixV() const
Definition: JacobiSVD.h:618
JacobiSVD(const MatrixType &matrix, unsigned int computationOptions=0)
Constructor performing the decomposition of given matrix.
Definition: JacobiSVD.h:561
bool computeU() const
Definition: JacobiSVD.h:637
JacobiSVD< PlainObject > jacobiSvd(unsigned int computationOptions=0) const
Definition: JacobiSVD.h:872


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Mon Jun 10 2019 12:34:46