mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-01-12 14:25:16 +08:00
Code Issues Packages Projects Releases Wiki Activity

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.

Browse Source 
Make Jacobi and divide-and-conquer fail fast and return info() == InvalidInput if the matrix contains NaN or +/-Inf.
This commit is contained in:
Rasmus Munk Larsen 2021-03-31 21:09:19 +00:00
parent b5a926a0f6
commit 5bbc9cea93
4 changed files with 72 additions and 19 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

32
Eigen/src/SVD/BDCSVD.h
View File

@ -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<MatrixType>& BDCSVD<MatrixType>::compute(const MatrixType& matrix, unsign
{
// 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 @@ BDCSVD<MatrixType>& BDCSVD<MatrixType>::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<m_diagSize; i++)
{
@ -394,7 +408,7 @@ void BDCSVD<MatrixType>::structured_update(Block<MatrixXr,Dynamic,Dynamic> 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<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 @@ void BDCSVD<MatrixType>::divide (Eigen::Index firstCol, Eigen::Index lastCol, Ei
{
// 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 @@ void BDCSVD<MatrixType>::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)
{

9
Eigen/src/SVD/JacobiSVD.h
View File

@ -585,6 +585,7 @@ template<typename _MatrixType, int QRPreconditioner> 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<MatrixType, QRPreconditioner>::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<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsig
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 */

41
Eigen/src/SVD/SVDBase.h
View File

@ -52,7 +52,7 @@ template<typename Derived> struct traits<SVDBase<Derived> >
* 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<RealScalar>(m_singularValues.coeff(0) * threshold(), (std::numeric_limits<RealScalar>::min)());
Index i = m_nonzeroSingularValues-1;
@ -224,6 +224,18 @@ public:
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 @@ 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<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 @@ 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<MatrixType>::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;

9
test/svd_common.h
View File

@ -298,7 +298,8 @@ EIGEN_DONT_INLINE Scalar zero() { return Scalar(0); }
// 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 @@ void svd_inf_nan()
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 @@ void svd_inf_nan()
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 @@ void svd_inf_nan()
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