From 8ba8d90063ab6b297b1ba912f806eaa71cbdf003 Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Fri, 8 Oct 2010 10:42:40 -0400 Subject: [PATCH] add option to compute thin U/V. By default nothing is computed. You have to ask explicitly for thin/full U/V if you want them. --- Eigen/src/Core/util/Constants.h | 8 +- Eigen/src/SVD/JacobiSVD.h | 141 ++++++++++++++++++++++++-------- test/jacobisvd.cpp | 105 +++++++++++++++++------- 3 files changed, 186 insertions(+), 68 deletions(-) diff --git a/Eigen/src/Core/util/Constants.h b/Eigen/src/Core/util/Constants.h index aacde6465..60da5c76a 100644 --- a/Eigen/src/Core/util/Constants.h +++ b/Eigen/src/Core/util/Constants.h @@ -242,10 +242,12 @@ enum AccessorLevels { enum DecompositionOptions { Pivoting = 0x01, // LDLT, NoPivoting = 0x02, // LDLT, - ComputeU = 0x10, // SVD, - ComputeV = 0x20, // SVD, + ComputeFullU = 0x04, // SVD, + ComputeThinU = 0x08, // SVD, + ComputeFullV = 0x10, // SVD, + ComputeThinV = 0x20, // SVD, EigenvaluesOnly = 0x40, // all eigen solvers - ComputeEigenvectors = 0x80, // all eigen solvers + ComputeEigenvectors = 0x80, // all eigen solvers EigVecMask = EigenvaluesOnly | ComputeEigenvectors, Ax_lBx = 0x100, ABx_lx = 0x200, diff --git a/Eigen/src/SVD/JacobiSVD.h b/Eigen/src/SVD/JacobiSVD.h index 60b4c6353..bce7719bf 100644 --- a/Eigen/src/SVD/JacobiSVD.h +++ b/Eigen/src/SVD/JacobiSVD.h @@ -81,10 +81,12 @@ struct ei_qr_preconditioner_impl matrix.cols()) { + ei_assert(!svd.m_computeThinU && "JacobiSVD: can't compute a thin U with the FullPivHouseholderQR preconditioner. " + "Use the ColPivHouseholderQR preconditioner instead."); 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(); + if(svd.m_computeFullU) svd.m_matrixU = qr.matrixQ(); + if(svd.computeV()) svd.m_matrixV = qr.colsPermutation(); return true; } return false; @@ -98,13 +100,15 @@ struct ei_qr_preconditioner_impl matrix.rows()) { + ei_assert(!svd.m_computeThinV && "JacobiSVD: can't compute a thin V with the FullPivHouseholderQR preconditioner. " + "Use the ColPivHouseholderQR preconditioner instead."); 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(); + if(svd.m_computeFullV) svd.m_matrixV = qr.matrixQ(); + if(svd.computeU()) svd.m_matrixU = qr.colsPermutation(); return true; } else return false; @@ -120,8 +124,12 @@ struct ei_qr_preconditioner_impl 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(); + if(svd.m_computeFullU) svd.m_matrixU = qr.householderQ(); + else if(svd.m_computeThinU) { + svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols()); + qr.householderQ().applyThisOnTheLeft(svd.m_matrixU); + } + if(svd.computeV()) svd.m_matrixV = qr.colsPermutation(); return true; } return false; @@ -140,8 +148,12 @@ struct ei_qr_preconditioner_impl 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(); + if(svd.m_computeFullV) svd.m_matrixV = qr.householderQ(); + else if(svd.m_computeThinV) { + svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows()); + qr.householderQ().applyThisOnTheLeft(svd.m_matrixV); + } + if(svd.computeU()) svd.m_matrixU = qr.colsPermutation(); return true; } else return false; @@ -157,8 +169,12 @@ struct ei_qr_preconditioner_impl 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()); + if(svd.m_computeFullU) svd.m_matrixU = qr.householderQ(); + else if(svd.m_computeThinU) { + svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols()); + qr.householderQ().applyThisOnTheLeft(svd.m_matrixU); + } + if(svd.computeV()) svd.m_matrixV.setIdentity(matrix.cols(), matrix.cols()); return true; } return false; @@ -177,8 +193,12 @@ struct ei_qr_preconditioner_impl 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()); + if(svd.m_computeFullV) svd.m_matrixV = qr.householderQ(); + else if(svd.m_computeThinV) { + svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows()); + qr.householderQ().applyThisOnTheLeft(svd.m_matrixV); + } + if(svd.computeU()) svd.m_matrixU.setIdentity(matrix.rows(), matrix.rows()); return true; } else return false; @@ -195,13 +215,21 @@ struct ei_qr_preconditioner_impl class JacobiSVD m_workMatrix(rows, cols), m_isInitialized(false) {} + /** \brief Constructor performing the decomposition of given matrix. + * + * \param matrix the matrix to decompose + * \param computationOptions optional parameter allowing to specify if you want full or thin U or V unitaries to be computed. + * By default, none is computed. This is a bit-field, the possible bits are ComputeFullU, ComputeThinU, + * ComputeFullV, ComputeThinV. + * + * Thin unitaries are not available with the default FullPivHouseholderQRPreconditioner, see class documentation for details. + * If you want thin unitaries, use another preconditioner, for example: + * \code + * JacobiSVD svd(matrix, ComputeThinU); + * \endcode + * + * Thin unitaries also are only available if your matrix type has a Dynamic number of columns (for example MatrixXf). + */ JacobiSVD(const MatrixType& matrix, unsigned int computationOptions = 0) : m_matrixU(matrix.rows(), matrix.rows()), m_matrixV(matrix.cols(), matrix.cols()), @@ -269,11 +312,27 @@ template class JacobiSVD compute(matrix, computationOptions); } + /** \brief Method performing the decomposition of given matrix. + * + * \param matrix the matrix to decompose + * \param computationOptions optional parameter allowing to specify if you want full or thin U or V unitaries to be computed. + * By default, none is computed. This is a bit-field, the possible bits are ComputeFullU, ComputeThinU, + * ComputeFullV, ComputeThinV. + * + * Thin unitaries are not available with the default FullPivHouseholderQRPreconditioner, see class documentation for details. + * If you want thin unitaries, use another preconditioner, for example: + * \code + * JacobiSVD svd(matrix, ComputeThinU); + * \endcode + * + * Thin unitaries also are only available if your matrix type has a Dynamic number of columns (for example MatrixXf). + */ JacobiSVD& compute(const MatrixType& matrix, unsigned int computationOptions = 0); const MatrixUType& matrixU() const { ei_assert(m_isInitialized && "JacobiSVD is not initialized."); + ei_assert(computeU() && "This JacobiSVD decomposition didn't compute U. Did you ask for it?"); return m_matrixU; } @@ -286,15 +345,21 @@ template class JacobiSVD const MatrixVType& matrixV() const { ei_assert(m_isInitialized && "JacobiSVD is not initialized."); + ei_assert(computeV() && "This JacobiSVD decomposition didn't compute V. Did you ask for it?"); return m_matrixV; } + inline bool computeU() const { return m_computeFullU || m_computeThinU; } + inline bool computeV() const { return m_computeFullV || m_computeThinV; } + protected: MatrixUType m_matrixU; MatrixVType m_matrixV; SingularValuesType m_singularValues; WorkMatrixType m_workMatrix; - bool m_isInitialized, m_computeU, m_computeV; + bool m_isInitialized; + bool m_computeFullU, m_computeThinU; + bool m_computeFullV, m_computeThinV; template friend struct ei_svd_precondition_2x2_block_to_be_real; @@ -326,28 +391,28 @@ struct ei_svd_precondition_2x2_block_to_be_real JacobiSVD& JacobiSVD::compute(const MatrixType& matrix, unsigned int computationOptions) { - m_computeU = computationOptions & ComputeU; - m_computeV = computationOptions & ComputeV; + m_computeFullU = computationOptions & ComputeFullU; + m_computeThinU = computationOptions & ComputeThinU; + m_computeFullV = computationOptions & ComputeFullV; + m_computeThinV = computationOptions & ComputeThinV; + ei_assert(!(m_computeFullU && m_computeThinU) && "JacobiSVD: you can't ask for both full and thin U"); + ei_assert(!(m_computeFullV && m_computeThinV) && "JacobiSVD: you can't ask for both full and thin V"); + ei_assert(EIGEN_IMPLIES(m_computeThinU || m_computeThinV, MatrixType::ColsAtCompileTime==Dynamic) && + "JacobiSVD: thin U and V are only available when your matrix has a dynamic number of columns."); Index rows = matrix.rows(); Index cols = matrix.cols(); Index diagSize = std::min(rows, cols); @@ -396,8 +467,10 @@ JacobiSVD::compute(const MatrixType& matrix, unsig && !ei_qr_preconditioner_impl::run(*this, matrix)) { m_workMatrix = matrix.block(0,0,diagSize,diagSize); - if(m_computeU) m_matrixU.setIdentity(rows,rows); - if(m_computeV) m_matrixV.setIdentity(cols,cols); + if(m_computeFullU) m_matrixU.setIdentity(rows,rows); + if(m_computeThinU) m_matrixU.setIdentity(rows,diagSize); + if(m_computeFullV) m_matrixV.setIdentity(cols,cols); + if(m_computeThinV) m_matrixV.setIdentity(diagSize,cols); } bool finished = false; @@ -418,10 +491,10 @@ JacobiSVD::compute(const MatrixType& matrix, unsig ei_real_2x2_jacobi_svd(m_workMatrix, p, q, &j_left, &j_right); m_workMatrix.applyOnTheLeft(p,q,j_left); - if(m_computeU) m_matrixU.applyOnTheRight(p,q,j_left.transpose()); + if(computeU()) m_matrixU.applyOnTheRight(p,q,j_left.transpose()); m_workMatrix.applyOnTheRight(p,q,j_right); - if(m_computeV) m_matrixV.applyOnTheRight(p,q,j_right); + if(computeV()) m_matrixV.applyOnTheRight(p,q,j_right); } } } @@ -431,7 +504,7 @@ JacobiSVD::compute(const MatrixType& matrix, unsig { RealScalar a = ei_abs(m_workMatrix.coeff(i,i)); m_singularValues.coeffRef(i) = a; - if(m_computeU && (a!=RealScalar(0))) m_matrixU.col(i) *= m_workMatrix.coeff(i,i)/a; + if(computeU() && (a!=RealScalar(0))) m_matrixU.col(i) *= m_workMatrix.coeff(i,i)/a; } for(Index i = 0; i < diagSize; i++) @@ -442,8 +515,8 @@ JacobiSVD::compute(const MatrixType& matrix, unsig { pos += i; std::swap(m_singularValues.coeffRef(i), m_singularValues.coeffRef(pos)); - if(m_computeU) m_matrixU.col(pos).swap(m_matrixU.col(i)); - if(m_computeV) m_matrixV.col(pos).swap(m_matrixV.col(i)); + if(computeU()) m_matrixU.col(pos).swap(m_matrixU.col(i)); + if(computeV()) m_matrixV.col(pos).swap(m_matrixV.col(i)); } } diff --git a/test/jacobisvd.cpp b/test/jacobisvd.cpp index 56b91c983..32b25f9d2 100644 --- a/test/jacobisvd.cpp +++ b/test/jacobisvd.cpp @@ -28,7 +28,7 @@ #include template -void svd_with_qr_preconditioner(const MatrixType& m = MatrixType(), bool pickrandom = true) +void jacobisvd_check_full(const MatrixType& m, const JacobiSVD& svd) { typedef typename MatrixType::Index Index; Index rows = m.rows(); @@ -46,33 +46,76 @@ void svd_with_qr_preconditioner(const MatrixType& m = MatrixType(), bool pickran typedef Matrix ColVectorType; typedef Matrix InputVectorType; - MatrixType a; - if(pickrandom) a = MatrixType::Random(rows,cols); - else a = m; - - JacobiSVD svd(a, ComputeU|ComputeV); MatrixType sigma = MatrixType::Zero(rows,cols); sigma.diagonal() = svd.singularValues().template cast(); MatrixUType u = svd.matrixU(); MatrixVType v = svd.matrixV(); - //std::cout << "a\n" << a << std::endl; - //std::cout << "b\n" << u * sigma * v.adjoint() << std::endl; - - VERIFY_IS_APPROX(a, u * sigma * v.adjoint()); + VERIFY_IS_APPROX(m, u * sigma * v.adjoint()); VERIFY_IS_UNITARY(u); VERIFY_IS_UNITARY(v); } -template -void svd(const MatrixType& m = MatrixType(), bool pickrandom = true) +template +void jacobisvd_compare_to_full(const MatrixType& m, + unsigned int computationOptions, + const JacobiSVD& referenceSvd) { - svd_with_qr_preconditioner(m, pickrandom); - svd_with_qr_preconditioner(m, pickrandom); - svd_with_qr_preconditioner(m, pickrandom); + typedef typename MatrixType::Index Index; + Index rows = m.rows(); + Index cols = m.cols(); + Index diagSize = std::min(rows, cols); + + JacobiSVD svd(m, computationOptions); + + VERIFY_IS_EQUAL(svd.singularValues(), referenceSvd.singularValues()); + if(computationOptions & ComputeFullU) + VERIFY_IS_EQUAL(svd.matrixU(), referenceSvd.matrixU()); + if(computationOptions & ComputeThinU) + VERIFY_IS_EQUAL(svd.matrixU(), referenceSvd.matrixU().leftCols(diagSize)); + if(computationOptions & ComputeFullV) + VERIFY_IS_EQUAL(svd.matrixV(), referenceSvd.matrixV()); + if(computationOptions & ComputeThinV) + VERIFY_IS_EQUAL(svd.matrixV(), referenceSvd.matrixV().leftCols(diagSize)); } -template void svd_verify_assert() +template +void jacobisvd_test_all_computation_options(const MatrixType& m) +{ + if (QRPreconditioner == NoQRPreconditioner && m.rows() != m.cols()) + return; + JacobiSVD fullSvd(m, ComputeFullU|ComputeFullV); + + jacobisvd_check_full(m, fullSvd); + + if(QRPreconditioner == FullPivHouseholderQRPreconditioner) + return; + + jacobisvd_compare_to_full(m, ComputeFullU, fullSvd); + jacobisvd_compare_to_full(m, ComputeFullV, fullSvd); + jacobisvd_compare_to_full(m, 0, fullSvd); + + if (MatrixType::ColsAtCompileTime == Dynamic) { + // thin U/V are only available with dynamic number of columns + jacobisvd_compare_to_full(m, ComputeFullU|ComputeThinV, fullSvd); + jacobisvd_compare_to_full(m, ComputeThinV, fullSvd); + jacobisvd_compare_to_full(m, ComputeThinU|ComputeFullV, fullSvd); + jacobisvd_compare_to_full(m, ComputeThinU , fullSvd); + jacobisvd_compare_to_full(m, ComputeThinU|ComputeThinV, fullSvd); + } +} + +template +void jacobisvd(const MatrixType& a = MatrixType(), bool pickrandom = true) +{ + MatrixType m = pickrandom ? MatrixType::Random(a.rows(), a.cols()) : a; + jacobisvd_test_all_computation_options(m); + jacobisvd_test_all_computation_options(m); + jacobisvd_test_all_computation_options(m); + jacobisvd_test_all_computation_options(m); +} + +template void jacobisvd_verify_assert() { MatrixType tmp; @@ -93,29 +136,29 @@ void test_jacobisvd() Matrix2cd m; m << 0, 1, 0, 1; - CALL_SUBTEST_1(( svd(m, false) )); + CALL_SUBTEST_1(( jacobisvd(m, false) )); m << 1, 0, 1, 0; - CALL_SUBTEST_1(( svd(m, false) )); + CALL_SUBTEST_1(( jacobisvd(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 >() )); - CALL_SUBTEST_6(( svd >(Matrix(10,2)) )); + CALL_SUBTEST_2(( jacobisvd(n, false) )); + CALL_SUBTEST_3(( jacobisvd() )); + CALL_SUBTEST_4(( jacobisvd() )); + CALL_SUBTEST_5(( jacobisvd >() )); + CALL_SUBTEST_6(( jacobisvd >(Matrix(10,2)) )); - CALL_SUBTEST_7(( svd(MatrixXf(50,50)) )); - CALL_SUBTEST_8(( svd(MatrixXcd(14,7)) )); + CALL_SUBTEST_7(( jacobisvd(MatrixXf(50,50)) )); + CALL_SUBTEST_8(( jacobisvd(MatrixXcd(14,7)) )); } - CALL_SUBTEST_9(( svd(MatrixXf(300,200)) )); - CALL_SUBTEST_10(( svd(MatrixXcd(100,150)) )); + CALL_SUBTEST_9(( jacobisvd(MatrixXf(300,200)) )); + CALL_SUBTEST_10(( jacobisvd(MatrixXcd(100,150)) )); - CALL_SUBTEST_3(( svd_verify_assert() )); - CALL_SUBTEST_3(( svd_verify_assert() )); - CALL_SUBTEST_9(( svd_verify_assert() )); - CALL_SUBTEST_11(( svd_verify_assert() )); + CALL_SUBTEST_3(( jacobisvd_verify_assert() )); + CALL_SUBTEST_3(( jacobisvd_verify_assert() )); + CALL_SUBTEST_9(( jacobisvd_verify_assert() )); + CALL_SUBTEST_11(( jacobisvd_verify_assert() )); // Test problem size constructors CALL_SUBTEST_12( JacobiSVD(10, 20) );