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.
diff --git a/Eigen/src/SVD/BDCSVD.h b/Eigen/src/SVD/BDCSVD.h index e0c4456..17f8e44 100644 --- a/Eigen/src/SVD/BDCSVD.h +++ b/Eigen/src/SVD/BDCSVD.h
@@ -208,6 +208,7 @@ 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 @@ { // FIXME this line involves temporaries JacobiSVD<MatrixType> 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<PropagateNaN>(); + 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 @@ 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<m_diagSize; i++) { @@ -394,7 +408,7 @@ //@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<typename MatrixType> -void BDCSVD<MatrixType>::divide (Eigen::Index firstCol, Eigen::Index lastCol, Eigen::Index firstRowW, Eigen::Index firstColW, Eigen::Index shift) +void BDCSVD<MatrixType>::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 @@ { // FIXME this line involves temporaries JacobiSVD<MatrixXr> 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 @@ // 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 2b68911..a22a2e5 100644 --- a/Eigen/src/SVD/JacobiSVD.h +++ b/Eigen/src/SVD/JacobiSVD.h
@@ -585,6 +585,7 @@ 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 @@ m_rows = rows; m_cols = cols; + m_info = Success; m_isInitialized = false; m_isAllocated = true; m_computationOptions = computationOptions; @@ -674,7 +676,12 @@ const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)(); // Scaling factor to reduce over/under-flows - RealScalar scale = matrix.cwiseAbs().maxCoeff(); + RealScalar scale = matrix.cwiseAbs().template maxCoeff<PropagateNaN>(); + 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 34d5c9d..13c8394 100644 --- a/Eigen/src/SVD/SVDBase.h +++ b/Eigen/src/SVD/SVDBase.h
@@ -52,7 +52,7 @@ * 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 @@ */ 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 @@ */ 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 @@ */ 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 @@ 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<RealScalar>(m_singularValues.coeff(0) * threshold(), (std::numeric_limits<RealScalar>::min)()); Index i = m_nonzeroSingularValues-1; @@ -224,6 +224,18 @@ solve(const MatrixBase<Rhs>& 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<typename RhsType, typename DstType> void _solve_impl(const RhsType &rhs, DstType &dst) const; @@ -233,26 +245,33 @@ #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<bool Transpose_, typename Rhs> 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 @@ * 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 @@ m_rows = rows; m_cols = cols; + m_info = Success; m_isInitialized = false; m_isAllocated = true; m_computationOptions = computationOptions;
diff --git a/test/svd_common.h b/test/svd_common.h index 5c0f2a0..bd62edc 100644 --- a/test/svd_common.h +++ b/test/svd_common.h
@@ -298,7 +298,8 @@ // workaround aggressive optimization in ICC template<typename T> EIGEN_DONT_INLINE T sub(T a, T b) { return a - b; } -// all this function does is verify we don't iterate infinitely on nan/inf values +// This function verifies we don't iterate infinitely on nan/inf values, +// and that info() returns InvalidInput. template<typename SvdType, typename MatrixType> void svd_inf_nan() { @@ -307,18 +308,22 @@ Scalar some_inf = Scalar(1) / zero<Scalar>(); VERIFY(sub(some_inf, some_inf) != sub(some_inf, some_inf)); svd.compute(MatrixType::Constant(10,10,some_inf), ComputeFullU | ComputeFullV); + VERIFY(svd.info() == InvalidInput); Scalar nan = std::numeric_limits<Scalar>::quiet_NaN(); VERIFY(nan != nan); svd.compute(MatrixType::Constant(10,10,nan), ComputeFullU | ComputeFullV); + VERIFY(svd.info() == InvalidInput); MatrixType m = MatrixType::Zero(10,10); m(internal::random<int>(0,9), internal::random<int>(0,9)) = some_inf; svd.compute(m, ComputeFullU | ComputeFullV); + VERIFY(svd.info() == InvalidInput); m = MatrixType::Zero(10,10); m(internal::random<int>(0,9), internal::random<int>(0,9)) = nan; svd.compute(m, ComputeFullU | ComputeFullV); + VERIFY(svd.info() == InvalidInput); // regression test for bug 791 m.resize(3,3); @@ -326,6 +331,7 @@ 0, -0.5, 0, nan, 0, 0; svd.compute(m, ComputeFullU | ComputeFullV); + VERIFY(svd.info() == InvalidInput); m.resize(4,4); m << 1, 0, 0, 0, @@ -333,6 +339,7 @@ 1, 0, 1, nan, 0, nan, nan, 0; svd.compute(m, ComputeFullU | ComputeFullV); + VERIFY(svd.info() == InvalidInput); } // Regression test for bug 286: JacobiSVD loops indefinitely with some