diff --git a/Eigen/src/Core/util/Constants.h b/Eigen/src/Core/util/Constants.h index 60fdcf2e2..aacde6465 100644 --- a/Eigen/src/Core/util/Constants.h +++ b/Eigen/src/Core/util/Constants.h @@ -209,15 +209,6 @@ enum { OnTheRight = 2 }; -// options for SVD decomposition -enum { - SkipU = 0x1, - SkipV = 0x2, - AtLeastAsManyRowsAsCols = 0x4, - AtLeastAsManyColsAsRows = 0x8, - Square = AtLeastAsManyRowsAsCols | AtLeastAsManyColsAsRows -}; - /* the following could as well be written: * enum NoChange_t { NoChange }; * but it feels dangerous to disambiguate overloaded functions on enum/integer types. @@ -252,7 +243,7 @@ enum DecompositionOptions { Pivoting = 0x01, // LDLT, NoPivoting = 0x02, // LDLT, ComputeU = 0x10, // SVD, - ComputeR = 0x20, // SVD, + ComputeV = 0x20, // SVD, EigenvaluesOnly = 0x40, // all eigen solvers ComputeEigenvectors = 0x80, // all eigen solvers EigVecMask = EigenvaluesOnly | ComputeEigenvectors, @@ -262,6 +253,13 @@ enum DecompositionOptions { GenEigMask = Ax_lBx | ABx_lx | BAx_lx }; +enum QRPreconditioners { + NoQRPreconditioner, + HouseholderQRPreconditioner, + ColPivHouseholderQRPreconditioner, + FullPivHouseholderQRPreconditioner +}; + /** \brief Enum for reporting the status of a computation. */ enum ComputationInfo { diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index 3f67a8f4f..e69957a31 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h @@ -170,7 +170,7 @@ template class HouseholderQR; template class ColPivHouseholderQR; template class FullPivHouseholderQR; template class SVD; -template class JacobiSVD; +template class JacobiSVD; template class LLT; template class LDLT; template class HouseholderSequence; diff --git a/Eigen/src/SVD/JacobiSVD.h b/Eigen/src/SVD/JacobiSVD.h index 8ee0269b1..60b4c6353 100644 --- a/Eigen/src/SVD/JacobiSVD.h +++ b/Eigen/src/SVD/JacobiSVD.h @@ -26,22 +26,167 @@ #define EIGEN_JACOBISVD_H // forward declarations (needed by ICC) -// the empty bodies are required by VC -template::IsComplex> +// the empty bodies are required by MSVC +template::IsComplex> struct ei_svd_precondition_2x2_block_to_be_real {}; -templateMatrixType::ColsAtCompileTime))> +template MatrixType::ColsAtCompileTime) > struct ei_svd_precondition_if_more_rows_than_cols; -templateMatrixType::RowsAtCompileTime))> +template MatrixType::RowsAtCompileTime) > struct ei_svd_precondition_if_more_cols_than_rows; + +/*** QR preconditioners (R-SVD) ***/ + +enum { PreconditionIfMoreColsThanRows, PreconditionIfMoreRowsThanCols }; + +template +struct ei_qr_preconditioner_should_do_anything +{ + enum { a = MatrixType::RowsAtCompileTime != Dynamic && + MatrixType::ColsAtCompileTime != Dynamic && + MatrixType::ColsAtCompileTime <= MatrixType::RowsAtCompileTime, + b = MatrixType::RowsAtCompileTime != Dynamic && + MatrixType::ColsAtCompileTime != Dynamic && + MatrixType::RowsAtCompileTime <= MatrixType::ColsAtCompileTime, + ret = !( (QRPreconditioner == NoQRPreconditioner) || + (Case == PreconditionIfMoreColsThanRows && bool(a)) || + (Case == PreconditionIfMoreRowsThanCols && bool(b)) ) + }; +}; + +template::ret +> struct ei_qr_preconditioner_impl {}; + +template +struct ei_qr_preconditioner_impl +{ + static bool run(JacobiSVD&, const MatrixType&) + { + return false; + } +}; + +template +struct ei_qr_preconditioner_impl +{ + static bool run(JacobiSVD& svd, const MatrixType& matrix) + { + if(matrix.rows() > matrix.cols()) + { + FullPivHouseholderQR qr(matrix); + svd.m_workMatrix = qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView(); + if(svd.m_computeU) svd.m_matrixU = qr.matrixQ(); + if(svd.m_computeV) svd.m_matrixV = qr.colsPermutation(); + return true; + } + return false; + } +}; + +template +struct ei_qr_preconditioner_impl +{ + static bool run(JacobiSVD& svd, const MatrixType& matrix) + { + if(matrix.cols() > matrix.rows()) + { + typedef Matrix + TransposeTypeWithSameStorageOrder; + FullPivHouseholderQR qr(matrix.adjoint()); + svd.m_workMatrix = qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView().adjoint(); + if(svd.m_computeV) svd.m_matrixV = qr.matrixQ(); + if(svd.m_computeU) svd.m_matrixU = qr.colsPermutation(); + return true; + } + else return false; + } +}; + +template +struct ei_qr_preconditioner_impl +{ + static bool run(JacobiSVD& svd, const MatrixType& matrix) + { + if(matrix.rows() > matrix.cols()) + { + ColPivHouseholderQR qr(matrix); + svd.m_workMatrix = qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView(); + if(svd.m_computeU) svd.m_matrixU = qr.householderQ(); + if(svd.m_computeV) svd.m_matrixV = qr.colsPermutation(); + return true; + } + return false; + } +}; + +template +struct ei_qr_preconditioner_impl +{ + static bool run(JacobiSVD& svd, const MatrixType& matrix) + { + if(matrix.cols() > matrix.rows()) + { + typedef Matrix + TransposeTypeWithSameStorageOrder; + ColPivHouseholderQR qr(matrix.adjoint()); + svd.m_workMatrix = qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView().adjoint(); + if(svd.m_computeV) svd.m_matrixV = qr.householderQ(); + if(svd.m_computeU) svd.m_matrixU = qr.colsPermutation(); + return true; + } + else return false; + } +}; + +template +struct ei_qr_preconditioner_impl +{ + static bool run(JacobiSVD& svd, const MatrixType& matrix) + { + if(matrix.rows() > matrix.cols()) + { + HouseholderQR qr(matrix); + svd.m_workMatrix = qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView(); + if(svd.m_computeU) svd.m_matrixU = qr.householderQ(); + if(svd.m_computeV) svd.m_matrixV.setIdentity(matrix.cols(), matrix.cols()); + return true; + } + return false; + } +}; + +template +struct ei_qr_preconditioner_impl +{ + static bool run(JacobiSVD& svd, const MatrixType& matrix) + { + if(matrix.cols() > matrix.rows()) + { + typedef Matrix + TransposeTypeWithSameStorageOrder; + HouseholderQR qr(matrix.adjoint()); + svd.m_workMatrix = qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView().adjoint(); + if(svd.m_computeV) svd.m_matrixV = qr.householderQ(); + if(svd.m_computeU) svd.m_matrixU.setIdentity(matrix.rows(), matrix.rows()); + return true; + } + else return false; + } +}; + + + /** \ingroup SVD_Module * * @@ -50,23 +195,23 @@ struct ei_svd_precondition_if_more_cols_than_rows; * \brief Jacobi SVD decomposition of a square matrix * * \param MatrixType the type of the matrix of which we are computing the SVD decomposition - * \param Options a bit field of flags offering the following options: \c SkipU and \c SkipV allow to skip the computation of - * the unitaries \a U and \a V respectively; \c AtLeastAsManyRowsAsCols and \c AtLeastAsManyRowsAsCols allows - * to hint the compiler to only generate the corresponding code paths; \c Square is equivantent to the combination of - * the latter two bits and is useful when you know that the matrix is square. Note that when this information can - * be automatically deduced from the matrix type (e.g. a Matrix3f is always square), Eigen does it for you. + * \param QRPreconditioner the type of QR decomposition that will be used internally for the R-SVD step + * for non-square matrices. The default, FullPivHouseholderQR, is safest but slow. + * Consider using ColPivHouseholderQR instead of greater speed while still being + * quite safe, or even HouseholderQR to get closer to the speed and unsafety of + * bidiagonalizing SVD implementations. Finally, if you don't need to handle non-square matrices, + * you don't need any QR decomposition; you can then pass the dummy type NoQRDecomposition, + * which will result in smaller executable size and shorter compilation times. * * \sa MatrixBase::jacobiSvd() */ -template class JacobiSVD +template class JacobiSVD { private: typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits::Real RealScalar; typedef typename MatrixType::Index Index; enum { - ComputeU = (Options & SkipU) == 0, - ComputeV = (Options & SkipV) == 0, RowsAtCompileTime = MatrixType::RowsAtCompileTime, ColsAtCompileTime = MatrixType::ColsAtCompileTime, DiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_DYNAMIC(RowsAtCompileTime,ColsAtCompileTime), @@ -76,21 +221,18 @@ template class JacobiSVD MatrixOptions = MatrixType::Options }; - typedef Matrix DummyMatrixType; - typedef typename ei_meta_if, - DummyMatrixType>::ret MatrixUType; - typedef typename ei_meta_if, - DummyMatrixType>::ret MatrixVType; + typedef Matrix + MatrixUType; + typedef Matrix + MatrixVType; typedef typename ei_plain_diag_type::type SingularValuesType; typedef typename ei_plain_row_type::type RowType; typedef typename ei_plain_col_type::type ColType; typedef Matrix - WorkMatrixType; + WorkMatrixType; public: @@ -114,19 +256,20 @@ template class JacobiSVD m_workMatrix(rows, cols), m_isInitialized(false) {} - JacobiSVD(const MatrixType& matrix) : m_matrixU(matrix.rows(), matrix.rows()), - m_matrixV(matrix.cols(), matrix.cols()), - m_singularValues(), - m_workMatrix(), - m_isInitialized(false) + JacobiSVD(const MatrixType& matrix, unsigned int computationOptions = 0) + : m_matrixU(matrix.rows(), matrix.rows()), + m_matrixV(matrix.cols(), matrix.cols()), + m_singularValues(), + m_workMatrix(), + m_isInitialized(false) { const Index minSize = std::min(matrix.rows(), matrix.cols()); m_singularValues.resize(minSize); m_workMatrix.resize(minSize, minSize); - compute(matrix); + compute(matrix, computationOptions); } - JacobiSVD& compute(const MatrixType& matrix); + JacobiSVD& compute(const MatrixType& matrix, unsigned int computationOptions = 0); const MatrixUType& matrixU() const { @@ -151,33 +294,30 @@ template class JacobiSVD MatrixVType m_matrixV; SingularValuesType m_singularValues; WorkMatrixType m_workMatrix; - bool m_isInitialized; + bool m_isInitialized, m_computeU, m_computeV; - template + template friend struct ei_svd_precondition_2x2_block_to_be_real; - template - friend struct ei_svd_precondition_if_more_rows_than_cols; - template - friend struct ei_svd_precondition_if_more_cols_than_rows; + template + friend struct ei_qr_preconditioner_impl; }; -template -struct ei_svd_precondition_2x2_block_to_be_real +template +struct ei_svd_precondition_2x2_block_to_be_real { - typedef JacobiSVD SVD; + typedef JacobiSVD SVD; typedef typename SVD::Index Index; - static void run(typename SVD::WorkMatrixType&, JacobiSVD&, Index, Index) {} + static void run(typename SVD::WorkMatrixType&, SVD&, Index, Index) {} }; -template -struct ei_svd_precondition_2x2_block_to_be_real +template +struct ei_svd_precondition_2x2_block_to_be_real { - typedef JacobiSVD SVD; + typedef JacobiSVD SVD; typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; typedef typename SVD::Index Index; - enum { ComputeU = SVD::ComputeU, ComputeV = SVD::ComputeV }; - static void run(typename SVD::WorkMatrixType& work_matrix, JacobiSVD& svd, Index p, Index q) + static void run(typename SVD::WorkMatrixType& work_matrix, SVD& svd, Index p, Index q) { Scalar z; PlanarRotation rot; @@ -186,28 +326,28 @@ struct ei_svd_precondition_2x2_block_to_be_real { z = ei_abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q); work_matrix.row(p) *= z; - if(ComputeU) svd.m_matrixU.col(p) *= ei_conj(z); + if(svd.m_computeU) svd.m_matrixU.col(p) *= ei_conj(z); z = ei_abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q); work_matrix.row(q) *= z; - if(ComputeU) svd.m_matrixU.col(q) *= ei_conj(z); + if(svd.m_computeU) svd.m_matrixU.col(q) *= ei_conj(z); } else { rot.c() = ei_conj(work_matrix.coeff(p,p)) / n; rot.s() = work_matrix.coeff(q,p) / n; work_matrix.applyOnTheLeft(p,q,rot); - if(ComputeU) svd.m_matrixU.applyOnTheRight(p,q,rot.adjoint()); + if(svd.m_computeU) svd.m_matrixU.applyOnTheRight(p,q,rot.adjoint()); if(work_matrix.coeff(p,q) != Scalar(0)) { Scalar z = ei_abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q); work_matrix.col(q) *= z; - if(ComputeV) svd.m_matrixV.col(q) *= z; + if(svd.m_computeV) svd.m_matrixV.col(q) *= z; } if(work_matrix.coeff(q,q) != Scalar(0)) { z = ei_abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q); work_matrix.row(q) *= z; - if(ComputeU) svd.m_matrixU.col(q) *= ei_conj(z); + if(svd.m_computeU) svd.m_matrixU.col(q) *= ei_conj(z); } } } @@ -240,98 +380,24 @@ void ei_real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q, *j_left = rot1 * j_right->transpose(); } -template -struct ei_svd_precondition_if_more_rows_than_cols -{ - typedef JacobiSVD SVD; - static bool run(const MatrixType&, typename SVD::WorkMatrixType&, JacobiSVD&) { return false; } -}; - -template -struct ei_svd_precondition_if_more_rows_than_cols -{ - typedef JacobiSVD SVD; - typedef typename MatrixType::Scalar Scalar; - typedef typename MatrixType::RealScalar RealScalar; - typedef typename MatrixType::Index Index; - enum { ComputeU = SVD::ComputeU, ComputeV = SVD::ComputeV }; - static bool run(const MatrixType& matrix, typename SVD::WorkMatrixType& work_matrix, SVD& svd) - { - Index rows = matrix.rows(); - Index cols = matrix.cols(); - Index diagSize = cols; - if(rows > cols) - { - FullPivHouseholderQR qr(matrix); - work_matrix = qr.matrixQR().block(0,0,diagSize,diagSize).template triangularView(); - if(ComputeU) svd.m_matrixU = qr.matrixQ(); - if(ComputeV) svd.m_matrixV = qr.colsPermutation(); - - return true; - } - else return false; - } -}; - -template -struct ei_svd_precondition_if_more_cols_than_rows -{ - typedef JacobiSVD SVD; - static bool run(const MatrixType&, typename SVD::WorkMatrixType&, JacobiSVD&) { return false; } -}; - -template -struct ei_svd_precondition_if_more_cols_than_rows -{ - typedef JacobiSVD SVD; - typedef typename MatrixType::Scalar Scalar; - typedef typename MatrixType::RealScalar RealScalar; - typedef typename MatrixType::Index Index; - enum { - ComputeU = SVD::ComputeU, - ComputeV = SVD::ComputeV, - RowsAtCompileTime = SVD::RowsAtCompileTime, - ColsAtCompileTime = SVD::ColsAtCompileTime, - MaxRowsAtCompileTime = SVD::MaxRowsAtCompileTime, - MaxColsAtCompileTime = SVD::MaxColsAtCompileTime, - MatrixOptions = SVD::MatrixOptions - }; - - static bool run(const MatrixType& matrix, typename SVD::WorkMatrixType& work_matrix, SVD& svd) - { - Index rows = matrix.rows(); - Index cols = matrix.cols(); - Index diagSize = rows; - if(cols > rows) - { - typedef Matrix - TransposeTypeWithSameStorageOrder; - FullPivHouseholderQR qr(matrix.adjoint()); - work_matrix = qr.matrixQR().block(0,0,diagSize,diagSize).template triangularView().adjoint(); - if(ComputeV) svd.m_matrixV = qr.matrixQ(); - if(ComputeU) svd.m_matrixU = qr.colsPermutation(); - return true; - } - else return false; - } -}; - -template -JacobiSVD& JacobiSVD::compute(const MatrixType& matrix) +template +JacobiSVD& +JacobiSVD::compute(const MatrixType& matrix, unsigned int computationOptions) { + m_computeU = computationOptions & ComputeU; + m_computeV = computationOptions & ComputeV; Index rows = matrix.rows(); Index cols = matrix.cols(); Index diagSize = std::min(rows, cols); m_singularValues.resize(diagSize); const RealScalar precision = 2 * NumTraits::epsilon(); - if(!ei_svd_precondition_if_more_rows_than_cols::run(matrix, m_workMatrix, *this) - && !ei_svd_precondition_if_more_cols_than_rows::run(matrix, m_workMatrix, *this)) + if(!ei_qr_preconditioner_impl::run(*this, matrix) + && !ei_qr_preconditioner_impl::run(*this, matrix)) { m_workMatrix = matrix.block(0,0,diagSize,diagSize); - if(ComputeU) m_matrixU.setIdentity(rows,rows); - if(ComputeV) m_matrixV.setIdentity(cols,cols); + if(m_computeU) m_matrixU.setIdentity(rows,rows); + if(m_computeV) m_matrixV.setIdentity(cols,cols); } bool finished = false; @@ -346,16 +412,16 @@ JacobiSVD& JacobiSVD::compute(const Ma > std::max(ei_abs(m_workMatrix.coeff(p,p)),ei_abs(m_workMatrix.coeff(q,q)))*precision) { finished = false; - ei_svd_precondition_2x2_block_to_be_real::run(m_workMatrix, *this, p, q); + ei_svd_precondition_2x2_block_to_be_real::run(m_workMatrix, *this, p, q); PlanarRotation j_left, j_right; ei_real_2x2_jacobi_svd(m_workMatrix, p, q, &j_left, &j_right); m_workMatrix.applyOnTheLeft(p,q,j_left); - if(ComputeU) m_matrixU.applyOnTheRight(p,q,j_left.transpose()); + if(m_computeU) m_matrixU.applyOnTheRight(p,q,j_left.transpose()); m_workMatrix.applyOnTheRight(p,q,j_right); - if(ComputeV) m_matrixV.applyOnTheRight(p,q,j_right); + if(m_computeV) m_matrixV.applyOnTheRight(p,q,j_right); } } } @@ -365,7 +431,7 @@ JacobiSVD& JacobiSVD::compute(const Ma { RealScalar a = ei_abs(m_workMatrix.coeff(i,i)); m_singularValues.coeffRef(i) = a; - if(ComputeU && (a!=RealScalar(0))) m_matrixU.col(i) *= m_workMatrix.coeff(i,i)/a; + if(m_computeU && (a!=RealScalar(0))) m_matrixU.col(i) *= m_workMatrix.coeff(i,i)/a; } for(Index i = 0; i < diagSize; i++) @@ -376,8 +442,8 @@ JacobiSVD& JacobiSVD::compute(const Ma { pos += i; std::swap(m_singularValues.coeffRef(i), m_singularValues.coeffRef(pos)); - if(ComputeU) m_matrixU.col(pos).swap(m_matrixU.col(i)); - if(ComputeV) m_matrixV.col(pos).swap(m_matrixV.col(i)); + if(m_computeU) m_matrixU.col(pos).swap(m_matrixU.col(i)); + if(m_computeV) m_matrixV.col(pos).swap(m_matrixV.col(i)); } } diff --git a/test/jacobisvd.cpp b/test/jacobisvd.cpp index 24dbc22a5..56b91c983 100644 --- a/test/jacobisvd.cpp +++ b/test/jacobisvd.cpp @@ -27,7 +27,8 @@ #include #include -template void svd(const MatrixType& m = MatrixType(), bool pickrandom = true) +template +void svd_with_qr_preconditioner(const MatrixType& m = MatrixType(), bool pickrandom = true) { typedef typename MatrixType::Index Index; Index rows = m.rows(); @@ -49,7 +50,7 @@ template void svd(const MatrixType& m if(pickrandom) a = MatrixType::Random(rows,cols); else a = m; - JacobiSVD svd(a); + JacobiSVD svd(a, ComputeU|ComputeV); MatrixType sigma = MatrixType::Zero(rows,cols); sigma.diagonal() = svd.singularValues().template cast(); MatrixUType u = svd.matrixU(); @@ -63,11 +64,19 @@ template void svd(const MatrixType& m VERIFY_IS_UNITARY(v); } +template +void svd(const MatrixType& m = MatrixType(), bool pickrandom = true) +{ + svd_with_qr_preconditioner(m, pickrandom); + svd_with_qr_preconditioner(m, pickrandom); + svd_with_qr_preconditioner(m, pickrandom); +} + template void svd_verify_assert() { MatrixType tmp; - SVD svd; + JacobiSVD svd; //VERIFY_RAISES_ASSERT(svd.solve(tmp, &tmp)) VERIFY_RAISES_ASSERT(svd.matrixU()) VERIFY_RAISES_ASSERT(svd.singularValues()) @@ -84,24 +93,24 @@ void test_jacobisvd() Matrix2cd m; m << 0, 1, 0, 1; - CALL_SUBTEST_1(( svd(m, false) )); + CALL_SUBTEST_1(( svd(m, false) )); m << 1, 0, 1, 0; - CALL_SUBTEST_1(( svd(m, false) )); + CALL_SUBTEST_1(( svd(m, false) )); Matrix2d n; n << 1, 1, 1, -1; - CALL_SUBTEST_2(( svd(n, false) )); - CALL_SUBTEST_3(( svd() )); - CALL_SUBTEST_4(( svd() )); - CALL_SUBTEST_5(( svd , AtLeastAsManyColsAsRows>() )); - CALL_SUBTEST_6(( svd , AtLeastAsManyRowsAsCols>(Matrix(10,2)) )); + CALL_SUBTEST_2(( svd(n, false) )); + CALL_SUBTEST_3(( svd() )); + CALL_SUBTEST_4(( svd() )); + CALL_SUBTEST_5(( svd >() )); + CALL_SUBTEST_6(( svd >(Matrix(10,2)) )); - CALL_SUBTEST_7(( svd(MatrixXf(50,50)) )); - CALL_SUBTEST_8(( svd(MatrixXcd(14,7)) )); + CALL_SUBTEST_7(( svd(MatrixXf(50,50)) )); + CALL_SUBTEST_8(( svd(MatrixXcd(14,7)) )); } - CALL_SUBTEST_9(( svd(MatrixXf(300,200)) )); - CALL_SUBTEST_10(( svd(MatrixXcd(100,150)) )); + CALL_SUBTEST_9(( svd(MatrixXf(300,200)) )); + CALL_SUBTEST_10(( svd(MatrixXcd(100,150)) )); CALL_SUBTEST_3(( svd_verify_assert() )); CALL_SUBTEST_3(( svd_verify_assert() ));