From 5bbc9cea93ef29cee2b8ffb2084d4ebca32600ba Mon Sep 17 00:00:00 2001 From: Rasmus Munk Larsen Date: Wed, 31 Mar 2021 21:09:19 +0000 Subject: Add an info() method to the SVDBase class to make it possible to tell the user that the computation failed, possibly due to invalid input. Make Jacobi and divide-and-conquer fail fast and return info() == InvalidInput if the matrix contains NaN or +/-Inf. --- Eigen/src/SVD/BDCSVD.h | 32 +++++++++++++++++++++++++------- Eigen/src/SVD/JacobiSVD.h | 9 ++++++++- Eigen/src/SVD/SVDBase.h | 41 +++++++++++++++++++++++++++++++---------- 3 files changed, 64 insertions(+), 18 deletions(-) (limited to 'Eigen') diff --git a/Eigen/src/SVD/BDCSVD.h b/Eigen/src/SVD/BDCSVD.h index e0c4456c7..17f8e4436 100644 --- a/Eigen/src/SVD/BDCSVD.h +++ b/Eigen/src/SVD/BDCSVD.h @@ -208,6 +208,7 @@ protected: using Base::m_computeThinV; using Base::m_matrixU; using Base::m_matrixV; + using Base::m_info; using Base::m_isInitialized; using Base::m_nonzeroSingularValues; @@ -256,16 +257,25 @@ BDCSVD& BDCSVD::compute(const MatrixType& matrix, unsign { // FIXME this line involves temporaries JacobiSVD jsvd(matrix,computationOptions); - if(computeU()) m_matrixU = jsvd.matrixU(); - if(computeV()) m_matrixV = jsvd.matrixV(); - m_singularValues = jsvd.singularValues(); - m_nonzeroSingularValues = jsvd.nonzeroSingularValues(); m_isInitialized = true; + m_info = jsvd.info(); + if (m_info == Success || m_info == NoConvergence) { + if(computeU()) m_matrixU = jsvd.matrixU(); + if(computeV()) m_matrixV = jsvd.matrixV(); + m_singularValues = jsvd.singularValues(); + m_nonzeroSingularValues = jsvd.nonzeroSingularValues(); + } return *this; } //**** step 0 - Copy the input matrix and apply scaling to reduce over/under-flows - RealScalar scale = matrix.cwiseAbs().maxCoeff(); + RealScalar scale = matrix.cwiseAbs().template maxCoeff(); + if (!(numext::isfinite)(scale)) { + m_isInitialized = true; + m_info = InvalidInput; + return *this; + } + if(scale==Literal(0)) scale = Literal(1); MatrixX copy; if (m_isTranspose) copy = matrix.adjoint()/scale; @@ -282,7 +292,11 @@ BDCSVD& BDCSVD::compute(const MatrixType& matrix, unsign m_computed.topRows(m_diagSize) = bid.bidiagonal().toDenseMatrix().transpose(); m_computed.template bottomRows<1>().setZero(); divide(0, m_diagSize - 1, 0, 0, 0); - + if (m_info != Success && m_info != NoConvergence) { + m_isInitialized = true; + return *this; + } + //**** step 3 - Copy singular values and vectors for (int i=0; i::structured_update(Block A, co //@param shift : Each time one takes the left submatrix, one must add 1 to the shift. Why? Because! We actually want the last column of the U submatrix // to become the first column (*coeff) and to shift all the other columns to the right. There are more details on the reference paper. template -void BDCSVD::divide (Eigen::Index firstCol, Eigen::Index lastCol, Eigen::Index firstRowW, Eigen::Index firstColW, Eigen::Index shift) +void BDCSVD::divide(Eigen::Index firstCol, Eigen::Index lastCol, Eigen::Index firstRowW, Eigen::Index firstColW, Eigen::Index shift) { // requires rows = cols + 1; using std::pow; @@ -414,6 +428,8 @@ void BDCSVD::divide (Eigen::Index firstCol, Eigen::Index lastCol, Ei { // FIXME this line involves temporaries JacobiSVD b(m_computed.block(firstCol, firstCol, n + 1, n), ComputeFullU | (m_compV ? ComputeFullV : 0)); + m_info = b.info(); + if (m_info != Success && m_info != NoConvergence) return; if (m_compU) m_naiveU.block(firstCol, firstCol, n + 1, n + 1).real() = b.matrixU(); else @@ -433,7 +449,9 @@ void BDCSVD::divide (Eigen::Index firstCol, Eigen::Index lastCol, Ei // and the divide of the right submatrice reads one column of the left submatrice. That's why we need to treat the // right submatrix before the left one. divide(k + 1 + firstCol, lastCol, k + 1 + firstRowW, k + 1 + firstColW, shift); + if (m_info != Success && m_info != NoConvergence) return; divide(firstCol, k - 1 + firstCol, firstRowW, firstColW + 1, shift + 1); + if (m_info != Success && m_info != NoConvergence) return; if (m_compU) { diff --git a/Eigen/src/SVD/JacobiSVD.h b/Eigen/src/SVD/JacobiSVD.h index 2b6891105..a22a2e5c3 100644 --- a/Eigen/src/SVD/JacobiSVD.h +++ b/Eigen/src/SVD/JacobiSVD.h @@ -585,6 +585,7 @@ template class JacobiSVD using Base::m_matrixU; using Base::m_matrixV; using Base::m_singularValues; + using Base::m_info; using Base::m_isInitialized; using Base::m_isAllocated; using Base::m_usePrescribedThreshold; @@ -625,6 +626,7 @@ void JacobiSVD::allocate(Eigen::Index rows, Eigen: m_rows = rows; m_cols = cols; + m_info = Success; m_isInitialized = false; m_isAllocated = true; m_computationOptions = computationOptions; @@ -674,7 +676,12 @@ JacobiSVD::compute(const MatrixType& matrix, unsig const RealScalar considerAsZero = (std::numeric_limits::min)(); // Scaling factor to reduce over/under-flows - RealScalar scale = matrix.cwiseAbs().maxCoeff(); + RealScalar scale = matrix.cwiseAbs().template maxCoeff(); + if (!(numext::isfinite)(scale)) { + m_isInitialized = true; + m_info = InvalidInput; + return *this; + } if(scale==RealScalar(0)) scale = RealScalar(1); /*** step 1. The R-SVD step: we use a QR decomposition to reduce to the case of a square matrix */ diff --git a/Eigen/src/SVD/SVDBase.h b/Eigen/src/SVD/SVDBase.h index 34d5c9dd3..13c8394ae 100644 --- a/Eigen/src/SVD/SVDBase.h +++ b/Eigen/src/SVD/SVDBase.h @@ -52,7 +52,7 @@ template struct traits > * 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, * 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. * - * If the input matrix has inf or nan coefficients, the result of the computation is undefined, but the computation is guaranteed to + * If the input matrix has inf or nan coefficients, the result of the computation is undefined, and \a info() will return \a InvalidInput, but the computation is guaranteed to * terminate in finite (and reasonable) time. * \sa class BDCSVD, class JacobiSVD */ @@ -97,7 +97,7 @@ public: */ const MatrixUType& matrixU() const { - eigen_assert(m_isInitialized && "SVD is not initialized."); + _check_compute_assertions(); eigen_assert(computeU() && "This SVD decomposition didn't compute U. Did you ask for it?"); return m_matrixU; } @@ -113,7 +113,7 @@ public: */ const MatrixVType& matrixV() const { - eigen_assert(m_isInitialized && "SVD is not initialized."); + _check_compute_assertions(); eigen_assert(computeV() && "This SVD decomposition didn't compute V. Did you ask for it?"); return m_matrixV; } @@ -125,14 +125,14 @@ public: */ const SingularValuesType& singularValues() const { - eigen_assert(m_isInitialized && "SVD is not initialized."); + _check_compute_assertions(); return m_singularValues; } /** \returns the number of singular values that are not exactly 0 */ Index nonzeroSingularValues() const { - eigen_assert(m_isInitialized && "SVD is not initialized."); + _check_compute_assertions(); return m_nonzeroSingularValues; } @@ -145,7 +145,7 @@ public: inline Index rank() const { using std::abs; - eigen_assert(m_isInitialized && "JacobiSVD is not initialized."); + _check_compute_assertions(); if(m_singularValues.size()==0) return 0; RealScalar premultiplied_threshold = numext::maxi(m_singularValues.coeff(0) * threshold(), (std::numeric_limits::min)()); Index i = m_nonzeroSingularValues-1; @@ -224,6 +224,18 @@ public: solve(const MatrixBase& b) const; #endif + + /** \brief Reports whether previous computation was successful. + * + * \returns \c Success if computation was successful. + */ + EIGEN_DEVICE_FUNC + ComputationInfo info() const + { + eigen_assert(m_isInitialized && "SVD is not initialized."); + return m_info; + } + #ifndef EIGEN_PARSED_BY_DOXYGEN template void _solve_impl(const RhsType &rhs, DstType &dst) const; @@ -233,26 +245,33 @@ public: #endif protected: - + static void check_template_parameters() { EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar); } + void _check_compute_assertions() const { + eigen_assert(m_isInitialized && "SVD is not initialized."); + eigen_assert(m_info != InvalidInput && "SVD failed due to invalid input."); + eigen_assert(m_info != NumericalIssue && "SVD failed due to invalid input."); + } + template void _check_solve_assertion(const Rhs& b) const { EIGEN_ONLY_USED_FOR_DEBUG(b); - eigen_assert(m_isInitialized && "SVD is not initialized."); + _check_compute_assertions(); eigen_assert(computeU() && computeV() && "SVDBase::solve(): Both unitaries U and V are required to be computed (thin unitaries suffice)."); eigen_assert((Transpose_?cols():rows())==b.rows() && "SVDBase::solve(): invalid number of rows of the right hand side matrix b"); } - + // return true if already allocated bool allocate(Index rows, Index cols, unsigned int computationOptions) ; MatrixUType m_matrixU; MatrixVType m_matrixV; SingularValuesType m_singularValues; + ComputationInfo m_info; bool m_isInitialized, m_isAllocated, m_usePrescribedThreshold; bool m_computeFullU, m_computeThinU; bool m_computeFullV, m_computeThinV; @@ -265,7 +284,8 @@ protected: * Default constructor of SVDBase */ SVDBase() - : m_isInitialized(false), + : m_info(Success), + m_isInitialized(false), m_isAllocated(false), m_usePrescribedThreshold(false), m_computeFullU(false), @@ -327,6 +347,7 @@ bool SVDBase::allocate(Index rows, Index cols, unsigned int computat m_rows = rows; m_cols = cols; + m_info = Success; m_isInitialized = false; m_isAllocated = true; m_computationOptions = computationOptions; -- cgit v1.2.3