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 // Eigen is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 3 of the License, or (at your option) any later version.
10 //
11 // Alternatively, you can redistribute it and/or
12 // modify it under the terms of the GNU General Public License as
13 // published by the Free Software Foundation; either version 2 of
14 // the License, or (at your option) any later version.
15 //
16 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
17 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
19 // GNU General Public License for more details.
20 //
21 // You should have received a copy of the GNU Lesser General Public
22 // License and a copy of the GNU General Public License along with
23 // Eigen. If not, see <http://www.gnu.org/licenses/>.
24 
25 #ifndef EIGEN_JACOBISVD_H
26 #define EIGEN_JACOBISVD_H
27 
28 namespace Eigen {
29 
30 namespace internal {
31 // forward declaration (needed by ICC)
32 // the empty body is required by MSVC
33 template<typename MatrixType, int QRPreconditioner,
35 struct svd_precondition_2x2_block_to_be_real {};
36 
37 /*** QR preconditioners (R-SVD)
38  ***
39  *** Their role is to reduce the problem of computing the SVD to the case of a square matrix.
40  *** This approach, known as R-SVD, is an optimization for rectangular-enough matrices, and is a requirement for
41  *** JacobiSVD which by itself is only able to work on square matrices.
42  ***/
43 
45 
46 template<typename MatrixType, int QRPreconditioner, int Case>
47 struct qr_preconditioner_should_do_anything
48 {
49  enum { a = MatrixType::RowsAtCompileTime != Dynamic &&
50  MatrixType::ColsAtCompileTime != Dynamic &&
51  MatrixType::ColsAtCompileTime <= MatrixType::RowsAtCompileTime,
52  b = MatrixType::RowsAtCompileTime != Dynamic &&
53  MatrixType::ColsAtCompileTime != Dynamic &&
54  MatrixType::RowsAtCompileTime <= MatrixType::ColsAtCompileTime,
55  ret = !( (QRPreconditioner == NoQRPreconditioner) ||
56  (Case == PreconditionIfMoreColsThanRows && bool(a)) ||
57  (Case == PreconditionIfMoreRowsThanCols && bool(b)) )
58  };
59 };
60 
61 template<typename MatrixType, int QRPreconditioner, int Case,
62  bool DoAnything = qr_preconditioner_should_do_anything<MatrixType, QRPreconditioner, Case>::ret
63 > struct qr_preconditioner_impl {};
64 
65 template<typename MatrixType, int QRPreconditioner, int Case>
66 class qr_preconditioner_impl<MatrixType, QRPreconditioner, Case, false>
67 {
68 public:
69  typedef typename MatrixType::Index Index;
70  void allocate(const JacobiSVD<MatrixType, QRPreconditioner>&) {}
71  bool run(JacobiSVD<MatrixType, QRPreconditioner>&, const MatrixType&)
72  {
73  return false;
74  }
75 };
76 
77 /*** preconditioner using FullPivHouseholderQR ***/
78 
79 template<typename MatrixType>
80 class qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
81 {
82 public:
83  typedef typename MatrixType::Index Index;
84  typedef typename MatrixType::Scalar Scalar;
85  enum
86  {
87  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
88  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime
89  };
90  typedef Matrix<Scalar, 1, RowsAtCompileTime, RowMajor, 1, MaxRowsAtCompileTime> WorkspaceType;
91 
92  void allocate(const JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd)
93  {
94  if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols())
95  {
96  m_qr = FullPivHouseholderQR<MatrixType>(svd.rows(), svd.cols());
97  }
98  if (svd.m_computeFullU) m_workspace.resize(svd.rows());
99  }
100 
101  bool run(JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
102  {
103  if(matrix.rows() > matrix.cols())
104  {
105  m_qr.compute(matrix);
106  svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>();
107  if(svd.m_computeFullU) m_qr.matrixQ().evalTo(svd.m_matrixU, m_workspace);
108  if(svd.computeV()) svd.m_matrixV = m_qr.colsPermutation();
109  return true;
110  }
111  return false;
112  }
113 private:
114  FullPivHouseholderQR<MatrixType> m_qr;
115  WorkspaceType m_workspace;
116 };
117 
118 template<typename MatrixType>
119 class qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true>
120 {
121 public:
122  typedef typename MatrixType::Index Index;
123  typedef typename MatrixType::Scalar Scalar;
124  enum
125  {
126  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
127  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
128  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
129  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
130  Options = MatrixType::Options
131  };
132  typedef Matrix<Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime>
133  TransposeTypeWithSameStorageOrder;
134 
135  void allocate(const JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd)
136  {
137  if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols())
138  {
139  m_qr = FullPivHouseholderQR<TransposeTypeWithSameStorageOrder>(svd.cols(), svd.rows());
140  }
141  m_adjoint.resize(svd.cols(), svd.rows());
142  if (svd.m_computeFullV) m_workspace.resize(svd.cols());
143  }
144 
145  bool run(JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
146  {
147  if(matrix.cols() > matrix.rows())
148  {
149  m_adjoint = matrix.adjoint();
150  m_qr.compute(m_adjoint);
151  svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint();
152  if(svd.m_computeFullV) m_qr.matrixQ().evalTo(svd.m_matrixV, m_workspace);
153  if(svd.computeU()) svd.m_matrixU = m_qr.colsPermutation();
154  return true;
155  }
156  else return false;
157  }
158 private:
159  FullPivHouseholderQR<TransposeTypeWithSameStorageOrder> m_qr;
160  TransposeTypeWithSameStorageOrder m_adjoint;
161  typename internal::plain_row_type<MatrixType>::type m_workspace;
162 };
163 
164 /*** preconditioner using ColPivHouseholderQR ***/
165 
166 template<typename MatrixType>
167 class qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
168 {
169 public:
170  typedef typename MatrixType::Index Index;
171 
172  void allocate(const JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd)
173  {
174  if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols())
175  {
176  m_qr = ColPivHouseholderQR<MatrixType>(svd.rows(), svd.cols());
177  }
178  if (svd.m_computeFullU) m_workspace.resize(svd.rows());
179  else if (svd.m_computeThinU) m_workspace.resize(svd.cols());
180  }
181 
182  bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
183  {
184  if(matrix.rows() > matrix.cols())
185  {
186  m_qr.compute(matrix);
187  svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>();
188  if(svd.m_computeFullU) m_qr.householderQ().evalTo(svd.m_matrixU, m_workspace);
189  else if(svd.m_computeThinU)
190  {
191  svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols());
192  m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_workspace);
193  }
194  if(svd.computeV()) svd.m_matrixV = m_qr.colsPermutation();
195  return true;
196  }
197  return false;
198  }
199 
200 private:
201  ColPivHouseholderQR<MatrixType> m_qr;
202  typename internal::plain_col_type<MatrixType>::type m_workspace;
203 };
204 
205 template<typename MatrixType>
206 class qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true>
207 {
208 public:
209  typedef typename MatrixType::Index Index;
210  typedef typename MatrixType::Scalar Scalar;
211  enum
212  {
213  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
214  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
215  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
216  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
217  Options = MatrixType::Options
218  };
219 
220  typedef Matrix<Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime>
221  TransposeTypeWithSameStorageOrder;
222 
223  void allocate(const JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd)
224  {
225  if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols())
226  {
227  m_qr = ColPivHouseholderQR<TransposeTypeWithSameStorageOrder>(svd.cols(), svd.rows());
228  }
229  if (svd.m_computeFullV) m_workspace.resize(svd.cols());
230  else if (svd.m_computeThinV) m_workspace.resize(svd.rows());
231  m_adjoint.resize(svd.cols(), svd.rows());
232  }
233 
234  bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
235  {
236  if(matrix.cols() > matrix.rows())
237  {
238  m_adjoint = matrix.adjoint();
239  m_qr.compute(m_adjoint);
240 
241  svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint();
242  if(svd.m_computeFullV) m_qr.householderQ().evalTo(svd.m_matrixV, m_workspace);
243  else if(svd.m_computeThinV)
244  {
245  svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows());
246  m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_workspace);
247  }
248  if(svd.computeU()) svd.m_matrixU = m_qr.colsPermutation();
249  return true;
250  }
251  else return false;
252  }
253 
254 private:
255  ColPivHouseholderQR<TransposeTypeWithSameStorageOrder> m_qr;
256  TransposeTypeWithSameStorageOrder m_adjoint;
257  typename internal::plain_row_type<MatrixType>::type m_workspace;
258 };
259 
260 /*** preconditioner using HouseholderQR ***/
261 
262 template<typename MatrixType>
263 class qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
264 {
265 public:
266  typedef typename MatrixType::Index Index;
267 
268  void allocate(const JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd)
269  {
270  if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols())
271  {
272  m_qr = HouseholderQR<MatrixType>(svd.rows(), svd.cols());
273  }
274  if (svd.m_computeFullU) m_workspace.resize(svd.rows());
275  else if (svd.m_computeThinU) m_workspace.resize(svd.cols());
276  }
277 
278  bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix)
279  {
280  if(matrix.rows() > matrix.cols())
281  {
282  m_qr.compute(matrix);
283  svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>();
284  if(svd.m_computeFullU) m_qr.householderQ().evalTo(svd.m_matrixU, m_workspace);
285  else if(svd.m_computeThinU)
286  {
287  svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols());
288  m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_workspace);
289  }
290  if(svd.computeV()) svd.m_matrixV.setIdentity(matrix.cols(), matrix.cols());
291  return true;
292  }
293  return false;
294  }
295 private:
296  HouseholderQR<MatrixType> m_qr;
297  typename internal::plain_col_type<MatrixType>::type m_workspace;
298 };
299 
300 template<typename MatrixType>
301 class qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true>
302 {
303 public:
304  typedef typename MatrixType::Index Index;
305  typedef typename MatrixType::Scalar Scalar;
306  enum
307  {
308  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
309  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
310  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
311  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
312  Options = MatrixType::Options
313  };
314 
315  typedef Matrix<Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime>
316  TransposeTypeWithSameStorageOrder;
317 
318  void allocate(const JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd)
319  {
320  if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols())
321  {
322  m_qr = HouseholderQR<TransposeTypeWithSameStorageOrder>(svd.cols(), svd.rows());
323  }
324  if (svd.m_computeFullV) m_workspace.resize(svd.cols());
325  else if (svd.m_computeThinV) m_workspace.resize(svd.rows());
326  m_adjoint.resize(svd.cols(), svd.rows());
327  }
328 
329  bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix)
330  {
331  if(matrix.cols() > matrix.rows())
332  {
333  m_adjoint = matrix.adjoint();
334  m_qr.compute(m_adjoint);
335 
336  svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint();
337  if(svd.m_computeFullV) m_qr.householderQ().evalTo(svd.m_matrixV, m_workspace);
338  else if(svd.m_computeThinV)
339  {
340  svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows());
341  m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_workspace);
342  }
343  if(svd.computeU()) svd.m_matrixU.setIdentity(matrix.rows(), matrix.rows());
344  return true;
345  }
346  else return false;
347  }
348 
349 private:
350  HouseholderQR<TransposeTypeWithSameStorageOrder> m_qr;
351  TransposeTypeWithSameStorageOrder m_adjoint;
352  typename internal::plain_row_type<MatrixType>::type m_workspace;
353 };
354 
355 /*** 2x2 SVD implementation
356  ***
357  *** JacobiSVD consists in performing a series of 2x2 SVD subproblems
358  ***/
359 
360 template<typename MatrixType, int QRPreconditioner>
361 struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, false>
362 {
363  typedef JacobiSVD<MatrixType, QRPreconditioner> SVD;
364  typedef typename SVD::Index Index;
365  static void run(typename SVD::WorkMatrixType&, SVD&, Index, Index) {}
366 };
367 
368 template<typename MatrixType, int QRPreconditioner>
369 struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, true>
370 {
371  typedef JacobiSVD<MatrixType, QRPreconditioner> SVD;
372  typedef typename MatrixType::Scalar Scalar;
373  typedef typename MatrixType::RealScalar RealScalar;
374  typedef typename SVD::Index Index;
375  static void run(typename SVD::WorkMatrixType& work_matrix, SVD& svd, Index p, Index q)
376  {
377  Scalar z;
378  JacobiRotation<Scalar> rot;
379  RealScalar n = sqrt(abs2(work_matrix.coeff(p,p)) + abs2(work_matrix.coeff(q,p)));
380  if(n==0)
381  {
382  z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q);
383  work_matrix.row(p) *= z;
384  if(svd.computeU()) svd.m_matrixU.col(p) *= conj(z);
385  z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q);
386  work_matrix.row(q) *= z;
387  if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z);
388  }
389  else
390  {
391  rot.c() = conj(work_matrix.coeff(p,p)) / n;
392  rot.s() = work_matrix.coeff(q,p) / n;
393  work_matrix.applyOnTheLeft(p,q,rot);
394  if(svd.computeU()) svd.m_matrixU.applyOnTheRight(p,q,rot.adjoint());
395  if(work_matrix.coeff(p,q) != Scalar(0))
396  {
397  Scalar z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q);
398  work_matrix.col(q) *= z;
399  if(svd.computeV()) svd.m_matrixV.col(q) *= z;
400  }
401  if(work_matrix.coeff(q,q) != Scalar(0))
402  {
403  z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q);
404  work_matrix.row(q) *= z;
405  if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z);
406  }
407  }
408  }
409 };
410 
411 template<typename MatrixType, typename RealScalar, typename Index>
412 void real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q,
415 {
417  m << real(matrix.coeff(p,p)), real(matrix.coeff(p,q)),
418  real(matrix.coeff(q,p)), real(matrix.coeff(q,q));
420  RealScalar t = m.coeff(0,0) + m.coeff(1,1);
421  RealScalar d = m.coeff(1,0) - m.coeff(0,1);
422  if(t == RealScalar(0))
423  {
424  rot1.c() = RealScalar(0);
425  rot1.s() = d > RealScalar(0) ? RealScalar(1) : RealScalar(-1);
426  }
427  else
428  {
429  RealScalar u = d / t;
430  rot1.c() = RealScalar(1) / sqrt(RealScalar(1) + abs2(u));
431  rot1.s() = rot1.c() * u;
432  }
433  m.applyOnTheLeft(0,1,rot1);
434  j_right->makeJacobi(m,0,1);
435  *j_left = rot1 * j_right->transpose();
436 }
437 
438 } // end namespace internal
439 
493 template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
494 {
495  public:
496 
497  typedef _MatrixType MatrixType;
498  typedef typename MatrixType::Scalar Scalar;
500  typedef typename MatrixType::Index Index;
501  enum {
508  MatrixOptions = MatrixType::Options
509  };
510 
517  typedef typename internal::plain_diag_type<MatrixType, RealScalar>::type SingularValuesType;
518  typedef typename internal::plain_row_type<MatrixType>::type RowType;
519  typedef typename internal::plain_col_type<MatrixType>::type ColType;
523 
530  : m_isInitialized(false),
531  m_isAllocated(false),
533  m_rows(-1), m_cols(-1)
534  {}
535 
536 
543  JacobiSVD(Index rows, Index cols, unsigned int computationOptions = 0)
544  : m_isInitialized(false),
545  m_isAllocated(false),
547  m_rows(-1), m_cols(-1)
548  {
549  allocate(rows, cols, computationOptions);
550  }
551 
562  JacobiSVD(const MatrixType& matrix, unsigned int computationOptions = 0)
563  : m_isInitialized(false),
564  m_isAllocated(false),
566  m_rows(-1), m_cols(-1)
567  {
568  compute(matrix, computationOptions);
569  }
570 
581  JacobiSVD& compute(const MatrixType& matrix, unsigned int computationOptions);
582 
589  JacobiSVD& compute(const MatrixType& matrix)
590  {
591  return compute(matrix, m_computationOptions);
592  }
593 
603  const MatrixUType& matrixU() const
604  {
605  eigen_assert(m_isInitialized && "JacobiSVD is not initialized.");
606  eigen_assert(computeU() && "This JacobiSVD decomposition didn't compute U. Did you ask for it?");
607  return m_matrixU;
608  }
609 
619  const MatrixVType& matrixV() const
620  {
621  eigen_assert(m_isInitialized && "JacobiSVD is not initialized.");
622  eigen_assert(computeV() && "This JacobiSVD decomposition didn't compute V. Did you ask for it?");
623  return m_matrixV;
624  }
625 
632  {
633  eigen_assert(m_isInitialized && "JacobiSVD is not initialized.");
634  return m_singularValues;
635  }
636 
638  inline bool computeU() const { return m_computeFullU || m_computeThinU; }
640  inline bool computeV() const { return m_computeFullV || m_computeThinV; }
641 
651  template<typename Rhs>
652  inline const internal::solve_retval<JacobiSVD, Rhs>
653  solve(const MatrixBase<Rhs>& b) const
654  {
655  eigen_assert(m_isInitialized && "JacobiSVD is not initialized.");
656  eigen_assert(computeU() && computeV() && "JacobiSVD::solve() requires both unitaries U and V to be computed (thin unitaries suffice).");
657  return internal::solve_retval<JacobiSVD, Rhs>(*this, b.derived());
658  }
659 
662  {
663  eigen_assert(m_isInitialized && "JacobiSVD is not initialized.");
665  }
666 
667  inline Index rows() const { return m_rows; }
668  inline Index cols() const { return m_cols; }
669 
670  private:
671  void allocate(Index rows, Index cols, unsigned int computationOptions);
672 
673  protected:
681  unsigned int m_computationOptions;
683 
684  template<typename __MatrixType, int _QRPreconditioner, bool _IsComplex>
685  friend struct internal::svd_precondition_2x2_block_to_be_real;
686  template<typename __MatrixType, int _QRPreconditioner, int _Case, bool _DoAnything>
687  friend struct internal::qr_preconditioner_impl;
688 
689  internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreColsThanRows> m_qr_precond_morecols;
690  internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreRowsThanCols> m_qr_precond_morerows;
691 };
692 
693 template<typename MatrixType, int QRPreconditioner>
694 void JacobiSVD<MatrixType, QRPreconditioner>::allocate(Index rows, Index cols, unsigned int computationOptions)
695 {
696  eigen_assert(rows >= 0 && cols >= 0);
697 
698  if (m_isAllocated &&
699  rows == m_rows &&
700  cols == m_cols &&
701  computationOptions == m_computationOptions)
702  {
703  return;
704  }
705 
706  m_rows = rows;
707  m_cols = cols;
708  m_isInitialized = false;
709  m_isAllocated = true;
710  m_computationOptions = computationOptions;
711  m_computeFullU = (computationOptions & ComputeFullU) != 0;
712  m_computeThinU = (computationOptions & ComputeThinU) != 0;
713  m_computeFullV = (computationOptions & ComputeFullV) != 0;
714  m_computeThinV = (computationOptions & ComputeThinV) != 0;
715  eigen_assert(!(m_computeFullU && m_computeThinU) && "JacobiSVD: you can't ask for both full and thin U");
716  eigen_assert(!(m_computeFullV && m_computeThinV) && "JacobiSVD: you can't ask for both full and thin V");
717  eigen_assert(EIGEN_IMPLIES(m_computeThinU || m_computeThinV, MatrixType::ColsAtCompileTime==Dynamic) &&
718  "JacobiSVD: thin U and V are only available when your matrix has a dynamic number of columns.");
719  if (QRPreconditioner == FullPivHouseholderQRPreconditioner)
720  {
721  eigen_assert(!(m_computeThinU || m_computeThinV) &&
722  "JacobiSVD: can't compute thin U or thin V with the FullPivHouseholderQR preconditioner. "
723  "Use the ColPivHouseholderQR preconditioner instead.");
724  }
725  m_diagSize = (std::min)(m_rows, m_cols);
726  m_singularValues.resize(m_diagSize);
727  m_matrixU.resize(m_rows, m_computeFullU ? m_rows
728  : m_computeThinU ? m_diagSize
729  : 0);
730  m_matrixV.resize(m_cols, m_computeFullV ? m_cols
731  : m_computeThinV ? m_diagSize
732  : 0);
733  m_workMatrix.resize(m_diagSize, m_diagSize);
734 
735  if(m_cols>m_rows) m_qr_precond_morecols.allocate(*this);
736  if(m_rows>m_cols) m_qr_precond_morerows.allocate(*this);
737 }
738 
739 template<typename MatrixType, int QRPreconditioner>
740 JacobiSVD<MatrixType, QRPreconditioner>&
741 JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsigned int computationOptions)
742 {
743  allocate(matrix.rows(), matrix.cols(), computationOptions);
744 
745  // currently we stop when we reach precision 2*epsilon as the last bit of precision can require an unreasonable number of iterations,
746  // only worsening the precision of U and V as we accumulate more rotations
747  const RealScalar precision = RealScalar(2) * NumTraits<Scalar>::epsilon();
748 
749  // limit for very small denormal numbers to be considered zero in order to avoid infinite loops (see bug 286)
750  const RealScalar considerAsZero = RealScalar(2) * std::numeric_limits<RealScalar>::denorm_min();
751 
752  /*** step 1. The R-SVD step: we use a QR decomposition to reduce to the case of a square matrix */
753 
754  if(!m_qr_precond_morecols.run(*this, matrix) && !m_qr_precond_morerows.run(*this, matrix))
755  {
756  m_workMatrix = matrix.block(0,0,m_diagSize,m_diagSize);
757  if(m_computeFullU) m_matrixU.setIdentity(m_rows,m_rows);
758  if(m_computeThinU) m_matrixU.setIdentity(m_rows,m_diagSize);
759  if(m_computeFullV) m_matrixV.setIdentity(m_cols,m_cols);
760  if(m_computeThinV) m_matrixV.setIdentity(m_cols, m_diagSize);
761  }
762 
763  /*** step 2. The main Jacobi SVD iteration. ***/
764 
765  bool finished = false;
766  while(!finished)
767  {
768  finished = true;
769 
770  // do a sweep: for all index pairs (p,q), perform SVD of the corresponding 2x2 sub-matrix
771 
772  for(Index p = 1; p < m_diagSize; ++p)
773  {
774  for(Index q = 0; q < p; ++q)
775  {
776  // if this 2x2 sub-matrix is not diagonal already...
777  // notice that this comparison will evaluate to false if any NaN is involved, ensuring that NaN's don't
778  // keep us iterating forever. Similarly, small denormal numbers are considered zero.
779  using std::max;
780  RealScalar threshold = (max)(considerAsZero, precision * (max)(internal::abs(m_workMatrix.coeff(p,p)),
781  internal::abs(m_workMatrix.coeff(q,q))));
782  if((max)(internal::abs(m_workMatrix.coeff(p,q)),internal::abs(m_workMatrix.coeff(q,p))) > threshold)
783  {
784  finished = false;
785 
786  // perform SVD decomposition of 2x2 sub-matrix corresponding to indices p,q to make it diagonal
787  internal::svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner>::run(m_workMatrix, *this, p, q);
788  JacobiRotation<RealScalar> j_left, j_right;
789  internal::real_2x2_jacobi_svd(m_workMatrix, p, q, &j_left, &j_right);
790 
791  // accumulate resulting Jacobi rotations
792  m_workMatrix.applyOnTheLeft(p,q,j_left);
793  if(computeU()) m_matrixU.applyOnTheRight(p,q,j_left.transpose());
794 
795  m_workMatrix.applyOnTheRight(p,q,j_right);
796  if(computeV()) m_matrixV.applyOnTheRight(p,q,j_right);
797  }
798  }
799  }
800  }
801 
802  /*** step 3. The work matrix is now diagonal, so ensure it's positive so its diagonal entries are the singular values ***/
803 
804  for(Index i = 0; i < m_diagSize; ++i)
805  {
806  RealScalar a = internal::abs(m_workMatrix.coeff(i,i));
807  m_singularValues.coeffRef(i) = a;
808  if(computeU() && (a!=RealScalar(0))) m_matrixU.col(i) *= m_workMatrix.coeff(i,i)/a;
809  }
810 
811  /*** step 4. Sort singular values in descending order and compute the number of nonzero singular values ***/
812 
813  m_nonzeroSingularValues = m_diagSize;
814  for(Index i = 0; i < m_diagSize; i++)
815  {
816  Index pos;
817  RealScalar maxRemainingSingularValue = m_singularValues.tail(m_diagSize-i).maxCoeff(&pos);
818  if(maxRemainingSingularValue == RealScalar(0))
819  {
820  m_nonzeroSingularValues = i;
821  break;
822  }
823  if(pos)
824  {
825  pos += i;
826  std::swap(m_singularValues.coeffRef(i), m_singularValues.coeffRef(pos));
827  if(computeU()) m_matrixU.col(pos).swap(m_matrixU.col(i));
828  if(computeV()) m_matrixV.col(pos).swap(m_matrixV.col(i));
829  }
830  }
831 
832  m_isInitialized = true;
833  return *this;
834 }
835 
836 namespace internal {
837 template<typename _MatrixType, int QRPreconditioner, typename Rhs>
838 struct solve_retval<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs>
839  : solve_retval_base<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs>
840 {
841  typedef JacobiSVD<_MatrixType, QRPreconditioner> JacobiSVDType;
842  EIGEN_MAKE_SOLVE_HELPERS(JacobiSVDType,Rhs)
843 
844  template<typename Dest> void evalTo(Dest& dst) const
845  {
846  eigen_assert(rhs().rows() == dec().rows());
847 
848  // A = U S V^*
849  // So A^{-1} = V S^{-1} U^*
850 
851  Index diagSize = (std::min)(dec().rows(), dec().cols());
852  typename JacobiSVDType::SingularValuesType invertedSingVals(diagSize);
853 
854  Index nonzeroSingVals = dec().nonzeroSingularValues();
855  invertedSingVals.head(nonzeroSingVals) = dec().singularValues().head(nonzeroSingVals).array().inverse();
856  invertedSingVals.tail(diagSize - nonzeroSingVals).setZero();
857 
858  dst = dec().matrixV().leftCols(diagSize)
859  * invertedSingVals.asDiagonal()
860  * dec().matrixU().leftCols(diagSize).adjoint()
861  * rhs();
862  }
863 };
864 } // end namespace internal
865 
873 template<typename Derived>
874 JacobiSVD<typename MatrixBase<Derived>::PlainObject>
875 MatrixBase<Derived>::jacobiSvd(unsigned int computationOptions) const
876 {
877  return JacobiSVD<PlainObject>(*this, computationOptions);
878 }
879 
880 } // end namespace Eigen
881 
882 #endif // EIGEN_JACOBISVD_H