JacobiSVD.h revision c981c48f5bc9aefeffc0bcb0cc3934c2fae179dd
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
13namespace Eigen {
14
15namespace internal {
16// forward declaration (needed by ICC)
17// the empty body is required by MSVC
18template<typename MatrixType, int QRPreconditioner,
19         bool IsComplex = NumTraits<typename MatrixType::Scalar>::IsComplex>
20struct svd_precondition_2x2_block_to_be_real {};
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
29enum { PreconditionIfMoreColsThanRows, PreconditionIfMoreRowsThanCols };
30
31template<typename MatrixType, int QRPreconditioner, int Case>
32struct qr_preconditioner_should_do_anything
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
46template<typename MatrixType, int QRPreconditioner, int Case,
47         bool DoAnything = qr_preconditioner_should_do_anything<MatrixType, QRPreconditioner, Case>::ret
48> struct qr_preconditioner_impl {};
49
50template<typename MatrixType, int QRPreconditioner, int Case>
51class qr_preconditioner_impl<MatrixType, QRPreconditioner, Case, false>
52{
53public:
54  typedef typename MatrixType::Index Index;
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
64template<typename MatrixType>
65class qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
66{
67public:
68  typedef typename MatrixType::Index Index;
69  typedef typename MatrixType::Scalar Scalar;
70  enum
71  {
72    RowsAtCompileTime = MatrixType::RowsAtCompileTime,
73    MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime
74  };
75  typedef Matrix<Scalar, 1, RowsAtCompileTime, RowMajor, 1, MaxRowsAtCompileTime> WorkspaceType;
76
77  void allocate(const JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd)
78  {
79    if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols())
80    {
81      m_qr = FullPivHouseholderQR<MatrixType>(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  }
98private:
99  FullPivHouseholderQR<MatrixType> m_qr;
100  WorkspaceType m_workspace;
101};
102
103template<typename MatrixType>
104class qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true>
105{
106public:
107  typedef typename MatrixType::Index Index;
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 = FullPivHouseholderQR<TransposeTypeWithSameStorageOrder>(svd.cols(), svd.rows());
125    }
126    m_adjoint.resize(svd.cols(), svd.rows());
127    if (svd.m_computeFullV) m_workspace.resize(svd.cols());
128  }
129
130  bool run(JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
131  {
132    if(matrix.cols() > matrix.rows())
133    {
134      m_adjoint = matrix.adjoint();
135      m_qr.compute(m_adjoint);
136      svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint();
137      if(svd.m_computeFullV) m_qr.matrixQ().evalTo(svd.m_matrixV, m_workspace);
138      if(svd.computeU()) svd.m_matrixU = m_qr.colsPermutation();
139      return true;
140    }
141    else return false;
142  }
143private:
144  FullPivHouseholderQR<TransposeTypeWithSameStorageOrder> m_qr;
145  TransposeTypeWithSameStorageOrder m_adjoint;
146  typename internal::plain_row_type<MatrixType>::type m_workspace;
147};
148
149/*** preconditioner using ColPivHouseholderQR ***/
150
151template<typename MatrixType>
152class qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
153{
154public:
155  typedef typename MatrixType::Index Index;
156
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 = ColPivHouseholderQR<MatrixType>(svd.rows(), svd.cols());
162    }
163    if (svd.m_computeFullU) m_workspace.resize(svd.rows());
164    else if (svd.m_computeThinU) m_workspace.resize(svd.cols());
165  }
166
167  bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
168  {
169    if(matrix.rows() > matrix.cols())
170    {
171      m_qr.compute(matrix);
172      svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>();
173      if(svd.m_computeFullU) m_qr.householderQ().evalTo(svd.m_matrixU, m_workspace);
174      else if(svd.m_computeThinU)
175      {
176        svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols());
177        m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_workspace);
178      }
179      if(svd.computeV()) svd.m_matrixV = m_qr.colsPermutation();
180      return true;
181    }
182    return false;
183  }
184
185private:
186  ColPivHouseholderQR<MatrixType> m_qr;
187  typename internal::plain_col_type<MatrixType>::type m_workspace;
188};
189
190template<typename MatrixType>
191class qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true>
192{
193public:
194  typedef typename MatrixType::Index Index;
195  typedef typename MatrixType::Scalar Scalar;
196  enum
197  {
198    RowsAtCompileTime = MatrixType::RowsAtCompileTime,
199    ColsAtCompileTime = MatrixType::ColsAtCompileTime,
200    MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
201    MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
202    Options = MatrixType::Options
203  };
204
205  typedef Matrix<Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime>
206          TransposeTypeWithSameStorageOrder;
207
208  void allocate(const JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd)
209  {
210    if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols())
211    {
212      m_qr = ColPivHouseholderQR<TransposeTypeWithSameStorageOrder>(svd.cols(), svd.rows());
213    }
214    if (svd.m_computeFullV) m_workspace.resize(svd.cols());
215    else if (svd.m_computeThinV) m_workspace.resize(svd.rows());
216    m_adjoint.resize(svd.cols(), svd.rows());
217  }
218
219  bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
220  {
221    if(matrix.cols() > matrix.rows())
222    {
223      m_adjoint = matrix.adjoint();
224      m_qr.compute(m_adjoint);
225
226      svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint();
227      if(svd.m_computeFullV) m_qr.householderQ().evalTo(svd.m_matrixV, m_workspace);
228      else if(svd.m_computeThinV)
229      {
230        svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows());
231        m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_workspace);
232      }
233      if(svd.computeU()) svd.m_matrixU = m_qr.colsPermutation();
234      return true;
235    }
236    else return false;
237  }
238
239private:
240  ColPivHouseholderQR<TransposeTypeWithSameStorageOrder> m_qr;
241  TransposeTypeWithSameStorageOrder m_adjoint;
242  typename internal::plain_row_type<MatrixType>::type m_workspace;
243};
244
245/*** preconditioner using HouseholderQR ***/
246
247template<typename MatrixType>
248class qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
249{
250public:
251  typedef typename MatrixType::Index Index;
252
253  void allocate(const JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd)
254  {
255    if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols())
256    {
257      m_qr = HouseholderQR<MatrixType>(svd.rows(), svd.cols());
258    }
259    if (svd.m_computeFullU) m_workspace.resize(svd.rows());
260    else if (svd.m_computeThinU) m_workspace.resize(svd.cols());
261  }
262
263  bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix)
264  {
265    if(matrix.rows() > matrix.cols())
266    {
267      m_qr.compute(matrix);
268      svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>();
269      if(svd.m_computeFullU) m_qr.householderQ().evalTo(svd.m_matrixU, m_workspace);
270      else if(svd.m_computeThinU)
271      {
272        svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols());
273        m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_workspace);
274      }
275      if(svd.computeV()) svd.m_matrixV.setIdentity(matrix.cols(), matrix.cols());
276      return true;
277    }
278    return false;
279  }
280private:
281  HouseholderQR<MatrixType> m_qr;
282  typename internal::plain_col_type<MatrixType>::type m_workspace;
283};
284
285template<typename MatrixType>
286class qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true>
287{
288public:
289  typedef typename MatrixType::Index Index;
290  typedef typename MatrixType::Scalar Scalar;
291  enum
292  {
293    RowsAtCompileTime = MatrixType::RowsAtCompileTime,
294    ColsAtCompileTime = MatrixType::ColsAtCompileTime,
295    MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
296    MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
297    Options = MatrixType::Options
298  };
299
300  typedef Matrix<Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime>
301          TransposeTypeWithSameStorageOrder;
302
303  void allocate(const JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd)
304  {
305    if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols())
306    {
307      m_qr = HouseholderQR<TransposeTypeWithSameStorageOrder>(svd.cols(), svd.rows());
308    }
309    if (svd.m_computeFullV) m_workspace.resize(svd.cols());
310    else if (svd.m_computeThinV) m_workspace.resize(svd.rows());
311    m_adjoint.resize(svd.cols(), svd.rows());
312  }
313
314  bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix)
315  {
316    if(matrix.cols() > matrix.rows())
317    {
318      m_adjoint = matrix.adjoint();
319      m_qr.compute(m_adjoint);
320
321      svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint();
322      if(svd.m_computeFullV) m_qr.householderQ().evalTo(svd.m_matrixV, m_workspace);
323      else if(svd.m_computeThinV)
324      {
325        svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows());
326        m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_workspace);
327      }
328      if(svd.computeU()) svd.m_matrixU.setIdentity(matrix.rows(), matrix.rows());
329      return true;
330    }
331    else return false;
332  }
333
334private:
335  HouseholderQR<TransposeTypeWithSameStorageOrder> m_qr;
336  TransposeTypeWithSameStorageOrder m_adjoint;
337  typename internal::plain_row_type<MatrixType>::type m_workspace;
338};
339
340/*** 2x2 SVD implementation
341 ***
342 *** JacobiSVD consists in performing a series of 2x2 SVD subproblems
343 ***/
344
345template<typename MatrixType, int QRPreconditioner>
346struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, false>
347{
348  typedef JacobiSVD<MatrixType, QRPreconditioner> SVD;
349  typedef typename SVD::Index Index;
350  static void run(typename SVD::WorkMatrixType&, SVD&, Index, Index) {}
351};
352
353template<typename MatrixType, int QRPreconditioner>
354struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, true>
355{
356  typedef JacobiSVD<MatrixType, QRPreconditioner> SVD;
357  typedef typename MatrixType::Scalar Scalar;
358  typedef typename MatrixType::RealScalar RealScalar;
359  typedef typename SVD::Index Index;
360  static void run(typename SVD::WorkMatrixType& work_matrix, SVD& svd, Index p, Index q)
361  {
362    Scalar z;
363    JacobiRotation<Scalar> rot;
364    RealScalar n = sqrt(abs2(work_matrix.coeff(p,p)) + abs2(work_matrix.coeff(q,p)));
365    if(n==0)
366    {
367      z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q);
368      work_matrix.row(p) *= z;
369      if(svd.computeU()) svd.m_matrixU.col(p) *= conj(z);
370      z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q);
371      work_matrix.row(q) *= z;
372      if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z);
373    }
374    else
375    {
376      rot.c() = conj(work_matrix.coeff(p,p)) / n;
377      rot.s() = work_matrix.coeff(q,p) / n;
378      work_matrix.applyOnTheLeft(p,q,rot);
379      if(svd.computeU()) svd.m_matrixU.applyOnTheRight(p,q,rot.adjoint());
380      if(work_matrix.coeff(p,q) != Scalar(0))
381      {
382        Scalar z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q);
383        work_matrix.col(q) *= z;
384        if(svd.computeV()) svd.m_matrixV.col(q) *= z;
385      }
386      if(work_matrix.coeff(q,q) != Scalar(0))
387      {
388        z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q);
389        work_matrix.row(q) *= z;
390        if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z);
391      }
392    }
393  }
394};
395
396template<typename MatrixType, typename RealScalar, typename Index>
397void real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q,
398                            JacobiRotation<RealScalar> *j_left,
399                            JacobiRotation<RealScalar> *j_right)
400{
401  Matrix<RealScalar,2,2> m;
402  m << real(matrix.coeff(p,p)), real(matrix.coeff(p,q)),
403       real(matrix.coeff(q,p)), real(matrix.coeff(q,q));
404  JacobiRotation<RealScalar> rot1;
405  RealScalar t = m.coeff(0,0) + m.coeff(1,1);
406  RealScalar d = m.coeff(1,0) - m.coeff(0,1);
407  if(t == RealScalar(0))
408  {
409    rot1.c() = RealScalar(0);
410    rot1.s() = d > RealScalar(0) ? RealScalar(1) : RealScalar(-1);
411  }
412  else
413  {
414    RealScalar u = d / t;
415    rot1.c() = RealScalar(1) / sqrt(RealScalar(1) + abs2(u));
416    rot1.s() = rot1.c() * u;
417  }
418  m.applyOnTheLeft(0,1,rot1);
419  j_right->makeJacobi(m,0,1);
420  *j_left  = rot1 * j_right->transpose();
421}
422
423} // end namespace internal
424
425/** \ingroup SVD_Module
426  *
427  *
428  * \class JacobiSVD
429  *
430  * \brief Two-sided Jacobi SVD decomposition of a rectangular matrix
431  *
432  * \param MatrixType the type of the matrix of which we are computing the SVD decomposition
433  * \param QRPreconditioner this optional parameter allows to specify the type of QR decomposition that will be used internally
434  *                        for the R-SVD step for non-square matrices. See discussion of possible values below.
435  *
436  * SVD decomposition consists in decomposing any n-by-p matrix \a A as a product
437  *   \f[ A = U S V^* \f]
438  * where \a U is a n-by-n unitary, \a V is a p-by-p unitary, and \a S is a n-by-p real positive matrix which is zero outside of its main diagonal;
439  * the diagonal entries of S are known as the \em singular \em values of \a A and the columns of \a U and \a V are known as the left
440  * and right \em singular \em vectors of \a A respectively.
441  *
442  * Singular values are always sorted in decreasing order.
443  *
444  * This JacobiSVD decomposition computes only the singular values by default. If you want \a U or \a V, you need to ask for them explicitly.
445  *
446  * You can ask for only \em thin \a U or \a V to be computed, meaning the following. In case of a rectangular n-by-p matrix, letting \a m be the
447  * smaller value among \a n and \a p, there are only \a m singular vectors; the remaining columns of \a U and \a V do not correspond to actual
448  * singular vectors. Asking for \em thin \a U or \a V means asking for only their \a m first columns to be formed. So \a U is then a n-by-m matrix,
449  * and \a V is then a p-by-m matrix. Notice that thin \a U and \a V are all you need for (least squares) solving.
450  *
451  * Here's an example demonstrating basic usage:
452  * \include JacobiSVD_basic.cpp
453  * Output: \verbinclude JacobiSVD_basic.out
454  *
455  * This JacobiSVD class is a two-sided Jacobi R-SVD decomposition, ensuring optimal reliability and accuracy. The downside is that it's slower than
456  * bidiagonalizing SVD algorithms for large square matrices; however its complexity is still \f$ O(n^2p) \f$ where \a n is the smaller dimension and
457  * \a p is the greater dimension, meaning that it is still of the same order of complexity as the faster bidiagonalizing R-SVD algorithms.
458  * In particular, like any R-SVD, it takes advantage of non-squareness in that its complexity is only linear in the greater dimension.
459  *
460  * If the input matrix has inf or nan coefficients, the result of the computation is undefined, but the computation is guaranteed to
461  * terminate in finite (and reasonable) time.
462  *
463  * The possible values for QRPreconditioner are:
464  * \li ColPivHouseholderQRPreconditioner is the default. In practice it's very safe. It uses column-pivoting QR.
465  * \li FullPivHouseholderQRPreconditioner, is the safest and slowest. It uses full-pivoting QR.
466  *     Contrary to other QRs, it doesn't allow computing thin unitaries.
467  * \li HouseholderQRPreconditioner is the fastest, and less safe and accurate than the pivoting variants. It uses non-pivoting QR.
468  *     This is very similar in safety and accuracy to the bidiagonalization process used by bidiagonalizing SVD algorithms (since bidiagonalization
469  *     is inherently non-pivoting). However the resulting SVD is still more reliable than bidiagonalizing SVDs because the Jacobi-based iterarive
470  *     process is more reliable than the optimized bidiagonal SVD iterations.
471  * \li NoQRPreconditioner allows not to use a QR preconditioner at all. This is useful if you know that you will only be computing
472  *     JacobiSVD decompositions of square matrices. Non-square matrices require a QR preconditioner. Using this option will result in
473  *     faster compilation and smaller executable code. It won't significantly speed up computation, since JacobiSVD is always checking
474  *     if QR preconditioning is needed before applying it anyway.
475  *
476  * \sa MatrixBase::jacobiSvd()
477  */
478template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
479{
480  public:
481
482    typedef _MatrixType MatrixType;
483    typedef typename MatrixType::Scalar Scalar;
484    typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
485    typedef typename MatrixType::Index Index;
486    enum {
487      RowsAtCompileTime = MatrixType::RowsAtCompileTime,
488      ColsAtCompileTime = MatrixType::ColsAtCompileTime,
489      DiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_DYNAMIC(RowsAtCompileTime,ColsAtCompileTime),
490      MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
491      MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
492      MaxDiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(MaxRowsAtCompileTime,MaxColsAtCompileTime),
493      MatrixOptions = MatrixType::Options
494    };
495
496    typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime,
497                   MatrixOptions, MaxRowsAtCompileTime, MaxRowsAtCompileTime>
498            MatrixUType;
499    typedef Matrix<Scalar, ColsAtCompileTime, ColsAtCompileTime,
500                   MatrixOptions, MaxColsAtCompileTime, MaxColsAtCompileTime>
501            MatrixVType;
502    typedef typename internal::plain_diag_type<MatrixType, RealScalar>::type SingularValuesType;
503    typedef typename internal::plain_row_type<MatrixType>::type RowType;
504    typedef typename internal::plain_col_type<MatrixType>::type ColType;
505    typedef Matrix<Scalar, DiagSizeAtCompileTime, DiagSizeAtCompileTime,
506                   MatrixOptions, MaxDiagSizeAtCompileTime, MaxDiagSizeAtCompileTime>
507            WorkMatrixType;
508
509    /** \brief Default Constructor.
510      *
511      * The default constructor is useful in cases in which the user intends to
512      * perform decompositions via JacobiSVD::compute(const MatrixType&).
513      */
514    JacobiSVD()
515      : m_isInitialized(false),
516        m_isAllocated(false),
517        m_computationOptions(0),
518        m_rows(-1), m_cols(-1)
519    {}
520
521
522    /** \brief Default Constructor with memory preallocation
523      *
524      * Like the default constructor but with preallocation of the internal data
525      * according to the specified problem size.
526      * \sa JacobiSVD()
527      */
528    JacobiSVD(Index rows, Index cols, unsigned int computationOptions = 0)
529      : m_isInitialized(false),
530        m_isAllocated(false),
531        m_computationOptions(0),
532        m_rows(-1), m_cols(-1)
533    {
534      allocate(rows, cols, computationOptions);
535    }
536
537    /** \brief Constructor performing the decomposition of given matrix.
538     *
539     * \param matrix the matrix to decompose
540     * \param computationOptions optional parameter allowing to specify if you want full or thin U or V unitaries to be computed.
541     *                           By default, none is computed. This is a bit-field, the possible bits are #ComputeFullU, #ComputeThinU,
542     *                           #ComputeFullV, #ComputeThinV.
543     *
544     * Thin unitaries are only available if your matrix type has a Dynamic number of columns (for example MatrixXf). They also are not
545     * available with the (non-default) FullPivHouseholderQR preconditioner.
546     */
547    JacobiSVD(const MatrixType& matrix, unsigned int computationOptions = 0)
548      : m_isInitialized(false),
549        m_isAllocated(false),
550        m_computationOptions(0),
551        m_rows(-1), m_cols(-1)
552    {
553      compute(matrix, computationOptions);
554    }
555
556    /** \brief Method performing the decomposition of given matrix using custom options.
557     *
558     * \param matrix the matrix to decompose
559     * \param computationOptions optional parameter allowing to specify if you want full or thin U or V unitaries to be computed.
560     *                           By default, none is computed. This is a bit-field, the possible bits are #ComputeFullU, #ComputeThinU,
561     *                           #ComputeFullV, #ComputeThinV.
562     *
563     * Thin unitaries are only available if your matrix type has a Dynamic number of columns (for example MatrixXf). They also are not
564     * available with the (non-default) FullPivHouseholderQR preconditioner.
565     */
566    JacobiSVD& compute(const MatrixType& matrix, unsigned int computationOptions);
567
568    /** \brief Method performing the decomposition of given matrix using current options.
569     *
570     * \param matrix the matrix to decompose
571     *
572     * This method uses the current \a computationOptions, as already passed to the constructor or to compute(const MatrixType&, unsigned int).
573     */
574    JacobiSVD& compute(const MatrixType& matrix)
575    {
576      return compute(matrix, m_computationOptions);
577    }
578
579    /** \returns the \a U matrix.
580     *
581     * For the SVD decomposition of a n-by-p matrix, letting \a m be the minimum of \a n and \a p,
582     * the U matrix is n-by-n if you asked for #ComputeFullU, and is n-by-m if you asked for #ComputeThinU.
583     *
584     * The \a m first columns of \a U are the left singular vectors of the matrix being decomposed.
585     *
586     * This method asserts that you asked for \a U to be computed.
587     */
588    const MatrixUType& matrixU() const
589    {
590      eigen_assert(m_isInitialized && "JacobiSVD is not initialized.");
591      eigen_assert(computeU() && "This JacobiSVD decomposition didn't compute U. Did you ask for it?");
592      return m_matrixU;
593    }
594
595    /** \returns the \a V matrix.
596     *
597     * For the SVD decomposition of a n-by-p matrix, letting \a m be the minimum of \a n and \a p,
598     * the V matrix is p-by-p if you asked for #ComputeFullV, and is p-by-m if you asked for ComputeThinV.
599     *
600     * The \a m first columns of \a V are the right singular vectors of the matrix being decomposed.
601     *
602     * This method asserts that you asked for \a V to be computed.
603     */
604    const MatrixVType& matrixV() const
605    {
606      eigen_assert(m_isInitialized && "JacobiSVD is not initialized.");
607      eigen_assert(computeV() && "This JacobiSVD decomposition didn't compute V. Did you ask for it?");
608      return m_matrixV;
609    }
610
611    /** \returns the vector of singular values.
612     *
613     * For the SVD decomposition of a n-by-p matrix, letting \a m be the minimum of \a n and \a p, the
614     * returned vector has size \a m.  Singular values are always sorted in decreasing order.
615     */
616    const SingularValuesType& singularValues() const
617    {
618      eigen_assert(m_isInitialized && "JacobiSVD is not initialized.");
619      return m_singularValues;
620    }
621
622    /** \returns true if \a U (full or thin) is asked for in this SVD decomposition */
623    inline bool computeU() const { return m_computeFullU || m_computeThinU; }
624    /** \returns true if \a V (full or thin) is asked for in this SVD decomposition */
625    inline bool computeV() const { return m_computeFullV || m_computeThinV; }
626
627    /** \returns a (least squares) solution of \f$ A x = b \f$ using the current SVD decomposition of A.
628      *
629      * \param b the right-hand-side of the equation to solve.
630      *
631      * \note Solving requires both U and V to be computed. Thin U and V are enough, there is no need for full U or V.
632      *
633      * \note SVD solving is implicitly least-squares. Thus, this method serves both purposes of exact solving and least-squares solving.
634      * In other words, the returned solution is guaranteed to minimize the Euclidean norm \f$ \Vert A x - b \Vert \f$.
635      */
636    template<typename Rhs>
637    inline const internal::solve_retval<JacobiSVD, Rhs>
638    solve(const MatrixBase<Rhs>& b) const
639    {
640      eigen_assert(m_isInitialized && "JacobiSVD is not initialized.");
641      eigen_assert(computeU() && computeV() && "JacobiSVD::solve() requires both unitaries U and V to be computed (thin unitaries suffice).");
642      return internal::solve_retval<JacobiSVD, Rhs>(*this, b.derived());
643    }
644
645    /** \returns the number of singular values that are not exactly 0 */
646    Index nonzeroSingularValues() const
647    {
648      eigen_assert(m_isInitialized && "JacobiSVD is not initialized.");
649      return m_nonzeroSingularValues;
650    }
651
652    inline Index rows() const { return m_rows; }
653    inline Index cols() const { return m_cols; }
654
655  private:
656    void allocate(Index rows, Index cols, unsigned int computationOptions);
657
658  protected:
659    MatrixUType m_matrixU;
660    MatrixVType m_matrixV;
661    SingularValuesType m_singularValues;
662    WorkMatrixType m_workMatrix;
663    bool m_isInitialized, m_isAllocated;
664    bool m_computeFullU, m_computeThinU;
665    bool m_computeFullV, m_computeThinV;
666    unsigned int m_computationOptions;
667    Index m_nonzeroSingularValues, m_rows, m_cols, m_diagSize;
668
669    template<typename __MatrixType, int _QRPreconditioner, bool _IsComplex>
670    friend struct internal::svd_precondition_2x2_block_to_be_real;
671    template<typename __MatrixType, int _QRPreconditioner, int _Case, bool _DoAnything>
672    friend struct internal::qr_preconditioner_impl;
673
674    internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreColsThanRows> m_qr_precond_morecols;
675    internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreRowsThanCols> m_qr_precond_morerows;
676};
677
678template<typename MatrixType, int QRPreconditioner>
679void JacobiSVD<MatrixType, QRPreconditioner>::allocate(Index rows, Index cols, unsigned int computationOptions)
680{
681  eigen_assert(rows >= 0 && cols >= 0);
682
683  if (m_isAllocated &&
684      rows == m_rows &&
685      cols == m_cols &&
686      computationOptions == m_computationOptions)
687  {
688    return;
689  }
690
691  m_rows = rows;
692  m_cols = cols;
693  m_isInitialized = false;
694  m_isAllocated = true;
695  m_computationOptions = computationOptions;
696  m_computeFullU = (computationOptions & ComputeFullU) != 0;
697  m_computeThinU = (computationOptions & ComputeThinU) != 0;
698  m_computeFullV = (computationOptions & ComputeFullV) != 0;
699  m_computeThinV = (computationOptions & ComputeThinV) != 0;
700  eigen_assert(!(m_computeFullU && m_computeThinU) && "JacobiSVD: you can't ask for both full and thin U");
701  eigen_assert(!(m_computeFullV && m_computeThinV) && "JacobiSVD: you can't ask for both full and thin V");
702  eigen_assert(EIGEN_IMPLIES(m_computeThinU || m_computeThinV, MatrixType::ColsAtCompileTime==Dynamic) &&
703              "JacobiSVD: thin U and V are only available when your matrix has a dynamic number of columns.");
704  if (QRPreconditioner == FullPivHouseholderQRPreconditioner)
705  {
706      eigen_assert(!(m_computeThinU || m_computeThinV) &&
707              "JacobiSVD: can't compute thin U or thin V with the FullPivHouseholderQR preconditioner. "
708              "Use the ColPivHouseholderQR preconditioner instead.");
709  }
710  m_diagSize = (std::min)(m_rows, m_cols);
711  m_singularValues.resize(m_diagSize);
712  m_matrixU.resize(m_rows, m_computeFullU ? m_rows
713                          : m_computeThinU ? m_diagSize
714                          : 0);
715  m_matrixV.resize(m_cols, m_computeFullV ? m_cols
716                          : m_computeThinV ? m_diagSize
717                          : 0);
718  m_workMatrix.resize(m_diagSize, m_diagSize);
719
720  if(m_cols>m_rows) m_qr_precond_morecols.allocate(*this);
721  if(m_rows>m_cols) m_qr_precond_morerows.allocate(*this);
722}
723
724template<typename MatrixType, int QRPreconditioner>
725JacobiSVD<MatrixType, QRPreconditioner>&
726JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsigned int computationOptions)
727{
728  allocate(matrix.rows(), matrix.cols(), computationOptions);
729
730  // currently we stop when we reach precision 2*epsilon as the last bit of precision can require an unreasonable number of iterations,
731  // only worsening the precision of U and V as we accumulate more rotations
732  const RealScalar precision = RealScalar(2) * NumTraits<Scalar>::epsilon();
733
734  // limit for very small denormal numbers to be considered zero in order to avoid infinite loops (see bug 286)
735  const RealScalar considerAsZero = RealScalar(2) * std::numeric_limits<RealScalar>::denorm_min();
736
737  /*** step 1. The R-SVD step: we use a QR decomposition to reduce to the case of a square matrix */
738
739  if(!m_qr_precond_morecols.run(*this, matrix) && !m_qr_precond_morerows.run(*this, matrix))
740  {
741    m_workMatrix = matrix.block(0,0,m_diagSize,m_diagSize);
742    if(m_computeFullU) m_matrixU.setIdentity(m_rows,m_rows);
743    if(m_computeThinU) m_matrixU.setIdentity(m_rows,m_diagSize);
744    if(m_computeFullV) m_matrixV.setIdentity(m_cols,m_cols);
745    if(m_computeThinV) m_matrixV.setIdentity(m_cols, m_diagSize);
746  }
747
748  /*** step 2. The main Jacobi SVD iteration. ***/
749
750  bool finished = false;
751  while(!finished)
752  {
753    finished = true;
754
755    // do a sweep: for all index pairs (p,q), perform SVD of the corresponding 2x2 sub-matrix
756
757    for(Index p = 1; p < m_diagSize; ++p)
758    {
759      for(Index q = 0; q < p; ++q)
760      {
761        // if this 2x2 sub-matrix is not diagonal already...
762        // notice that this comparison will evaluate to false if any NaN is involved, ensuring that NaN's don't
763        // keep us iterating forever. Similarly, small denormal numbers are considered zero.
764        using std::max;
765        RealScalar threshold = (max)(considerAsZero, precision * (max)(internal::abs(m_workMatrix.coeff(p,p)),
766                                                                       internal::abs(m_workMatrix.coeff(q,q))));
767        if((max)(internal::abs(m_workMatrix.coeff(p,q)),internal::abs(m_workMatrix.coeff(q,p))) > threshold)
768        {
769          finished = false;
770
771          // perform SVD decomposition of 2x2 sub-matrix corresponding to indices p,q to make it diagonal
772          internal::svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner>::run(m_workMatrix, *this, p, q);
773          JacobiRotation<RealScalar> j_left, j_right;
774          internal::real_2x2_jacobi_svd(m_workMatrix, p, q, &j_left, &j_right);
775
776          // accumulate resulting Jacobi rotations
777          m_workMatrix.applyOnTheLeft(p,q,j_left);
778          if(computeU()) m_matrixU.applyOnTheRight(p,q,j_left.transpose());
779
780          m_workMatrix.applyOnTheRight(p,q,j_right);
781          if(computeV()) m_matrixV.applyOnTheRight(p,q,j_right);
782        }
783      }
784    }
785  }
786
787  /*** step 3. The work matrix is now diagonal, so ensure it's positive so its diagonal entries are the singular values ***/
788
789  for(Index i = 0; i < m_diagSize; ++i)
790  {
791    RealScalar a = internal::abs(m_workMatrix.coeff(i,i));
792    m_singularValues.coeffRef(i) = a;
793    if(computeU() && (a!=RealScalar(0))) m_matrixU.col(i) *= m_workMatrix.coeff(i,i)/a;
794  }
795
796  /*** step 4. Sort singular values in descending order and compute the number of nonzero singular values ***/
797
798  m_nonzeroSingularValues = m_diagSize;
799  for(Index i = 0; i < m_diagSize; i++)
800  {
801    Index pos;
802    RealScalar maxRemainingSingularValue = m_singularValues.tail(m_diagSize-i).maxCoeff(&pos);
803    if(maxRemainingSingularValue == RealScalar(0))
804    {
805      m_nonzeroSingularValues = i;
806      break;
807    }
808    if(pos)
809    {
810      pos += i;
811      std::swap(m_singularValues.coeffRef(i), m_singularValues.coeffRef(pos));
812      if(computeU()) m_matrixU.col(pos).swap(m_matrixU.col(i));
813      if(computeV()) m_matrixV.col(pos).swap(m_matrixV.col(i));
814    }
815  }
816
817  m_isInitialized = true;
818  return *this;
819}
820
821namespace internal {
822template<typename _MatrixType, int QRPreconditioner, typename Rhs>
823struct solve_retval<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs>
824  : solve_retval_base<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs>
825{
826  typedef JacobiSVD<_MatrixType, QRPreconditioner> JacobiSVDType;
827  EIGEN_MAKE_SOLVE_HELPERS(JacobiSVDType,Rhs)
828
829  template<typename Dest> void evalTo(Dest& dst) const
830  {
831    eigen_assert(rhs().rows() == dec().rows());
832
833    // A = U S V^*
834    // So A^{-1} = V S^{-1} U^*
835
836    Index diagSize = (std::min)(dec().rows(), dec().cols());
837    typename JacobiSVDType::SingularValuesType invertedSingVals(diagSize);
838
839    Index nonzeroSingVals = dec().nonzeroSingularValues();
840    invertedSingVals.head(nonzeroSingVals) = dec().singularValues().head(nonzeroSingVals).array().inverse();
841    invertedSingVals.tail(diagSize - nonzeroSingVals).setZero();
842
843    dst = dec().matrixV().leftCols(diagSize)
844        * invertedSingVals.asDiagonal()
845        * dec().matrixU().leftCols(diagSize).adjoint()
846        * rhs();
847  }
848};
849} // end namespace internal
850
851/** \svd_module
852  *
853  * \return the singular value decomposition of \c *this computed by two-sided
854  * Jacobi transformations.
855  *
856  * \sa class JacobiSVD
857  */
858template<typename Derived>
859JacobiSVD<typename MatrixBase<Derived>::PlainObject>
860MatrixBase<Derived>::jacobiSvd(unsigned int computationOptions) const
861{
862  return JacobiSVD<PlainObject>(*this, computationOptions);
863}
864
865} // end namespace Eigen
866
867#endif // EIGEN_JACOBISVD_H
868