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 // 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,
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 
31 
32 template<typename MatrixType, int QRPreconditioner, int Case>
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,
50 
51 template<typename MatrixType, int QRPreconditioner, int Case>
52 class qr_preconditioner_impl<MatrixType, QRPreconditioner, Case, false>
53 {
54 public:
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  };
75 
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 
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:
100  QRType m_qr;
101  WorkspaceType m_workspace;
102 };
103 
104 template<typename MatrixType>
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 
118  typedef typename internal::make_proper_matrix_type<
119  Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime
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;
149  TransposeTypeWithSameStorageOrder m_adjoint;
151 };
152 
153 /*** preconditioner using ColPivHouseholderQR ***/
154 
155 template<typename MatrixType>
156 class qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
157 {
158 public:
160  {
161  if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols())
162  {
163  m_qr.~QRType();
164  ::new (&m_qr) QRType(svd.rows(), svd.cols());
165  }
166  if (svd.m_computeFullU) m_workspace.resize(svd.rows());
167  else if (svd.m_computeThinU) m_workspace.resize(svd.cols());
168  }
169 
171  {
172  if(matrix.rows() > matrix.cols())
173  {
174  m_qr.compute(matrix);
175  svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>();
176  if(svd.m_computeFullU) m_qr.householderQ().evalTo(svd.m_matrixU, m_workspace);
177  else if(svd.m_computeThinU)
178  {
179  svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols());
180  m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_workspace);
181  }
182  if(svd.computeV()) svd.m_matrixV = m_qr.colsPermutation();
183  return true;
184  }
185  return false;
186  }
187 
188 private:
190  QRType m_qr;
192 };
193 
194 template<typename MatrixType>
196 {
197 public:
198  typedef typename MatrixType::Scalar Scalar;
199  enum
200  {
201  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
202  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
203  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
204  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
205  Options = MatrixType::Options
206  };
207 
208  typedef typename internal::make_proper_matrix_type<
209  Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime
211 
213  {
214  if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols())
215  {
216  m_qr.~QRType();
217  ::new (&m_qr) QRType(svd.cols(), svd.rows());
218  }
219  if (svd.m_computeFullV) m_workspace.resize(svd.cols());
220  else if (svd.m_computeThinV) m_workspace.resize(svd.rows());
221  m_adjoint.resize(svd.cols(), svd.rows());
222  }
223 
225  {
226  if(matrix.cols() > matrix.rows())
227  {
228  m_adjoint = matrix.adjoint();
229  m_qr.compute(m_adjoint);
230 
231  svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint();
232  if(svd.m_computeFullV) m_qr.householderQ().evalTo(svd.m_matrixV, m_workspace);
233  else if(svd.m_computeThinV)
234  {
235  svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows());
236  m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_workspace);
237  }
238  if(svd.computeU()) svd.m_matrixU = m_qr.colsPermutation();
239  return true;
240  }
241  else return false;
242  }
243 
244 private:
246  QRType m_qr;
247  TransposeTypeWithSameStorageOrder m_adjoint;
249 };
250 
251 /*** preconditioner using HouseholderQR ***/
252 
253 template<typename MatrixType>
254 class qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
255 {
256 public:
258  {
259  if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols())
260  {
261  m_qr.~QRType();
262  ::new (&m_qr) QRType(svd.rows(), svd.cols());
263  }
264  if (svd.m_computeFullU) m_workspace.resize(svd.rows());
265  else if (svd.m_computeThinU) m_workspace.resize(svd.cols());
266  }
267 
269  {
270  if(matrix.rows() > matrix.cols())
271  {
272  m_qr.compute(matrix);
273  svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>();
274  if(svd.m_computeFullU) m_qr.householderQ().evalTo(svd.m_matrixU, m_workspace);
275  else if(svd.m_computeThinU)
276  {
277  svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols());
278  m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_workspace);
279  }
280  if(svd.computeV()) svd.m_matrixV.setIdentity(matrix.cols(), matrix.cols());
281  return true;
282  }
283  return false;
284  }
285 private:
287  QRType m_qr;
289 };
290 
291 template<typename MatrixType>
293 {
294 public:
295  typedef typename MatrixType::Scalar Scalar;
296  enum
297  {
298  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
299  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
300  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
301  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
302  Options = MatrixType::Options
303  };
304 
305  typedef typename internal::make_proper_matrix_type<
306  Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime
308 
310  {
311  if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols())
312  {
313  m_qr.~QRType();
314  ::new (&m_qr) QRType(svd.cols(), svd.rows());
315  }
316  if (svd.m_computeFullV) m_workspace.resize(svd.cols());
317  else if (svd.m_computeThinV) m_workspace.resize(svd.rows());
318  m_adjoint.resize(svd.cols(), svd.rows());
319  }
320 
322  {
323  if(matrix.cols() > matrix.rows())
324  {
325  m_adjoint = matrix.adjoint();
326  m_qr.compute(m_adjoint);
327 
328  svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint();
329  if(svd.m_computeFullV) m_qr.householderQ().evalTo(svd.m_matrixV, m_workspace);
330  else if(svd.m_computeThinV)
331  {
332  svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows());
333  m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_workspace);
334  }
335  if(svd.computeU()) svd.m_matrixU.setIdentity(matrix.rows(), matrix.rows());
336  return true;
337  }
338  else return false;
339  }
340 
341 private:
343  QRType m_qr;
344  TransposeTypeWithSameStorageOrder m_adjoint;
346 };
347 
348 /*** 2x2 SVD implementation
349  ***
350  *** JacobiSVD consists in performing a series of 2x2 SVD subproblems
351  ***/
352 
353 template<typename MatrixType, int QRPreconditioner>
354 struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, false>
355 {
358  static bool run(typename SVD::WorkMatrixType&, SVD&, Index, Index, RealScalar&) { return true; }
359 };
360 
361 template<typename MatrixType, int QRPreconditioner>
362 struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, true>
363 {
365  typedef typename MatrixType::Scalar Scalar;
367  static bool run(typename SVD::WorkMatrixType& work_matrix, SVD& svd, Index p, Index q, RealScalar& maxDiagEntry)
368  {
369  using std::sqrt;
370  using std::abs;
371  Scalar z;
373  RealScalar n = sqrt(numext::abs2(work_matrix.coeff(p,p)) + numext::abs2(work_matrix.coeff(q,p)));
374 
375  const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)();
376  const RealScalar precision = NumTraits<Scalar>::epsilon();
377 
378  if(n==0)
379  {
380  // make sure first column is zero
381  work_matrix.coeffRef(p,p) = work_matrix.coeffRef(q,p) = Scalar(0);
382 
383  if(abs(numext::imag(work_matrix.coeff(p,q)))>considerAsZero)
384  {
385  // 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
386  z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q);
387  work_matrix.row(p) *= z;
388  if(svd.computeU()) svd.m_matrixU.col(p) *= conj(z);
389  }
390  if(abs(numext::imag(work_matrix.coeff(q,q)))>considerAsZero)
391  {
392  z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q);
393  work_matrix.row(q) *= z;
394  if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z);
395  }
396  // otherwise the second row is already zero, so we have nothing to do.
397  }
398  else
399  {
400  rot.c() = conj(work_matrix.coeff(p,p)) / n;
401  rot.s() = work_matrix.coeff(q,p) / n;
402  work_matrix.applyOnTheLeft(p,q,rot);
403  if(svd.computeU()) svd.m_matrixU.applyOnTheRight(p,q,rot.adjoint());
404  if(abs(numext::imag(work_matrix.coeff(p,q)))>considerAsZero)
405  {
406  z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q);
407  work_matrix.col(q) *= z;
408  if(svd.computeV()) svd.m_matrixV.col(q) *= z;
409  }
410  if(abs(numext::imag(work_matrix.coeff(q,q)))>considerAsZero)
411  {
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>(maxDiagEntry,numext::maxi<RealScalar>(abs(work_matrix.coeff(p,p)), abs(work_matrix.coeff(q,q))));
420  // and check whether the 2x2 block is already diagonal
421  RealScalar threshold = numext::maxi<RealScalar>(considerAsZero, precision * maxDiagEntry);
422  return abs(work_matrix.coeff(p,q))>threshold || abs(work_matrix.coeff(q,p)) > threshold;
423  }
424 };
425 
426 template<typename _MatrixType, int QRPreconditioner>
427 struct traits<JacobiSVD<_MatrixType,QRPreconditioner> >
428  : traits<_MatrixType>
429 {
430  typedef _MatrixType MatrixType;
431 };
432 
433 } // end namespace internal
434 
488 template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
489  : public SVDBase<JacobiSVD<_MatrixType,QRPreconditioner> >
490 {
492  public:
493 
494  typedef _MatrixType MatrixType;
495  typedef typename MatrixType::Scalar Scalar;
497  enum {
498  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
499  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
500  DiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_DYNAMIC(RowsAtCompileTime,ColsAtCompileTime),
501  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
502  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
503  MaxDiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(MaxRowsAtCompileTime,MaxColsAtCompileTime),
504  MatrixOptions = MatrixType::Options
505  };
506 
507  typedef typename Base::MatrixUType MatrixUType;
508  typedef typename Base::MatrixVType MatrixVType;
510 
513  typedef Matrix<Scalar, DiagSizeAtCompileTime, DiagSizeAtCompileTime,
514  MatrixOptions, MaxDiagSizeAtCompileTime, MaxDiagSizeAtCompileTime>
516 
523  {}
524 
525 
532  JacobiSVD(Index rows, Index cols, unsigned int computationOptions = 0)
533  {
534  allocate(rows, cols, computationOptions);
535  }
536 
547  explicit JacobiSVD(const MatrixType& matrix, unsigned int computationOptions = 0)
548  {
549  compute(matrix, computationOptions);
550  }
551 
562  JacobiSVD& compute(const MatrixType& matrix, unsigned int computationOptions);
563 
570  JacobiSVD& compute(const MatrixType& matrix)
571  {
572  return compute(matrix, m_computationOptions);
573  }
574 
575  using Base::computeU;
576  using Base::computeV;
577  using Base::rows;
578  using Base::cols;
579  using Base::rank;
580 
581  private:
582  void allocate(Index rows, Index cols, unsigned int computationOptions);
583 
584  protected:
585  using Base::m_matrixU;
586  using Base::m_matrixV;
587  using Base::m_singularValues;
588  using Base::m_info;
589  using Base::m_isInitialized;
590  using Base::m_isAllocated;
591  using Base::m_usePrescribedThreshold;
592  using Base::m_computeFullU;
593  using Base::m_computeThinU;
594  using Base::m_computeFullV;
595  using Base::m_computeThinV;
596  using Base::m_computationOptions;
597  using Base::m_nonzeroSingularValues;
598  using Base::m_rows;
599  using Base::m_cols;
600  using Base::m_diagSize;
601  using Base::m_prescribedThreshold;
603 
604  template<typename __MatrixType, int _QRPreconditioner, bool _IsComplex>
606  template<typename __MatrixType, int _QRPreconditioner, int _Case, bool _DoAnything>
608 
611  MatrixType m_scaledMatrix;
612 };
613 
614 template<typename MatrixType, int QRPreconditioner>
616 {
617  eigen_assert(rows >= 0 && cols >= 0);
618 
619  if (m_isAllocated &&
620  rows == m_rows &&
621  cols == m_cols &&
622  computationOptions == m_computationOptions)
623  {
624  return;
625  }
626 
627  m_rows = rows;
628  m_cols = cols;
629  m_info = Success;
630  m_isInitialized = false;
631  m_isAllocated = true;
632  m_computationOptions = computationOptions;
633  m_computeFullU = (computationOptions & ComputeFullU) != 0;
634  m_computeThinU = (computationOptions & ComputeThinU) != 0;
635  m_computeFullV = (computationOptions & ComputeFullV) != 0;
636  m_computeThinV = (computationOptions & ComputeThinV) != 0;
637  eigen_assert(!(m_computeFullU && m_computeThinU) && "JacobiSVD: you can't ask for both full and thin U");
638  eigen_assert(!(m_computeFullV && m_computeThinV) && "JacobiSVD: you can't ask for both full and thin V");
639  eigen_assert(EIGEN_IMPLIES(m_computeThinU || m_computeThinV, MatrixType::ColsAtCompileTime==Dynamic) &&
640  "JacobiSVD: thin U and V are only available when your matrix has a dynamic number of columns.");
641  if (QRPreconditioner == FullPivHouseholderQRPreconditioner)
642  {
643  eigen_assert(!(m_computeThinU || m_computeThinV) &&
644  "JacobiSVD: can't compute thin U or thin V with the FullPivHouseholderQR preconditioner. "
645  "Use the ColPivHouseholderQR preconditioner instead.");
646  }
647  m_diagSize = (std::min)(m_rows, m_cols);
648  m_singularValues.resize(m_diagSize);
649  if(RowsAtCompileTime==Dynamic)
650  m_matrixU.resize(m_rows, m_computeFullU ? m_rows
651  : m_computeThinU ? m_diagSize
652  : 0);
653  if(ColsAtCompileTime==Dynamic)
654  m_matrixV.resize(m_cols, m_computeFullV ? m_cols
655  : m_computeThinV ? m_diagSize
656  : 0);
657  m_workMatrix.resize(m_diagSize, m_diagSize);
658 
659  if(m_cols>m_rows) m_qr_precond_morecols.allocate(*this);
660  if(m_rows>m_cols) m_qr_precond_morerows.allocate(*this);
661  if(m_rows!=m_cols) m_scaledMatrix.resize(rows,cols);
662 }
663 
664 template<typename MatrixType, int QRPreconditioner>
667 {
668  using std::abs;
669  allocate(matrix.rows(), matrix.cols(), computationOptions);
670 
671  // currently we stop when we reach precision 2*epsilon as the last bit of precision can require an unreasonable number of iterations,
672  // only worsening the precision of U and V as we accumulate more rotations
674 
675  // limit for denormal numbers to be considered zero in order to avoid infinite loops (see bug 286)
676  const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)();
677 
678  // Scaling factor to reduce over/under-flows
679  RealScalar scale = matrix.cwiseAbs().template maxCoeff<PropagateNaN>();
680  if (!(numext::isfinite)(scale)) {
681  m_isInitialized = true;
682  m_info = InvalidInput;
683  return *this;
684  }
685  if(scale==RealScalar(0)) scale = RealScalar(1);
686 
687  /*** step 1. The R-SVD step: we use a QR decomposition to reduce to the case of a square matrix */
688 
689  if(m_rows!=m_cols)
690  {
691  m_scaledMatrix = matrix / scale;
692  m_qr_precond_morecols.run(*this, m_scaledMatrix);
693  m_qr_precond_morerows.run(*this, m_scaledMatrix);
694  }
695  else
696  {
697  m_workMatrix = matrix.block(0,0,m_diagSize,m_diagSize) / scale;
698  if(m_computeFullU) m_matrixU.setIdentity(m_rows,m_rows);
699  if(m_computeThinU) m_matrixU.setIdentity(m_rows,m_diagSize);
700  if(m_computeFullV) m_matrixV.setIdentity(m_cols,m_cols);
701  if(m_computeThinV) m_matrixV.setIdentity(m_cols, m_diagSize);
702  }
703 
704  /*** step 2. The main Jacobi SVD iteration. ***/
705  RealScalar maxDiagEntry = m_workMatrix.cwiseAbs().diagonal().maxCoeff();
706 
707  bool finished = false;
708  while(!finished)
709  {
710  finished = true;
711 
712  // do a sweep: for all index pairs (p,q), perform SVD of the corresponding 2x2 sub-matrix
713 
714  for(Index p = 1; p < m_diagSize; ++p)
715  {
716  for(Index q = 0; q < p; ++q)
717  {
718  // if this 2x2 sub-matrix is not diagonal already...
719  // notice that this comparison will evaluate to false if any NaN is involved, ensuring that NaN's don't
720  // keep us iterating forever. Similarly, small denormal numbers are considered zero.
721  RealScalar threshold = numext::maxi<RealScalar>(considerAsZero, precision * maxDiagEntry);
722  if(abs(m_workMatrix.coeff(p,q))>threshold || abs(m_workMatrix.coeff(q,p)) > threshold)
723  {
724  finished = false;
725  // perform SVD decomposition of 2x2 sub-matrix corresponding to indices p,q to make it diagonal
726  // the complex to real operation returns true if the updated 2x2 block is not already diagonal
728  {
729  JacobiRotation<RealScalar> j_left, j_right;
730  internal::real_2x2_jacobi_svd(m_workMatrix, p, q, &j_left, &j_right);
731 
732  // accumulate resulting Jacobi rotations
733  m_workMatrix.applyOnTheLeft(p,q,j_left);
734  if(computeU()) m_matrixU.applyOnTheRight(p,q,j_left.transpose());
735 
736  m_workMatrix.applyOnTheRight(p,q,j_right);
737  if(computeV()) m_matrixV.applyOnTheRight(p,q,j_right);
738 
739  // keep track of the largest diagonal coefficient
740  maxDiagEntry = numext::maxi<RealScalar>(maxDiagEntry,numext::maxi<RealScalar>(abs(m_workMatrix.coeff(p,p)), abs(m_workMatrix.coeff(q,q))));
741  }
742  }
743  }
744  }
745  }
746 
747  /*** step 3. The work matrix is now diagonal, so ensure it's positive so its diagonal entries are the singular values ***/
748 
749  for(Index i = 0; i < m_diagSize; ++i)
750  {
751  // For a complex matrix, some diagonal coefficients might note have been
752  // treated by svd_precondition_2x2_block_to_be_real, and the imaginary part
753  // of some diagonal entry might not be null.
754  if(NumTraits<Scalar>::IsComplex && abs(numext::imag(m_workMatrix.coeff(i,i)))>considerAsZero)
755  {
756  RealScalar a = abs(m_workMatrix.coeff(i,i));
757  m_singularValues.coeffRef(i) = abs(a);
758  if(computeU()) m_matrixU.col(i) *= m_workMatrix.coeff(i,i)/a;
759  }
760  else
761  {
762  // m_workMatrix.coeff(i,i) is already real, no difficulty:
763  RealScalar a = numext::real(m_workMatrix.coeff(i,i));
764  m_singularValues.coeffRef(i) = abs(a);
765  if(computeU() && (a<RealScalar(0))) m_matrixU.col(i) = -m_matrixU.col(i);
766  }
767  }
768 
769  m_singularValues *= scale;
770 
771  /*** step 4. Sort singular values in descending order and compute the number of nonzero singular values ***/
772 
773  m_nonzeroSingularValues = m_diagSize;
774  for(Index i = 0; i < m_diagSize; i++)
775  {
776  Index pos;
777  RealScalar maxRemainingSingularValue = m_singularValues.tail(m_diagSize-i).maxCoeff(&pos);
778  if(maxRemainingSingularValue == RealScalar(0))
779  {
780  m_nonzeroSingularValues = i;
781  break;
782  }
783  if(pos)
784  {
785  pos += i;
786  std::swap(m_singularValues.coeffRef(i), m_singularValues.coeffRef(pos));
787  if(computeU()) m_matrixU.col(pos).swap(m_matrixU.col(i));
788  if(computeV()) m_matrixV.col(pos).swap(m_matrixV.col(i));
789  }
790  }
791 
792  m_isInitialized = true;
793  return *this;
794 }
795 
803 template<typename Derived>
805 MatrixBase<Derived>::jacobiSvd(unsigned int computationOptions) const
806 {
807  return JacobiSVD<PlainObject>(*this, computationOptions);
808 }
809 
810 } // end namespace Eigen
811 
812 #endif // EIGEN_JACOBISVD_H
Matrix< Scalar, DiagSizeAtCompileTime, DiagSizeAtCompileTime, MatrixOptions, MaxDiagSizeAtCompileTime, MaxDiagSizeAtCompileTime > WorkMatrixType
Definition: JacobiSVD.h:515
internal::make_proper_matrix_type< Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime >::type TransposeTypeWithSameStorageOrder
Definition: JacobiSVD.h:307
bool run(JacobiSVD< MatrixType, QRPreconditioner > &, const MatrixType &)
Definition: JacobiSVD.h:56
int EIGEN_BLAS_FUNC() rot(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *pc, RealScalar *ps)
SCALAR Scalar
Definition: bench_gemm.cpp:46
cout<< "Here is the matrix m:"<< endl<< m<< endl;JacobiSVD< MatrixXf > svd(m, ComputeThinU|ComputeThinV)
const AutoDiffScalar< DerType > & conj(const AutoDiffScalar< DerType > &x)
internal::plain_col_type< MatrixType >::type ColType
Definition: JacobiSVD.h:512
static bool run(typename SVD::WorkMatrixType &, SVD &, Index, Index, RealScalar &)
Definition: JacobiSVD.h:358
float real
Definition: datatypes.h:10
Scalar * b
Definition: benchVecAdd.cpp:17
void allocate(const JacobiSVD< MatrixType, QRPreconditioner > &)
Definition: JacobiSVD.h:55
EIGEN_DEVICE_FUNC Scalar & c()
Definition: Jacobi.h:47
#define min(a, b)
Definition: datatypes.h:19
int n
SVDBase< JacobiSVD > Base
Definition: JacobiSVD.h:491
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
Rotation given by a cosine-sine pair.
MatrixXf MatrixType
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:232
static bool run(typename SVD::WorkMatrixType &work_matrix, SVD &svd, Index p, Index q, RealScalar &maxDiagEntry)
Definition: JacobiSVD.h:367
internal::make_proper_matrix_type< Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime >::type TransposeTypeWithSameStorageOrder
Definition: JacobiSVD.h:120
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:39
#define EIGEN_IMPLIES(a, b)
Definition: Macros.h:1315
#define EIGEN_SIZE_MIN_PREFER_FIXED(a, b)
Definition: Macros.h:1302
static double epsilon
Definition: testRot3.cpp:37
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar & coeffRef(Index rowId, Index colId)
MatrixType::Scalar Scalar
Definition: JacobiSVD.h:495
bool run(JacobiSVD< MatrixType, ColPivHouseholderQRPreconditioner > &svd, const MatrixType &matrix)
Definition: JacobiSVD.h:170
void real_2x2_jacobi_svd(const MatrixType &matrix, Index p, Index q, JacobiRotation< RealScalar > *j_left, JacobiRotation< RealScalar > *j_right)
Definition: RealSvd2x2.h:19
internal::make_proper_matrix_type< Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime >::type TransposeTypeWithSameStorageOrder
Definition: JacobiSVD.h:210
Base class of SVD algorithms.
JacobiSVD< PlainObject > jacobiSvd(unsigned int computationOptions=0) const
Definition: JacobiSVD.h:805
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar & coeff(Index rowId, Index colId) const
JacobiSVD()
Default Constructor.
Definition: JacobiSVD.h:522
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE bool() isfinite(const Eigen::bfloat16 &h)
Definition: BFloat16.h:671
#define eigen_assert(x)
Definition: Macros.h:1037
internal::qr_preconditioner_impl< MatrixType, QRPreconditioner, internal::PreconditionIfMoreColsThanRows > m_qr_precond_morecols
Definition: JacobiSVD.h:609
_MatrixType MatrixType
Definition: JacobiSVD.h:494
EIGEN_DEVICE_FUNC const Scalar & q
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:47
void swap(GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &a, GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &b)
EIGEN_DEVICE_FUNC JacobiRotation adjoint() const
Definition: Jacobi.h:67
DenseIndex ret
Base::MatrixVType MatrixVType
Definition: JacobiSVD.h:508
bool run(JacobiSVD< MatrixType, HouseholderQRPreconditioner > &svd, const MatrixType &matrix)
Definition: JacobiSVD.h:268
bool run(JacobiSVD< MatrixType, FullPivHouseholderQRPreconditioner > &svd, const MatrixType &matrix)
Definition: JacobiSVD.h:86
bool run(JacobiSVD< MatrixType, HouseholderQRPreconditioner > &svd, const MatrixType &matrix)
Definition: JacobiSVD.h:321
MatrixType m_scaledMatrix
Definition: JacobiSVD.h:611
EIGEN_DEVICE_FUNC Scalar & s()
Definition: Jacobi.h:49
JacobiSVD & compute(const MatrixType &matrix)
Method performing the decomposition of given matrix using current options.
Definition: JacobiSVD.h:570
EIGEN_DEVICE_FUNC JacobiRotation transpose() const
Definition: Jacobi.h:63
cout precision(2)
void allocate(Index rows, Index cols, unsigned int computationOptions)
Definition: JacobiSVD.h:615
Two-sided Jacobi SVD decomposition of a rectangular matrix.
float * p
internal::plain_row_type< MatrixType >::type RowType
Definition: JacobiSVD.h:511
NumTraits< typename MatrixType::Scalar >::Real RealScalar
Definition: JacobiSVD.h:496
JacobiSVD & compute(const MatrixType &matrix, unsigned int computationOptions)
Method performing the decomposition of given matrix using custom options.
Definition: JacobiSVD.h:666
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy y set format x g set format y g set format x2 g set format y2 g set format z g set angles radians set nogrid set key title set key left top Right noreverse box linetype linewidth samplen spacing width set nolabel set noarrow set nologscale set logscale x set set pointsize set encoding default set nopolar set noparametric set set set set surface set nocontour set clabel set mapping cartesian set nohidden3d set cntrparam order set cntrparam linear set cntrparam levels auto set cntrparam points set size set set xzeroaxis lt lw set x2zeroaxis lt lw set yzeroaxis lt lw set y2zeroaxis lt lw set tics in set ticslevel set tics scale
bool run(JacobiSVD< MatrixType, FullPivHouseholderQRPreconditioner > &svd, const MatrixType &matrix)
Definition: JacobiSVD.h:133
EIGEN_DEVICE_FUNC const ImagReturnType imag() const
Base::SingularValuesType SingularValuesType
Definition: JacobiSVD.h:509
Jet< T, N > sqrt(const Jet< T, N > &f)
Definition: jet.h:418
internal::qr_preconditioner_impl< MatrixType, QRPreconditioner, internal::PreconditionIfMoreRowsThanCols > m_qr_precond_morerows
Definition: JacobiSVD.h:610
const int Dynamic
Definition: Constants.h:22
Generic expression where a coefficient-wise unary operator is applied to an expression.
Definition: CwiseUnaryOp.h:55
EIGEN_DONT_INLINE void compute(Solver &solver, const MatrixType &A)
Map< Matrix< T, Dynamic, Dynamic, ColMajor >, 0, OuterStride<> > matrix(T *data, int rows, int cols, int stride)
#define EIGEN_SIZE_MIN_PREFER_DYNAMIC(a, b)
Definition: Macros.h:1294
#define abs(x)
Definition: datatypes.h:17
JacobiSVD(Index rows, Index cols, unsigned int computationOptions=0)
Default Constructor with memory preallocation.
Definition: JacobiSVD.h:532
EIGEN_DEVICE_FUNC bool abs2(bool x)
bool run(JacobiSVD< MatrixType, ColPivHouseholderQRPreconditioner > &svd, const MatrixType &matrix)
Definition: JacobiSVD.h:224
WorkMatrixType m_workMatrix
Definition: JacobiSVD.h:602
JacobiSVD(const MatrixType &matrix, unsigned int computationOptions=0)
Constructor performing the decomposition of given matrix.
Definition: JacobiSVD.h:547
Definition: pytypes.h:1370
Base::MatrixUType MatrixUType
Definition: JacobiSVD.h:507


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:34:29