From 7aa6fd362558718304937ddfd60232c9802d69be Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Thu, 3 Sep 2009 02:53:51 -0400 Subject: [PATCH] big reorganization in JacobiSVD: - R-SVD preconditioning now done with meta selectors to avoid compiling useless code - SVD options now honored, with options to hint "at least as many rows as cols" etc... - fix compilation in bad cases (rectangular and fixed-size) - the check for termination is now done on the fly, no more goto (should have done that earlier!) --- Eigen/src/Core/util/Constants.h | 9 ++ Eigen/src/QR/FullPivotingHouseholderQR.h | 6 +- Eigen/src/SVD/JacobiSVD.h | 184 ++++++++++++++++------- test/jacobisvd.cpp | 37 +++-- 4 files changed, 156 insertions(+), 80 deletions(-) diff --git a/Eigen/src/Core/util/Constants.h b/Eigen/src/Core/util/Constants.h index 21e7f62c3..affc1d478 100644 --- a/Eigen/src/Core/util/Constants.h +++ b/Eigen/src/Core/util/Constants.h @@ -238,6 +238,15 @@ 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. diff --git a/Eigen/src/QR/FullPivotingHouseholderQR.h b/Eigen/src/QR/FullPivotingHouseholderQR.h index 52405170b..0d542cf7a 100644 --- a/Eigen/src/QR/FullPivotingHouseholderQR.h +++ b/Eigen/src/QR/FullPivotingHouseholderQR.h @@ -67,12 +67,10 @@ template class FullPivotingHouseholderQR * The default constructor is useful in cases in which the user intends to * perform decompositions via FullPivotingHouseholderQR::compute(const MatrixType&). */ - FullPivotingHouseholderQR() : m_qr(), m_hCoeffs(), m_isInitialized(false) {} + FullPivotingHouseholderQR() : m_isInitialized(false) {} FullPivotingHouseholderQR(const MatrixType& matrix) - : m_qr(matrix.rows(), matrix.cols()), - m_hCoeffs(std::min(matrix.rows(),matrix.cols())), - m_isInitialized(false) + : m_isInitialized(false) { compute(matrix); } diff --git a/Eigen/src/SVD/JacobiSVD.h b/Eigen/src/SVD/JacobiSVD.h index ca359832b..2801ee077 100644 --- a/Eigen/src/SVD/JacobiSVD.h +++ b/Eigen/src/SVD/JacobiSVD.h @@ -33,8 +33,11 @@ * \brief Jacobi SVD decomposition of a square matrix * * \param MatrixType the type of the matrix of which we are computing the SVD decomposition - * \param ComputeU whether the U matrix should be computed - * \param ComputeV whether the V matrix should be computed + * \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. * * \sa MatrixBase::jacobiSvd() */ @@ -44,14 +47,14 @@ template class JacobiSVD typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits::Real RealScalar; enum { - ComputeU = 1, - ComputeV = 1, + ComputeU = (Options & SkipU) == 0, + ComputeV = (Options & SkipV) == 0, RowsAtCompileTime = MatrixType::RowsAtCompileTime, ColsAtCompileTime = MatrixType::ColsAtCompileTime, - DiagSizeAtCompileTime = EIGEN_ENUM_MIN(RowsAtCompileTime,ColsAtCompileTime), + DiagSizeAtCompileTime = EIGEN_SIZE_MIN(RowsAtCompileTime,ColsAtCompileTime), MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, - MaxDiagSizeAtCompileTime = EIGEN_ENUM_MIN(MaxRowsAtCompileTime,MaxColsAtCompileTime), + MaxDiagSizeAtCompileTime = EIGEN_SIZE_MIN(MaxRowsAtCompileTime,MaxColsAtCompileTime), MatrixOptions = MatrixType::Options }; @@ -65,9 +68,12 @@ template class JacobiSVD MatrixOptions, MaxColsAtCompileTime, MaxColsAtCompileTime>, DummyMatrixType>::ret MatrixVType; typedef Matrix SingularValuesType; + MatrixOptions, MaxDiagSizeAtCompileTime, 1> SingularValuesType; typedef Matrix RowType; typedef Matrix ColType; + typedef Matrix + WorkMatrixType; public: @@ -92,7 +98,7 @@ template class JacobiSVD return m_singularValues; } - const MatrixUType& matrixV() const + const MatrixVType& matrixV() const { ei_assert(m_isInitialized && "JacobiSVD is not initialized."); return m_matrixV; @@ -106,12 +112,17 @@ template class JacobiSVD 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::IsComplex> struct ei_svd_precondition_2x2_block_to_be_real { - static void run(MatrixType&, JacobiSVD&, int, int) {} + typedef JacobiSVD SVD; + static void run(typename SVD::WorkMatrixType&, JacobiSVD&, int, int) {} }; template @@ -120,9 +131,8 @@ struct ei_svd_precondition_2x2_block_to_be_real typedef JacobiSVD SVD; typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; - enum { ComputeU = SVD::ComputeU, ComputeV = SVD::ComputeV }; - static void run(MatrixType& work_matrix, JacobiSVD& svd, int p, int q) + static void run(typename SVD::WorkMatrixType& work_matrix, JacobiSVD& svd, int p, int q) { Scalar z; PlanarRotation rot; @@ -185,10 +195,94 @@ void ei_real_2x2_jacobi_svd(const MatrixType& matrix, int p, int q, *j_left = rot1 * j_right->transpose(); } +templateMatrixType::ColsAtCompileTime)> +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; + enum { ComputeU = SVD::ComputeU, ComputeV = SVD::ComputeV }; + static bool run(const MatrixType& matrix, typename SVD::WorkMatrixType& work_matrix, SVD& svd) + { + int rows = matrix.rows(); + int cols = matrix.cols(); + int diagSize = cols; + if(rows > cols) + { + FullPivotingHouseholderQR qr(matrix); + work_matrix = qr.matrixQR().block(0,0,diagSize,diagSize).template triangularView(); + if(ComputeU) svd.m_matrixU = qr.matrixQ(); + if(ComputeV) + for(int i = 0; i < cols; i++) + svd.m_matrixV.coeffRef(qr.colsPermutation().coeff(i),i) = Scalar(1); + return true; + } + else return false; + } +}; + +templateMatrixType::RowsAtCompileTime)> +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; + 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) + { + int rows = matrix.rows(); + int cols = matrix.cols(); + int diagSize = rows; + if(cols > rows) + { + typedef Matrix + TransposeTypeWithSameStorageOrder; + FullPivotingHouseholderQR qr(matrix.adjoint()); + work_matrix = qr.matrixQR().block(0,0,diagSize,diagSize).template triangularView().adjoint(); + if(ComputeV) svd.m_matrixV = qr.matrixQ(); + if(ComputeU) + for(int i = 0; i < rows; i++) + svd.m_matrixU.coeffRef(qr.colsPermutation().coeff(i),i) = Scalar(1); + return true; + } + else return false; + } +}; + template JacobiSVD& JacobiSVD::compute(const MatrixType& matrix) { - MatrixType work_matrix; + WorkMatrixType work_matrix; int rows = matrix.rows(); int cols = matrix.cols(); int diagSize = std::min(rows, cols); @@ -197,65 +291,41 @@ JacobiSVD& JacobiSVD::compute(const Ma m_singularValues.resize(diagSize); const RealScalar precision = 2 * epsilon(); - if(rows > cols) + if(!ei_svd_precondition_if_more_rows_than_cols::run(matrix, work_matrix, *this) + && !ei_svd_precondition_if_more_cols_than_rows::run(matrix, work_matrix, *this)) { - FullPivotingHouseholderQR qr(matrix); - work_matrix = qr.matrixQR().block(0,0,diagSize,diagSize).template triangularView(); - if(ComputeU) m_matrixU = qr.matrixQ(); - if(ComputeV) - for(int i = 0; i < cols; i++) - m_matrixV.coeffRef(qr.colsPermutation().coeff(i),i) = Scalar(1); - } - else if(rows < cols) - { - FullPivotingHouseholderQR qr(MatrixType(matrix.adjoint())); - work_matrix = qr.matrixQR().block(0,0,diagSize,diagSize).template triangularView().adjoint(); - if(ComputeV) m_matrixV = qr.matrixQ(); - if(ComputeU) - for(int i = 0; i < rows; i++) - m_matrixU.coeffRef(qr.colsPermutation().coeff(i),i) = Scalar(1); - - } - else - { - work_matrix = matrix; + work_matrix = matrix.block(0,0,diagSize,diagSize); if(ComputeU) m_matrixU.diagonal().setOnes(); if(ComputeV) m_matrixV.diagonal().setOnes(); } -sweep_again: - for(int p = 1; p < diagSize; ++p) + + bool finished = false; + while(!finished) { - for(int q = 0; q < p; ++q) + finished = true; + for(int p = 1; p < diagSize; ++p) { - if(std::max(ei_abs(work_matrix.coeff(p,q)),ei_abs(work_matrix.coeff(q,p))) - > std::max(ei_abs(work_matrix.coeff(p,p)),ei_abs(work_matrix.coeff(q,q)))*precision) + for(int q = 0; q < p; ++q) { - ei_svd_precondition_2x2_block_to_be_real::run(work_matrix, *this, p, q); + if(std::max(ei_abs(work_matrix.coeff(p,q)),ei_abs(work_matrix.coeff(q,p))) + > std::max(ei_abs(work_matrix.coeff(p,p)),ei_abs(work_matrix.coeff(q,q)))*precision) + { + finished = false; + ei_svd_precondition_2x2_block_to_be_real::run(work_matrix, *this, p, q); - PlanarRotation j_left, j_right; - ei_real_2x2_jacobi_svd(work_matrix, p, q, &j_left, &j_right); + PlanarRotation j_left, j_right; + ei_real_2x2_jacobi_svd(work_matrix, p, q, &j_left, &j_right); - work_matrix.applyOnTheLeft(p,q,j_left); - if(ComputeU) m_matrixU.applyOnTheRight(p,q,j_left.transpose()); + work_matrix.applyOnTheLeft(p,q,j_left); + if(ComputeU) m_matrixU.applyOnTheRight(p,q,j_left.transpose()); - work_matrix.applyOnTheRight(p,q,j_right); - if(ComputeV) m_matrixV.applyOnTheRight(p,q,j_right); + work_matrix.applyOnTheRight(p,q,j_right); + if(ComputeV) m_matrixV.applyOnTheRight(p,q,j_right); + } } } } - RealScalar biggestOnDiag = work_matrix.diagonal().cwise().abs().maxCoeff(); - RealScalar maxAllowedOffDiag = biggestOnDiag * precision; - for(int p = 0; p < diagSize; ++p) - { - for(int q = 0; q < p; ++q) - if(ei_abs(work_matrix.coeff(p,q)) > maxAllowedOffDiag) - goto sweep_again; - for(int q = p+1; q < diagSize; ++q) - if(ei_abs(work_matrix.coeff(p,q)) > maxAllowedOffDiag) - goto sweep_again; - } - for(int i = 0; i < diagSize; ++i) { RealScalar a = ei_abs(work_matrix.coeff(i,i)); diff --git a/test/jacobisvd.cpp b/test/jacobisvd.cpp index 8b4c7584e..5940b8497 100644 --- a/test/jacobisvd.cpp +++ b/test/jacobisvd.cpp @@ -27,7 +27,7 @@ #include #include -template void svd(const MatrixType& m, bool pickrandom = true) +template void svd(const MatrixType& m = MatrixType(), bool pickrandom = true) { int rows = m.rows(); int cols = m.cols(); @@ -48,7 +48,7 @@ template void svd(const MatrixType& m, bool pickrandom = tr if(pickrandom) a = MatrixType::Random(rows,cols); else a = m; - JacobiSVD svd(a); + JacobiSVD svd(a); MatrixType sigma = MatrixType::Zero(rows,cols); sigma.diagonal() = svd.singularValues().template cast(); MatrixUType u = svd.matrixU(); @@ -80,28 +80,27 @@ void test_jacobisvd() Matrix2cd m; m << 0, 1, 0, 1; - CALL_SUBTEST( svd(m, false) ); + CALL_SUBTEST(( svd(m, false) )); m << 1, 0, 1, 0; - CALL_SUBTEST( svd(m, false) ); + CALL_SUBTEST(( svd(m, false) )); Matrix2d n; n << 1, 1, 1, -1; - CALL_SUBTEST( svd(n, false) ); - CALL_SUBTEST( svd(Matrix3f()) ); - CALL_SUBTEST( svd(Matrix4d()) ); - CALL_SUBTEST( svd(MatrixXf(50,50)) ); - CALL_SUBTEST( svd(MatrixXcd(14,7)) ); - CALL_SUBTEST( svd(MatrixXd(10,50)) ); - - CALL_SUBTEST( svd(MatrixXcf(3,3)) ); - CALL_SUBTEST( svd(MatrixXd(30,30)) ); + CALL_SUBTEST(( svd(n, false) )); + CALL_SUBTEST(( svd() )); + CALL_SUBTEST(( svd() )); + CALL_SUBTEST(( svd , AtLeastAsManyColsAsRows>() )); + CALL_SUBTEST(( svd , AtLeastAsManyRowsAsCols>(Matrix(10,2)) )); + + CALL_SUBTEST(( svd(MatrixXf(50,50)) )); + CALL_SUBTEST(( svd(MatrixXcd(14,7)) )); } - CALL_SUBTEST( svd(MatrixXf(300,200)) ); - CALL_SUBTEST( svd(MatrixXcd(100,150)) ); + CALL_SUBTEST(( svd(MatrixXf(300,200)) )); + CALL_SUBTEST(( svd(MatrixXcd(100,150)) )); - CALL_SUBTEST( svd_verify_assert() ); - CALL_SUBTEST( svd_verify_assert() ); - CALL_SUBTEST( svd_verify_assert() ); - CALL_SUBTEST( svd_verify_assert() ); + CALL_SUBTEST(( svd_verify_assert() )); + CALL_SUBTEST(( svd_verify_assert() )); + CALL_SUBTEST(( svd_verify_assert() )); + CALL_SUBTEST(( svd_verify_assert() )); }