Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-04-19 09:06:45 

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2009-2010 Benoit Jacob <jacob.benoit.1@gmail.com>
0005 // Copyright (C) 2013-2014 Gael Guennebaud <gael.guennebaud@inria.fr>
0006 //
0007 // This Source Code Form is subject to the terms of the Mozilla
0008 // Public License v. 2.0. If a copy of the MPL was not distributed
0009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0010 
0011 #ifndef EIGEN_JACOBISVD_H
0012 #define EIGEN_JACOBISVD_H
0013 
0014 namespace RivetEigen { 
0015 
0016 namespace internal {
0017 // forward declaration (needed by ICC)
0018 // the empty body is required by MSVC
0019 template<typename MatrixType, int QRPreconditioner,
0020          bool IsComplex = NumTraits<typename MatrixType::Scalar>::IsComplex>
0021 struct svd_precondition_2x2_block_to_be_real {};
0022 
0023 /*** QR preconditioners (R-SVD)
0024  ***
0025  *** Their role is to reduce the problem of computing the SVD to the case of a square matrix.
0026  *** This approach, known as R-SVD, is an optimization for rectangular-enough matrices, and is a requirement for
0027  *** JacobiSVD which by itself is only able to work on square matrices.
0028  ***/
0029 
0030 enum { PreconditionIfMoreColsThanRows, PreconditionIfMoreRowsThanCols };
0031 
0032 template<typename MatrixType, int QRPreconditioner, int Case>
0033 struct qr_preconditioner_should_do_anything
0034 {
0035   enum { a = MatrixType::RowsAtCompileTime != Dynamic &&
0036              MatrixType::ColsAtCompileTime != Dynamic &&
0037              MatrixType::ColsAtCompileTime <= MatrixType::RowsAtCompileTime,
0038          b = MatrixType::RowsAtCompileTime != Dynamic &&
0039              MatrixType::ColsAtCompileTime != Dynamic &&
0040              MatrixType::RowsAtCompileTime <= MatrixType::ColsAtCompileTime,
0041          ret = !( (QRPreconditioner == NoQRPreconditioner) ||
0042                   (Case == PreconditionIfMoreColsThanRows && bool(a)) ||
0043                   (Case == PreconditionIfMoreRowsThanCols && bool(b)) )
0044   };
0045 };
0046 
0047 template<typename MatrixType, int QRPreconditioner, int Case,
0048          bool DoAnything = qr_preconditioner_should_do_anything<MatrixType, QRPreconditioner, Case>::ret
0049 > struct qr_preconditioner_impl {};
0050 
0051 template<typename MatrixType, int QRPreconditioner, int Case>
0052 class qr_preconditioner_impl<MatrixType, QRPreconditioner, Case, false>
0053 {
0054 public:
0055   void allocate(const JacobiSVD<MatrixType, QRPreconditioner>&) {}
0056   bool run(JacobiSVD<MatrixType, QRPreconditioner>&, const MatrixType&)
0057   {
0058     return false;
0059   }
0060 };
0061 
0062 /*** preconditioner using FullPivHouseholderQR ***/
0063 
0064 template<typename MatrixType>
0065 class qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
0066 {
0067 public:
0068   typedef typename MatrixType::Scalar Scalar;
0069   enum
0070   {
0071     RowsAtCompileTime = MatrixType::RowsAtCompileTime,
0072     MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime
0073   };
0074   typedef Matrix<Scalar, 1, RowsAtCompileTime, RowMajor, 1, MaxRowsAtCompileTime> WorkspaceType;
0075 
0076   void allocate(const JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd)
0077   {
0078     if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols())
0079     {
0080       m_qr.~QRType();
0081       ::new (&m_qr) QRType(svd.rows(), svd.cols());
0082     }
0083     if (svd.m_computeFullU) m_workspace.resize(svd.rows());
0084   }
0085 
0086   bool run(JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
0087   {
0088     if(matrix.rows() > matrix.cols())
0089     {
0090       m_qr.compute(matrix);
0091       svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>();
0092       if(svd.m_computeFullU) m_qr.matrixQ().evalTo(svd.m_matrixU, m_workspace);
0093       if(svd.computeV()) svd.m_matrixV = m_qr.colsPermutation();
0094       return true;
0095     }
0096     return false;
0097   }
0098 private:
0099   typedef FullPivHouseholderQR<MatrixType> QRType;
0100   QRType m_qr;
0101   WorkspaceType m_workspace;
0102 };
0103 
0104 template<typename MatrixType>
0105 class qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true>
0106 {
0107 public:
0108   typedef typename MatrixType::Scalar Scalar;
0109   enum
0110   {
0111     RowsAtCompileTime = MatrixType::RowsAtCompileTime,
0112     ColsAtCompileTime = MatrixType::ColsAtCompileTime,
0113     MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
0114     MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
0115     Options = MatrixType::Options
0116   };
0117 
0118   typedef typename internal::make_proper_matrix_type<
0119     Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime
0120   >::type TransposeTypeWithSameStorageOrder;
0121 
0122   void allocate(const JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd)
0123   {
0124     if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols())
0125     {
0126       m_qr.~QRType();
0127       ::new (&m_qr) QRType(svd.cols(), svd.rows());
0128     }
0129     m_adjoint.resize(svd.cols(), svd.rows());
0130     if (svd.m_computeFullV) m_workspace.resize(svd.cols());
0131   }
0132 
0133   bool run(JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
0134   {
0135     if(matrix.cols() > matrix.rows())
0136     {
0137       m_adjoint = matrix.adjoint();
0138       m_qr.compute(m_adjoint);
0139       svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint();
0140       if(svd.m_computeFullV) m_qr.matrixQ().evalTo(svd.m_matrixV, m_workspace);
0141       if(svd.computeU()) svd.m_matrixU = m_qr.colsPermutation();
0142       return true;
0143     }
0144     else return false;
0145   }
0146 private:
0147   typedef FullPivHouseholderQR<TransposeTypeWithSameStorageOrder> QRType;
0148   QRType m_qr;
0149   TransposeTypeWithSameStorageOrder m_adjoint;
0150   typename internal::plain_row_type<MatrixType>::type m_workspace;
0151 };
0152 
0153 /*** preconditioner using ColPivHouseholderQR ***/
0154 
0155 template<typename MatrixType>
0156 class qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
0157 {
0158 public:
0159   void allocate(const JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd)
0160   {
0161     if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols())
0162     {
0163       m_qr.~QRType();
0164       ::new (&m_qr) QRType(svd.rows(), svd.cols());
0165     }
0166     if (svd.m_computeFullU) m_workspace.resize(svd.rows());
0167     else if (svd.m_computeThinU) m_workspace.resize(svd.cols());
0168   }
0169 
0170   bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
0171   {
0172     if(matrix.rows() > matrix.cols())
0173     {
0174       m_qr.compute(matrix);
0175       svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>();
0176       if(svd.m_computeFullU) m_qr.householderQ().evalTo(svd.m_matrixU, m_workspace);
0177       else if(svd.m_computeThinU)
0178       {
0179         svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols());
0180         m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_workspace);
0181       }
0182       if(svd.computeV()) svd.m_matrixV = m_qr.colsPermutation();
0183       return true;
0184     }
0185     return false;
0186   }
0187 
0188 private:
0189   typedef ColPivHouseholderQR<MatrixType> QRType;
0190   QRType m_qr;
0191   typename internal::plain_col_type<MatrixType>::type m_workspace;
0192 };
0193 
0194 template<typename MatrixType>
0195 class qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true>
0196 {
0197 public:
0198   typedef typename MatrixType::Scalar Scalar;
0199   enum
0200   {
0201     RowsAtCompileTime = MatrixType::RowsAtCompileTime,
0202     ColsAtCompileTime = MatrixType::ColsAtCompileTime,
0203     MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
0204     MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
0205     Options = MatrixType::Options
0206   };
0207 
0208   typedef typename internal::make_proper_matrix_type<
0209     Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime
0210   >::type TransposeTypeWithSameStorageOrder;
0211 
0212   void allocate(const JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd)
0213   {
0214     if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols())
0215     {
0216       m_qr.~QRType();
0217       ::new (&m_qr) QRType(svd.cols(), svd.rows());
0218     }
0219     if (svd.m_computeFullV) m_workspace.resize(svd.cols());
0220     else if (svd.m_computeThinV) m_workspace.resize(svd.rows());
0221     m_adjoint.resize(svd.cols(), svd.rows());
0222   }
0223 
0224   bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
0225   {
0226     if(matrix.cols() > matrix.rows())
0227     {
0228       m_adjoint = matrix.adjoint();
0229       m_qr.compute(m_adjoint);
0230 
0231       svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint();
0232       if(svd.m_computeFullV) m_qr.householderQ().evalTo(svd.m_matrixV, m_workspace);
0233       else if(svd.m_computeThinV)
0234       {
0235         svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows());
0236         m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_workspace);
0237       }
0238       if(svd.computeU()) svd.m_matrixU = m_qr.colsPermutation();
0239       return true;
0240     }
0241     else return false;
0242   }
0243 
0244 private:
0245   typedef ColPivHouseholderQR<TransposeTypeWithSameStorageOrder> QRType;
0246   QRType m_qr;
0247   TransposeTypeWithSameStorageOrder m_adjoint;
0248   typename internal::plain_row_type<MatrixType>::type m_workspace;
0249 };
0250 
0251 /*** preconditioner using HouseholderQR ***/
0252 
0253 template<typename MatrixType>
0254 class qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
0255 {
0256 public:
0257   void allocate(const JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd)
0258   {
0259     if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols())
0260     {
0261       m_qr.~QRType();
0262       ::new (&m_qr) QRType(svd.rows(), svd.cols());
0263     }
0264     if (svd.m_computeFullU) m_workspace.resize(svd.rows());
0265     else if (svd.m_computeThinU) m_workspace.resize(svd.cols());
0266   }
0267 
0268   bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix)
0269   {
0270     if(matrix.rows() > matrix.cols())
0271     {
0272       m_qr.compute(matrix);
0273       svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>();
0274       if(svd.m_computeFullU) m_qr.householderQ().evalTo(svd.m_matrixU, m_workspace);
0275       else if(svd.m_computeThinU)
0276       {
0277         svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols());
0278         m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_workspace);
0279       }
0280       if(svd.computeV()) svd.m_matrixV.setIdentity(matrix.cols(), matrix.cols());
0281       return true;
0282     }
0283     return false;
0284   }
0285 private:
0286   typedef HouseholderQR<MatrixType> QRType;
0287   QRType m_qr;
0288   typename internal::plain_col_type<MatrixType>::type m_workspace;
0289 };
0290 
0291 template<typename MatrixType>
0292 class qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true>
0293 {
0294 public:
0295   typedef typename MatrixType::Scalar Scalar;
0296   enum
0297   {
0298     RowsAtCompileTime = MatrixType::RowsAtCompileTime,
0299     ColsAtCompileTime = MatrixType::ColsAtCompileTime,
0300     MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
0301     MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
0302     Options = MatrixType::Options
0303   };
0304 
0305   typedef typename internal::make_proper_matrix_type<
0306     Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime
0307   >::type TransposeTypeWithSameStorageOrder;
0308 
0309   void allocate(const JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd)
0310   {
0311     if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols())
0312     {
0313       m_qr.~QRType();
0314       ::new (&m_qr) QRType(svd.cols(), svd.rows());
0315     }
0316     if (svd.m_computeFullV) m_workspace.resize(svd.cols());
0317     else if (svd.m_computeThinV) m_workspace.resize(svd.rows());
0318     m_adjoint.resize(svd.cols(), svd.rows());
0319   }
0320 
0321   bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix)
0322   {
0323     if(matrix.cols() > matrix.rows())
0324     {
0325       m_adjoint = matrix.adjoint();
0326       m_qr.compute(m_adjoint);
0327 
0328       svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint();
0329       if(svd.m_computeFullV) m_qr.householderQ().evalTo(svd.m_matrixV, m_workspace);
0330       else if(svd.m_computeThinV)
0331       {
0332         svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows());
0333         m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_workspace);
0334       }
0335       if(svd.computeU()) svd.m_matrixU.setIdentity(matrix.rows(), matrix.rows());
0336       return true;
0337     }
0338     else return false;
0339   }
0340 
0341 private:
0342   typedef HouseholderQR<TransposeTypeWithSameStorageOrder> QRType;
0343   QRType m_qr;
0344   TransposeTypeWithSameStorageOrder m_adjoint;
0345   typename internal::plain_row_type<MatrixType>::type m_workspace;
0346 };
0347 
0348 /*** 2x2 SVD implementation
0349  ***
0350  *** JacobiSVD consists in performing a series of 2x2 SVD subproblems
0351  ***/
0352 
0353 template<typename MatrixType, int QRPreconditioner>
0354 struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, false>
0355 {
0356   typedef JacobiSVD<MatrixType, QRPreconditioner> SVD;
0357   typedef typename MatrixType::RealScalar RealScalar;
0358   static bool run(typename SVD::WorkMatrixType&, SVD&, Index, Index, RealScalar&) { return true; }
0359 };
0360 
0361 template<typename MatrixType, int QRPreconditioner>
0362 struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, true>
0363 {
0364   typedef JacobiSVD<MatrixType, QRPreconditioner> SVD;
0365   typedef typename MatrixType::Scalar Scalar;
0366   typedef typename MatrixType::RealScalar RealScalar;
0367   static bool run(typename SVD::WorkMatrixType& work_matrix, SVD& svd, Index p, Index q, RealScalar& maxDiagEntry)
0368   {
0369     using std::sqrt;
0370     using std::abs;
0371     Scalar z;
0372     JacobiRotation<Scalar> rot;
0373     RealScalar n = sqrt(numext::abs2(work_matrix.coeff(p,p)) + numext::abs2(work_matrix.coeff(q,p)));
0374 
0375     const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)();
0376     const RealScalar precision = NumTraits<Scalar>::epsilon();
0377 
0378     if(n==0)
0379     {
0380       // make sure first column is zero
0381       work_matrix.coeffRef(p,p) = work_matrix.coeffRef(q,p) = Scalar(0);
0382 
0383       if(abs(numext::imag(work_matrix.coeff(p,q)))>considerAsZero)
0384       {
0385         // 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
0386         z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q);
0387         work_matrix.row(p) *= z;
0388         if(svd.computeU()) svd.m_matrixU.col(p) *= conj(z);
0389       }
0390       if(abs(numext::imag(work_matrix.coeff(q,q)))>considerAsZero)
0391       {
0392         z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q);
0393         work_matrix.row(q) *= z;
0394         if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z);
0395       }
0396       // otherwise the second row is already zero, so we have nothing to do.
0397     }
0398     else
0399     {
0400       rot.c() = conj(work_matrix.coeff(p,p)) / n;
0401       rot.s() = work_matrix.coeff(q,p) / n;
0402       work_matrix.applyOnTheLeft(p,q,rot);
0403       if(svd.computeU()) svd.m_matrixU.applyOnTheRight(p,q,rot.adjoint());
0404       if(abs(numext::imag(work_matrix.coeff(p,q)))>considerAsZero)
0405       {
0406         z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q);
0407         work_matrix.col(q) *= z;
0408         if(svd.computeV()) svd.m_matrixV.col(q) *= z;
0409       }
0410       if(abs(numext::imag(work_matrix.coeff(q,q)))>considerAsZero)
0411       {
0412         z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q);
0413         work_matrix.row(q) *= z;
0414         if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z);
0415       }
0416     }
0417 
0418     // update largest diagonal entry
0419     maxDiagEntry = numext::maxi<RealScalar>(maxDiagEntry,numext::maxi<RealScalar>(abs(work_matrix.coeff(p,p)), abs(work_matrix.coeff(q,q))));
0420     // and check whether the 2x2 block is already diagonal
0421     RealScalar threshold = numext::maxi<RealScalar>(considerAsZero, precision * maxDiagEntry);
0422     return abs(work_matrix.coeff(p,q))>threshold || abs(work_matrix.coeff(q,p)) > threshold;
0423   }
0424 };
0425 
0426 template<typename _MatrixType, int QRPreconditioner> 
0427 struct traits<JacobiSVD<_MatrixType,QRPreconditioner> >
0428         : traits<_MatrixType>
0429 {
0430   typedef _MatrixType MatrixType;
0431 };
0432 
0433 } // end namespace internal
0434 
0435 /** \ingroup SVD_Module
0436   *
0437   *
0438   * \class JacobiSVD
0439   *
0440   * \brief Two-sided Jacobi SVD decomposition of a rectangular matrix
0441   *
0442   * \tparam _MatrixType the type of the matrix of which we are computing the SVD decomposition
0443   * \tparam QRPreconditioner this optional parameter allows to specify the type of QR decomposition that will be used internally
0444   *                        for the R-SVD step for non-square matrices. See discussion of possible values below.
0445   *
0446   * SVD decomposition consists in decomposing any n-by-p matrix \a A as a product
0447   *   \f[ A = U S V^* \f]
0448   * 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;
0449   * 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
0450   * and right \em singular \em vectors of \a A respectively.
0451   *
0452   * Singular values are always sorted in decreasing order.
0453   *
0454   * 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.
0455   *
0456   * 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
0457   * 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
0458   * 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,
0459   * 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.
0460   *
0461   * Here's an example demonstrating basic usage:
0462   * \include JacobiSVD_basic.cpp
0463   * Output: \verbinclude JacobiSVD_basic.out
0464   *
0465   * This JacobiSVD class is a two-sided Jacobi R-SVD decomposition, ensuring optimal reliability and accuracy. The downside is that it's slower than
0466   * 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
0467   * \a p is the greater dimension, meaning that it is still of the same order of complexity as the faster bidiagonalizing R-SVD algorithms.
0468   * In particular, like any R-SVD, it takes advantage of non-squareness in that its complexity is only linear in the greater dimension.
0469   *
0470   * If the input matrix has inf or nan coefficients, the result of the computation is undefined, but the computation is guaranteed to
0471   * terminate in finite (and reasonable) time.
0472   *
0473   * The possible values for QRPreconditioner are:
0474   * \li ColPivHouseholderQRPreconditioner is the default. In practice it's very safe. It uses column-pivoting QR.
0475   * \li FullPivHouseholderQRPreconditioner, is the safest and slowest. It uses full-pivoting QR.
0476   *     Contrary to other QRs, it doesn't allow computing thin unitaries.
0477   * \li HouseholderQRPreconditioner is the fastest, and less safe and accurate than the pivoting variants. It uses non-pivoting QR.
0478   *     This is very similar in safety and accuracy to the bidiagonalization process used by bidiagonalizing SVD algorithms (since bidiagonalization
0479   *     is inherently non-pivoting). However the resulting SVD is still more reliable than bidiagonalizing SVDs because the Jacobi-based iterarive
0480   *     process is more reliable than the optimized bidiagonal SVD iterations.
0481   * \li NoQRPreconditioner allows not to use a QR preconditioner at all. This is useful if you know that you will only be computing
0482   *     JacobiSVD decompositions of square matrices. Non-square matrices require a QR preconditioner. Using this option will result in
0483   *     faster compilation and smaller executable code. It won't significantly speed up computation, since JacobiSVD is always checking
0484   *     if QR preconditioning is needed before applying it anyway.
0485   *
0486   * \sa MatrixBase::jacobiSvd()
0487   */
0488 template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
0489  : public SVDBase<JacobiSVD<_MatrixType,QRPreconditioner> >
0490 {
0491     typedef SVDBase<JacobiSVD> Base;
0492   public:
0493 
0494     typedef _MatrixType MatrixType;
0495     typedef typename MatrixType::Scalar Scalar;
0496     typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
0497     enum {
0498       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
0499       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
0500       DiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_DYNAMIC(RowsAtCompileTime,ColsAtCompileTime),
0501       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
0502       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
0503       MaxDiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(MaxRowsAtCompileTime,MaxColsAtCompileTime),
0504       MatrixOptions = MatrixType::Options
0505     };
0506 
0507     typedef typename Base::MatrixUType MatrixUType;
0508     typedef typename Base::MatrixVType MatrixVType;
0509     typedef typename Base::SingularValuesType SingularValuesType;
0510     
0511     typedef typename internal::plain_row_type<MatrixType>::type RowType;
0512     typedef typename internal::plain_col_type<MatrixType>::type ColType;
0513     typedef Matrix<Scalar, DiagSizeAtCompileTime, DiagSizeAtCompileTime,
0514                    MatrixOptions, MaxDiagSizeAtCompileTime, MaxDiagSizeAtCompileTime>
0515             WorkMatrixType;
0516 
0517     /** \brief Default Constructor.
0518       *
0519       * The default constructor is useful in cases in which the user intends to
0520       * perform decompositions via JacobiSVD::compute(const MatrixType&).
0521       */
0522     JacobiSVD()
0523     {}
0524 
0525 
0526     /** \brief Default Constructor with memory preallocation
0527       *
0528       * Like the default constructor but with preallocation of the internal data
0529       * according to the specified problem size.
0530       * \sa JacobiSVD()
0531       */
0532     JacobiSVD(Index rows, Index cols, unsigned int computationOptions = 0)
0533     {
0534       allocate(rows, cols, computationOptions);
0535     }
0536 
0537     /** \brief Constructor performing the decomposition of given matrix.
0538      *
0539      * \param matrix the matrix to decompose
0540      * \param computationOptions optional parameter allowing to specify if you want full or thin U or V unitaries to be computed.
0541      *                           By default, none is computed. This is a bit-field, the possible bits are #ComputeFullU, #ComputeThinU,
0542      *                           #ComputeFullV, #ComputeThinV.
0543      *
0544      * Thin unitaries are only available if your matrix type has a Dynamic number of columns (for example MatrixXf). They also are not
0545      * available with the (non-default) FullPivHouseholderQR preconditioner.
0546      */
0547     explicit JacobiSVD(const MatrixType& matrix, unsigned int computationOptions = 0)
0548     {
0549       compute(matrix, computationOptions);
0550     }
0551 
0552     /** \brief Method performing the decomposition of given matrix using custom options.
0553      *
0554      * \param matrix the matrix to decompose
0555      * \param computationOptions optional parameter allowing to specify if you want full or thin U or V unitaries to be computed.
0556      *                           By default, none is computed. This is a bit-field, the possible bits are #ComputeFullU, #ComputeThinU,
0557      *                           #ComputeFullV, #ComputeThinV.
0558      *
0559      * Thin unitaries are only available if your matrix type has a Dynamic number of columns (for example MatrixXf). They also are not
0560      * available with the (non-default) FullPivHouseholderQR preconditioner.
0561      */
0562     JacobiSVD& compute(const MatrixType& matrix, unsigned int computationOptions);
0563 
0564     /** \brief Method performing the decomposition of given matrix using current options.
0565      *
0566      * \param matrix the matrix to decompose
0567      *
0568      * This method uses the current \a computationOptions, as already passed to the constructor or to compute(const MatrixType&, unsigned int).
0569      */
0570     JacobiSVD& compute(const MatrixType& matrix)
0571     {
0572       return compute(matrix, m_computationOptions);
0573     }
0574 
0575     using Base::computeU;
0576     using Base::computeV;
0577     using Base::rows;
0578     using Base::cols;
0579     using Base::rank;
0580 
0581   private:
0582     void allocate(Index rows, Index cols, unsigned int computationOptions);
0583 
0584   protected:
0585     using Base::m_matrixU;
0586     using Base::m_matrixV;
0587     using Base::m_singularValues;
0588     using Base::m_info;
0589     using Base::m_isInitialized;
0590     using Base::m_isAllocated;
0591     using Base::m_usePrescribedThreshold;
0592     using Base::m_computeFullU;
0593     using Base::m_computeThinU;
0594     using Base::m_computeFullV;
0595     using Base::m_computeThinV;
0596     using Base::m_computationOptions;
0597     using Base::m_nonzeroSingularValues;
0598     using Base::m_rows;
0599     using Base::m_cols;
0600     using Base::m_diagSize;
0601     using Base::m_prescribedThreshold;
0602     WorkMatrixType m_workMatrix;
0603 
0604     template<typename __MatrixType, int _QRPreconditioner, bool _IsComplex>
0605     friend struct internal::svd_precondition_2x2_block_to_be_real;
0606     template<typename __MatrixType, int _QRPreconditioner, int _Case, bool _DoAnything>
0607     friend struct internal::qr_preconditioner_impl;
0608 
0609     internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreColsThanRows> m_qr_precond_morecols;
0610     internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreRowsThanCols> m_qr_precond_morerows;
0611     MatrixType m_scaledMatrix;
0612 };
0613 
0614 template<typename MatrixType, int QRPreconditioner>
0615 void JacobiSVD<MatrixType, QRPreconditioner>::allocate(RivetEigen::Index rows, RivetEigen::Index cols, unsigned int computationOptions)
0616 {
0617   eigen_assert(rows >= 0 && cols >= 0);
0618 
0619   if (m_isAllocated &&
0620       rows == m_rows &&
0621       cols == m_cols &&
0622       computationOptions == m_computationOptions)
0623   {
0624     return;
0625   }
0626 
0627   m_rows = rows;
0628   m_cols = cols;
0629   m_info = Success;
0630   m_isInitialized = false;
0631   m_isAllocated = true;
0632   m_computationOptions = computationOptions;
0633   m_computeFullU = (computationOptions & ComputeFullU) != 0;
0634   m_computeThinU = (computationOptions & ComputeThinU) != 0;
0635   m_computeFullV = (computationOptions & ComputeFullV) != 0;
0636   m_computeThinV = (computationOptions & ComputeThinV) != 0;
0637   eigen_assert(!(m_computeFullU && m_computeThinU) && "JacobiSVD: you can't ask for both full and thin U");
0638   eigen_assert(!(m_computeFullV && m_computeThinV) && "JacobiSVD: you can't ask for both full and thin V");
0639   eigen_assert(EIGEN_IMPLIES(m_computeThinU || m_computeThinV, MatrixType::ColsAtCompileTime==Dynamic) &&
0640               "JacobiSVD: thin U and V are only available when your matrix has a dynamic number of columns.");
0641   if (QRPreconditioner == FullPivHouseholderQRPreconditioner)
0642   {
0643       eigen_assert(!(m_computeThinU || m_computeThinV) &&
0644               "JacobiSVD: can't compute thin U or thin V with the FullPivHouseholderQR preconditioner. "
0645               "Use the ColPivHouseholderQR preconditioner instead.");
0646   }
0647   m_diagSize = (std::min)(m_rows, m_cols);
0648   m_singularValues.resize(m_diagSize);
0649   if(RowsAtCompileTime==Dynamic)
0650     m_matrixU.resize(m_rows, m_computeFullU ? m_rows
0651                             : m_computeThinU ? m_diagSize
0652                             : 0);
0653   if(ColsAtCompileTime==Dynamic)
0654     m_matrixV.resize(m_cols, m_computeFullV ? m_cols
0655                             : m_computeThinV ? m_diagSize
0656                             : 0);
0657   m_workMatrix.resize(m_diagSize, m_diagSize);
0658   
0659   if(m_cols>m_rows)   m_qr_precond_morecols.allocate(*this);
0660   if(m_rows>m_cols)   m_qr_precond_morerows.allocate(*this);
0661   if(m_rows!=m_cols)  m_scaledMatrix.resize(rows,cols);
0662 }
0663 
0664 template<typename MatrixType, int QRPreconditioner>
0665 JacobiSVD<MatrixType, QRPreconditioner>&
0666 JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsigned int computationOptions)
0667 {
0668   using std::abs;
0669   allocate(matrix.rows(), matrix.cols(), computationOptions);
0670 
0671   // currently we stop when we reach precision 2*epsilon as the last bit of precision can require an unreasonable number of iterations,
0672   // only worsening the precision of U and V as we accumulate more rotations
0673   const RealScalar precision = RealScalar(2) * NumTraits<Scalar>::epsilon();
0674 
0675   // limit for denormal numbers to be considered zero in order to avoid infinite loops (see bug 286)
0676   const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)();
0677 
0678   // Scaling factor to reduce over/under-flows
0679   RealScalar scale = matrix.cwiseAbs().template maxCoeff<PropagateNaN>();
0680   if (!(numext::isfinite)(scale)) {
0681     m_isInitialized = true;
0682     m_info = InvalidInput;
0683     return *this;
0684   }
0685   if(scale==RealScalar(0)) scale = RealScalar(1);
0686   
0687   /*** step 1. The R-SVD step: we use a QR decomposition to reduce to the case of a square matrix */
0688 
0689   if(m_rows!=m_cols)
0690   {
0691     m_scaledMatrix = matrix / scale;
0692     m_qr_precond_morecols.run(*this, m_scaledMatrix);
0693     m_qr_precond_morerows.run(*this, m_scaledMatrix);
0694   }
0695   else
0696   {
0697     m_workMatrix = matrix.block(0,0,m_diagSize,m_diagSize) / scale;
0698     if(m_computeFullU) m_matrixU.setIdentity(m_rows,m_rows);
0699     if(m_computeThinU) m_matrixU.setIdentity(m_rows,m_diagSize);
0700     if(m_computeFullV) m_matrixV.setIdentity(m_cols,m_cols);
0701     if(m_computeThinV) m_matrixV.setIdentity(m_cols, m_diagSize);
0702   }
0703 
0704   /*** step 2. The main Jacobi SVD iteration. ***/
0705   RealScalar maxDiagEntry = m_workMatrix.cwiseAbs().diagonal().maxCoeff();
0706 
0707   bool finished = false;
0708   while(!finished)
0709   {
0710     finished = true;
0711 
0712     // do a sweep: for all index pairs (p,q), perform SVD of the corresponding 2x2 sub-matrix
0713 
0714     for(Index p = 1; p < m_diagSize; ++p)
0715     {
0716       for(Index q = 0; q < p; ++q)
0717       {
0718         // if this 2x2 sub-matrix is not diagonal already...
0719         // notice that this comparison will evaluate to false if any NaN is involved, ensuring that NaN's don't
0720         // keep us iterating forever. Similarly, small denormal numbers are considered zero.
0721         RealScalar threshold = numext::maxi<RealScalar>(considerAsZero, precision * maxDiagEntry);
0722         if(abs(m_workMatrix.coeff(p,q))>threshold || abs(m_workMatrix.coeff(q,p)) > threshold)
0723         {
0724           finished = false;
0725           // perform SVD decomposition of 2x2 sub-matrix corresponding to indices p,q to make it diagonal
0726           // the complex to real operation returns true if the updated 2x2 block is not already diagonal
0727           if(internal::svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner>::run(m_workMatrix, *this, p, q, maxDiagEntry))
0728           {
0729             JacobiRotation<RealScalar> j_left, j_right;
0730             internal::real_2x2_jacobi_svd(m_workMatrix, p, q, &j_left, &j_right);
0731 
0732             // accumulate resulting Jacobi rotations
0733             m_workMatrix.applyOnTheLeft(p,q,j_left);
0734             if(computeU()) m_matrixU.applyOnTheRight(p,q,j_left.transpose());
0735 
0736             m_workMatrix.applyOnTheRight(p,q,j_right);
0737             if(computeV()) m_matrixV.applyOnTheRight(p,q,j_right);
0738 
0739             // keep track of the largest diagonal coefficient
0740             maxDiagEntry = numext::maxi<RealScalar>(maxDiagEntry,numext::maxi<RealScalar>(abs(m_workMatrix.coeff(p,p)), abs(m_workMatrix.coeff(q,q))));
0741           }
0742         }
0743       }
0744     }
0745   }
0746 
0747   /*** step 3. The work matrix is now diagonal, so ensure it's positive so its diagonal entries are the singular values ***/
0748 
0749   for(Index i = 0; i < m_diagSize; ++i)
0750   {
0751     // For a complex matrix, some diagonal coefficients might note have been
0752     // treated by svd_precondition_2x2_block_to_be_real, and the imaginary part
0753     // of some diagonal entry might not be null.
0754     if(NumTraits<Scalar>::IsComplex && abs(numext::imag(m_workMatrix.coeff(i,i)))>considerAsZero)
0755     {
0756       RealScalar a = abs(m_workMatrix.coeff(i,i));
0757       m_singularValues.coeffRef(i) = abs(a);
0758       if(computeU()) m_matrixU.col(i) *= m_workMatrix.coeff(i,i)/a;
0759     }
0760     else
0761     {
0762       // m_workMatrix.coeff(i,i) is already real, no difficulty:
0763       RealScalar a = numext::real(m_workMatrix.coeff(i,i));
0764       m_singularValues.coeffRef(i) = abs(a);
0765       if(computeU() && (a<RealScalar(0))) m_matrixU.col(i) = -m_matrixU.col(i);
0766     }
0767   }
0768   
0769   m_singularValues *= scale;
0770 
0771   /*** step 4. Sort singular values in descending order and compute the number of nonzero singular values ***/
0772 
0773   m_nonzeroSingularValues = m_diagSize;
0774   for(Index i = 0; i < m_diagSize; i++)
0775   {
0776     Index pos;
0777     RealScalar maxRemainingSingularValue = m_singularValues.tail(m_diagSize-i).maxCoeff(&pos);
0778     if(maxRemainingSingularValue == RealScalar(0))
0779     {
0780       m_nonzeroSingularValues = i;
0781       break;
0782     }
0783     if(pos)
0784     {
0785       pos += i;
0786       std::swap(m_singularValues.coeffRef(i), m_singularValues.coeffRef(pos));
0787       if(computeU()) m_matrixU.col(pos).swap(m_matrixU.col(i));
0788       if(computeV()) m_matrixV.col(pos).swap(m_matrixV.col(i));
0789     }
0790   }
0791 
0792   m_isInitialized = true;
0793   return *this;
0794 }
0795 
0796 /** \svd_module
0797   *
0798   * \return the singular value decomposition of \c *this computed by two-sided
0799   * Jacobi transformations.
0800   *
0801   * \sa class JacobiSVD
0802   */
0803 template<typename Derived>
0804 JacobiSVD<typename MatrixBase<Derived>::PlainObject>
0805 MatrixBase<Derived>::jacobiSvd(unsigned int computationOptions) const
0806 {
0807   return JacobiSVD<PlainObject>(*this, computationOptions);
0808 }
0809 
0810 } // end namespace RivetEigen
0811 
0812 #endif // EIGEN_JACOBISVD_H
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!