From 8f0e80fe304ad34b520a869ae58bbff8b64c63f2 Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Thu, 14 Oct 2010 10:14:43 -0400 Subject: [PATCH] JacobiSVD: * fix preallocating constructors, allocate U and V of the right size for computation options * complete documentation and internal comments * improve unit test, test inf/nan values --- Eigen/src/SVD/JacobiSVD.h | 202 ++++++++++++++++++++++--------- doc/snippets/JacobiSVD_basic.cpp | 9 ++ test/jacobisvd.cpp | 36 +++++- 3 files changed, 186 insertions(+), 61 deletions(-) create mode 100644 doc/snippets/JacobiSVD_basic.cpp diff --git a/Eigen/src/SVD/JacobiSVD.h b/Eigen/src/SVD/JacobiSVD.h index e7156a3fb..e50c9ed15 100644 --- a/Eigen/src/SVD/JacobiSVD.h +++ b/Eigen/src/SVD/JacobiSVD.h @@ -25,14 +25,19 @@ #ifndef EIGEN_JACOBISVD_H #define EIGEN_JACOBISVD_H -// forward declarations (needed by ICC) -// the empty bodies are required by MSVC +// forward declaration (needed by ICC) +// the empty body is required by MSVC template::IsComplex> struct ei_svd_precondition_2x2_block_to_be_real {}; -/*** QR preconditioners (R-SVD) ***/ +/*** QR preconditioners (R-SVD) + *** + *** Their role is to reduce the problem of computing the SVD to the case of a square matrix. + *** This approach, known as R-SVD, is an optimization for rectangular-enough matrices, and is a requirement for + *** JacobiSVD which by itself is only able to work on square matrices. + ***/ enum { PreconditionIfMoreColsThanRows, PreconditionIfMoreRowsThanCols }; @@ -40,11 +45,11 @@ template struct ei_qr_preconditioner_should_do_anything { enum { a = MatrixType::RowsAtCompileTime != Dynamic && - MatrixType::ColsAtCompileTime != Dynamic && - MatrixType::ColsAtCompileTime <= MatrixType::RowsAtCompileTime, + MatrixType::ColsAtCompileTime != Dynamic && + MatrixType::ColsAtCompileTime <= MatrixType::RowsAtCompileTime, b = MatrixType::RowsAtCompileTime != Dynamic && - MatrixType::ColsAtCompileTime != Dynamic && - MatrixType::RowsAtCompileTime <= MatrixType::ColsAtCompileTime, + MatrixType::ColsAtCompileTime != Dynamic && + MatrixType::RowsAtCompileTime <= MatrixType::ColsAtCompileTime, ret = !( (QRPreconditioner == NoQRPreconditioner) || (Case == PreconditionIfMoreColsThanRows && bool(a)) || (Case == PreconditionIfMoreRowsThanCols && bool(b)) ) @@ -64,6 +69,8 @@ struct ei_qr_preconditioner_impl } }; +/*** preconditioner using FullPivHouseholderQR ***/ + template struct ei_qr_preconditioner_impl { @@ -101,6 +108,8 @@ struct ei_qr_preconditioner_impl struct ei_qr_preconditioner_impl { @@ -146,6 +155,8 @@ struct ei_qr_preconditioner_impl struct ei_qr_preconditioner_impl { @@ -198,12 +209,36 @@ struct ei_qr_preconditioner_impl class JacobiSVD /** \brief Default Constructor with memory preallocation * * Like the default constructor but with preallocation of the internal data - * according to the specified problem \a size. + * according to the specified problem size. * \sa JacobiSVD() */ - JacobiSVD(Index rows, Index cols) : m_matrixU(rows, rows), - m_matrixV(cols, cols), - m_singularValues(std::min(rows, cols)), - m_workMatrix(rows, cols), - m_isInitialized(false) {} + JacobiSVD(Index rows, Index cols, unsigned int computationOptions = 0) + { + allocate(rows, cols, computationOptions); + } /** \brief Constructor performing the decomposition of given matrix. * @@ -286,15 +320,7 @@ template class JacobiSVD * 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()), - 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, computationOptions); } @@ -315,6 +341,15 @@ template class JacobiSVD */ JacobiSVD& compute(const MatrixType& matrix, unsigned int computationOptions = 0); + /** \returns the \a U matrix. + * + * For the SVD decomposition of a n-by-p matrix, letting \a m be the minimum of \a n and \a p, + * the U matrix is n-by-n if you asked for ComputeFullU, and is n-by-m if you asked for ComputeThinU. + * + * The \a m first columns of \a U are the left singular vectors of the matrix being decomposed. + * + * This method asserts that you asked for \a U to be computed. + */ const MatrixUType& matrixU() const { ei_assert(m_isInitialized && "JacobiSVD is not initialized."); @@ -322,12 +357,15 @@ template class JacobiSVD return m_matrixU; } - const SingularValuesType& singularValues() const - { - ei_assert(m_isInitialized && "JacobiSVD is not initialized."); - return m_singularValues; - } - + /** \returns the \a V matrix. + * + * For the SVD decomposition of a n-by-p matrix, letting \a m be the minimum of \a n and \a p, + * the V matrix is p-by-p if you asked for ComputeFullV, and is p-by-m if you asked for ComputeThinV. + * + * The \a m first columns of \a V are the right singular vectors of the matrix being decomposed. + * + * This method asserts that you asked for \a V to be computed. + */ const MatrixVType& matrixV() const { ei_assert(m_isInitialized && "JacobiSVD is not initialized."); @@ -335,14 +373,27 @@ template class JacobiSVD return m_matrixV; } + /** \returns the vector of singular values. + * + * For the SVD decomposition of a n-by-p matrix, letting \a m be the minimum of \a n and \a p, the + * returned vector has size \a m. + */ + const SingularValuesType& singularValues() const + { + ei_assert(m_isInitialized && "JacobiSVD is not initialized."); + return m_singularValues; + } + + /** \returns true if \a U (full or thin) is asked for in this SVD decomposition */ inline bool computeU() const { return m_computeFullU || m_computeThinU; } + /** \returns true if \a V (full or thin) is asked for in this SVD decomposition */ inline bool computeV() const { return m_computeFullV || m_computeThinV; } /** \returns a (least squares) solution of \f$ A x = b \f$ using the current SVD decomposition of A. * * \param b the right-hand-side of the equation to solve. * - * \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, + * \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. * * \note SVD solving is implicitly least-squares. Thus, this method serves both purposes of exact solving and least-squares solving. * In other words, the returned solution is guaranteed to minimize the Euclidean norm \f$ \Vert A x - b \Vert \f$. @@ -356,6 +407,7 @@ template class JacobiSVD return ei_solve_retval(*this, b.derived()); } + /** \returns the number of singular values that are not exactly 0 */ Index nonzeroSingularValues() const { ei_assert(m_isInitialized && "JacobiSVD is not initialized."); @@ -365,6 +417,9 @@ template class JacobiSVD inline Index rows() const { return m_rows; } inline Index cols() const { return m_cols; } + private: + void allocate(Index rows, Index cols, unsigned int computationOptions = 0); + protected: MatrixUType m_matrixU; MatrixVType m_matrixV; @@ -373,7 +428,7 @@ template class JacobiSVD bool m_isInitialized; bool m_computeFullU, m_computeThinU; bool m_computeFullV, m_computeThinV; - Index m_nonzeroSingularValues, m_rows, m_cols; + Index m_nonzeroSingularValues, m_rows, m_cols, m_diagSize; template friend struct ei_svd_precondition_2x2_block_to_be_real; @@ -381,6 +436,38 @@ template class JacobiSVD friend struct ei_qr_preconditioner_impl; }; +template +void JacobiSVD::allocate(Index rows, Index cols, unsigned int computationOptions) +{ + m_rows = rows; + m_cols = cols; + m_isInitialized = false; + 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."); + if (QRPreconditioner == FullPivHouseholderQRPreconditioner) + { + ei_assert(!(m_computeThinU || m_computeThinV) && + "JacobiSVD: can't compute thin U or thin V with the FullPivHouseholderQR preconditioner. " + "Use the ColPivHouseholderQR preconditioner instead."); + } + m_diagSize = std::min(m_rows, m_cols); + m_singularValues.resize(m_diagSize); + m_matrixU.resize(m_rows, m_computeFullU ? m_rows + : m_computeThinU ? m_diagSize + : 0); + m_matrixV.resize(m_cols, m_computeFullV ? m_cols + : m_computeThinV ? m_diagSize + : 0); + m_workMatrix.resize(m_diagSize, m_diagSize); +} + + template struct ei_svd_precondition_2x2_block_to_be_real { @@ -463,53 +550,51 @@ template JacobiSVD& JacobiSVD::compute(const MatrixType& matrix, unsigned int computationOptions) { - 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."); - if (QRPreconditioner == FullPivHouseholderQRPreconditioner) - { - ei_assert(!(m_computeThinU || m_computeThinV) && - "JacobiSVD: can't compute thin U or thin V with the FullPivHouseholderQR preconditioner. " - "Use the ColPivHouseholderQR preconditioner instead."); - } - m_rows = matrix.rows(); - m_cols = matrix.cols(); - Index diagSize = std::min(m_rows, m_cols); - m_singularValues.resize(diagSize); + allocate(matrix.rows(), matrix.cols(), computationOptions); + + // currently we stop when we reach precision 2*epsilon as the last bit of precision can require an unreasonable number of iterations, + // only worsening the precision of U and V as we accumulate more rotations const RealScalar precision = 2 * NumTraits::epsilon(); + /*** step 1. The R-SVD step: we use a QR decomposition to reduce to the case of a square matrix */ + if(!ei_qr_preconditioner_impl::run(*this, matrix) && !ei_qr_preconditioner_impl::run(*this, matrix)) { - m_workMatrix = matrix.block(0,0,diagSize,diagSize); + m_workMatrix = matrix.block(0,0,m_diagSize,m_diagSize); if(m_computeFullU) m_matrixU.setIdentity(m_rows,m_rows); - if(m_computeThinU) m_matrixU.setIdentity(m_rows,diagSize); + if(m_computeThinU) m_matrixU.setIdentity(m_rows,m_diagSize); if(m_computeFullV) m_matrixV.setIdentity(m_cols,m_cols); - if(m_computeThinV) m_matrixV.setIdentity(diagSize,m_cols); + if(m_computeThinV) m_matrixV.setIdentity(m_cols, m_diagSize); } + /*** step 2. The main Jacobi SVD iteration. ***/ + bool finished = false; while(!finished) { finished = true; - for(Index p = 1; p < diagSize; ++p) + + // do a sweep: for all index pairs (p,q), perform SVD of the corresponding 2x2 sub-matrix + + for(Index p = 1; p < m_diagSize; ++p) { for(Index q = 0; q < p; ++q) { + // if this 2x2 sub-matrix is not diagonal already... + // notice that this comparison will evaluate to false if any NaN is involved, ensuring that NaN's don't + // keep us iterating forever. if(std::max(ei_abs(m_workMatrix.coeff(p,q)),ei_abs(m_workMatrix.coeff(q,p))) > 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); + // perform SVD decomposition of 2x2 sub-matrix corresponding to indices p,q to make it diagonal + 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); + // accumulate resulting Jacobi rotations m_workMatrix.applyOnTheLeft(p,q,j_left); if(computeU()) m_matrixU.applyOnTheRight(p,q,j_left.transpose()); @@ -520,19 +605,22 @@ JacobiSVD::compute(const MatrixType& matrix, unsig } } - for(Index i = 0; i < diagSize; ++i) + /*** step 3. The work matrix is now diagonal, so ensure it's positive so its diagonal entries are the singular values ***/ + + for(Index i = 0; i < m_diagSize; ++i) { 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; } - m_nonzeroSingularValues = diagSize; + /*** step 4. Sort singular values in descending order and compute the number of nonzero singular values ***/ - for(Index i = 0; i < diagSize; i++) + m_nonzeroSingularValues = m_diagSize; + for(Index i = 0; i < m_diagSize; i++) { Index pos; - RealScalar maxRemainingSingularValue = m_singularValues.tail(diagSize-i).maxCoeff(&pos); + RealScalar maxRemainingSingularValue = m_singularValues.tail(m_diagSize-i).maxCoeff(&pos); if(maxRemainingSingularValue == RealScalar(0)) { m_nonzeroSingularValues = i; diff --git a/doc/snippets/JacobiSVD_basic.cpp b/doc/snippets/JacobiSVD_basic.cpp new file mode 100644 index 000000000..ab24b9bca --- /dev/null +++ b/doc/snippets/JacobiSVD_basic.cpp @@ -0,0 +1,9 @@ +MatrixXf m = MatrixXf::Random(3,2); +cout << "Here is the matrix m:" << endl << m << endl; +JacobiSVD svd(m, ComputeThinU | ComputeThinV); +cout << "Its singular values are:" << endl << svd.singularValues() << endl; +cout << "Its left singular vectors are the columns of the thin U matrix:" << endl << svd.matrixU() << endl; +cout << "Its right singular vectors are the columns of the thin V matrix:" << endl << svd.matrixV() << endl; +Vector3f rhs(1, 0, 0); +cout << "Now consider this rhs vector:" << endl << rhs << endl; +cout << "A least-squares solution of m*x = rhs is:" << endl << svd.solve(rhs) << endl; diff --git a/test/jacobisvd.cpp b/test/jacobisvd.cpp index a6dbcf2e8..7d8c36308 100644 --- a/test/jacobisvd.cpp +++ b/test/jacobisvd.cpp @@ -196,6 +196,23 @@ template void jacobisvd_verify_assert(const MatrixType& m) } } +template +void jacobisvd_inf_nan() +{ + JacobiSVD svd; + typedef typename MatrixType::Scalar Scalar; + Scalar some_inf = Scalar(1) / Scalar(0); + svd.compute(MatrixType::Constant(10,10,some_inf), ComputeFullU | ComputeFullV); + Scalar some_nan = Scalar(0) / Scalar(0); + svd.compute(MatrixType::Constant(10,10,some_nan), ComputeFullU | ComputeFullV); + MatrixType m = MatrixType::Zero(10,10); + m(ei_random(0,9), ei_random(0,9)) = some_inf; + svd.compute(m, ComputeFullU | ComputeFullV); + m = MatrixType::Zero(10,10); + m(ei_random(0,9), ei_random(0,9)) = some_nan; + svd.compute(m, ComputeFullU | ComputeFullV); +} + void test_jacobisvd() { CALL_SUBTEST_3(( jacobisvd_verify_assert(Matrix3f()) )); @@ -211,10 +228,15 @@ void test_jacobisvd() m << 1, 0, 1, 0; CALL_SUBTEST_1(( jacobisvd(m, false) )); + Matrix2d n; - n << 1, 1, - 1, -1; + n << 0, 0, + 0, 0; CALL_SUBTEST_2(( jacobisvd(n, false) )); + n << 0, 0, + 0, 1; + CALL_SUBTEST_2(( jacobisvd(n, false) )); + CALL_SUBTEST_3(( jacobisvd() )); CALL_SUBTEST_4(( jacobisvd() )); CALL_SUBTEST_5(( jacobisvd >() )); @@ -226,8 +248,14 @@ void test_jacobisvd() CALL_SUBTEST_8(( jacobisvd(MatrixXcd(r,c)) )); (void) r; (void) c; + + // Test on inf/nan matrix + CALL_SUBTEST_7( jacobisvd_inf_nan() ); } - CALL_SUBTEST_7(( jacobisvd(MatrixXf(ei_random(200, 300), ei_random(200, 300))) )); - CALL_SUBTEST_8(( jacobisvd(MatrixXcd(ei_random(100, 120), ei_random(100, 120))) )); + CALL_SUBTEST_7(( jacobisvd(MatrixXf(ei_random(100, 150), ei_random(100, 150))) )); + CALL_SUBTEST_8(( jacobisvd(MatrixXcd(ei_random(80, 100), ei_random(80, 100))) )); + + // Test problem size constructors + CALL_SUBTEST_7( JacobiSVD(10,10) ); }