mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-03-07 18:27:40 +08:00
Code Issues Packages Projects Releases Wiki Activity

D&C SVD: add scaling to avoid overflow, fix handling of fixed size matrices

Browse Source
This commit is contained in:
Gael Guennebaud 2014-10-06 19:35:57 +02:00
parent d44d432baa
commit dbdd8b0883
3 changed files with 28 additions and 23 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
test/svd_common.h
View File

@ -301,7 +301,7 @@ void svd_underoverflow()
0, 5.60844e-313;
SVD_DEFAULT(Matrix2d) svd;
svd.compute(M,ComputeFullU|ComputeFullV);
svd_check_full(M,svd);
CALL_SUBTEST( svd_check_full(M,svd) );
// Check all 2x2 matrices made with the following coefficients:
VectorXd value_set(9);
@ -312,7 +312,7 @@ void svd_underoverflow()
{
M << value_set(id(0)), value_set(id(1)), value_set(id(2)), value_set(id(3));
svd.compute(M,ComputeFullU|ComputeFullV);
svd_check_full(M,svd);
CALL_SUBTEST( svd_check_full(M,svd) );
id(k)++;
if(id(k)>=value_set.size())
@ -336,7 +336,7 @@ void svd_underoverflow()
SVD_DEFAULT(Matrix3d) svd3;
svd3.compute(M3,ComputeFullU|ComputeFullV); // just check we don't loop indefinitely
svd_check_full(M3,svd3);
CALL_SUBTEST( svd_check_full(M3,svd3) );
}
// void jacobisvd(const MatrixType& a = MatrixType(), bool pickrandom = true)

35
unsupported/Eigen/src/BDCSVD/BDCSVD.h
View File

@ -236,26 +236,29 @@ BDCSVD<MatrixType>& BDCSVD<MatrixType>::compute(const MatrixType& matrix, unsign
{
allocate(matrix.rows(), matrix.cols(), computationOptions);
using std::abs;
//**** step 1 Bidiagonalization
MatrixType copy;
if (m_isTranspose) copy = matrix.adjoint();
else copy = matrix;
//**** step 0 - Copy the input matrix and apply scaling to reduce over/under-flows
RealScalar scale = matrix.cwiseAbs().maxCoeff();
if(scale==RealScalar(0)) scale = RealScalar(1);
MatrixX copy;
if (m_isTranspose) copy = matrix.adjoint()/scale;
else copy = matrix/scale;
//**** step 1 - Bidiagonalization
internal::UpperBidiagonalization<MatrixX> bid(copy);
//**** step 2 Divide
//**** step 2 - Divide & Conquer
m_naiveU.setZero();
m_naiveV.setZero();
m_computed.topRows(m_diagSize) = bid.bidiagonal().toDenseMatrix().transpose();
m_computed.template bottomRows<1>().setZero();
divide(0, m_diagSize - 1, 0, 0, 0);
//**** step 3 copy
//**** step 3 - Copy singular values and vectors
for (int i=0; i<m_diagSize; i++)
{
RealScalar a = abs(m_computed.coeff(i, i));
m_singularValues.coeffRef(i) = a;
m_singularValues.coeffRef(i) = a * scale;
if (a == 0)
{
m_nonzeroSingularValues = i;
@ -268,11 +271,13 @@ BDCSVD<MatrixType>& BDCSVD<MatrixType>::compute(const MatrixType& matrix, unsign
break;
}
}
// std::cout << "m_naiveU\n" << m_naiveU << "\n\n";
// std::cout << "m_naiveV\n" << m_naiveV << "\n\n";
#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
std::cout << "m_naiveU\n" << m_naiveU << "\n\n";
std::cout << "m_naiveV\n" << m_naiveV << "\n\n";
#endif
if(m_isTranspose) copyUV(bid.householderV(), bid.householderU(), m_naiveV, m_naiveU);
else copyUV(bid.householderU(), bid.householderV(), m_naiveU, m_naiveV);
// std::cout << "DONE\n";
m_isInitialized = true;
return *this;
}// end compute
@ -569,13 +574,13 @@ void BDCSVD<MatrixType>::computeSVDofM(Index firstCol, Index n, MatrixXr& U, Vec
computeSingVecs(zhat, diag, singVals, shifts, mus, U, V);
#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
std::cout << "U^T U: " << (U.transpose() * U - MatrixType(MatrixType::Identity(U.cols(),U.cols()))).norm() << "\n";
std::cout << "V^T V: " << (V.transpose() * V - MatrixType(MatrixType::Identity(V.cols(),V.cols()))).norm() << "\n";
std::cout << "U^T U: " << (U.transpose() * U - MatrixXr(MatrixXr::Identity(U.cols(),U.cols()))).norm() << "\n";
std::cout << "V^T V: " << (V.transpose() * V - MatrixXr(MatrixXr::Identity(V.cols(),V.cols()))).norm() << "\n";
#endif
#ifdef EIGEN_BDCSVD_SANITY_CHECKS
assert((U.transpose() * U - MatrixType(MatrixType::Identity(U.cols(),U.cols()))).norm() < 1e-14 * n);
assert((V.transpose() * V - MatrixType(MatrixType::Identity(V.cols(),V.cols()))).norm() < 1e-14 * n);
assert((U.transpose() * U - MatrixXr(MatrixXr::Identity(U.cols(),U.cols()))).norm() < 1e-14 * n);
assert((V.transpose() * V - MatrixXr(MatrixXr::Identity(V.cols(),V.cols()))).norm() < 1e-14 * n);
assert(m_naiveU.allFinite());
assert(m_naiveV.allFinite());
assert(m_computed.allFinite());

10
unsupported/test/bdcsvd.cpp
View File

@ -70,13 +70,13 @@ void test_bdcsvd()
CALL_SUBTEST_7(( svd_verify_assert<BDCSVD<MatrixXf> >(MatrixXf(10,12)) ));
CALL_SUBTEST_8(( svd_verify_assert<BDCSVD<MatrixXcd> >(MatrixXcd(7,5)) ));
// svd_all_trivial_2x2(bdcsvd<Matrix2cd>);
// svd_all_trivial_2x2(bdcsvd<Matrix2d>);
CALL_SUBTEST_1(( svd_all_trivial_2x2(bdcsvd<Matrix2cd>) ));
CALL_SUBTEST_1(( svd_all_trivial_2x2(bdcsvd<Matrix2d>) ));
for(int i = 0; i < g_repeat; i++) {
// CALL_SUBTEST_3(( bdcsvd<Matrix3f>() ));
// CALL_SUBTEST_4(( bdcsvd<Matrix4d>() ));
// CALL_SUBTEST_5(( bdcsvd<Matrix<float,3,5> >() ));
CALL_SUBTEST_3(( bdcsvd<Matrix3f>() ));
CALL_SUBTEST_4(( bdcsvd<Matrix4d>() ));
CALL_SUBTEST_5(( bdcsvd<Matrix<float,3,5> >() ));
int r = internal::random<int>(1, EIGEN_TEST_MAX_SIZE/2),
c = internal::random<int>(1, EIGEN_TEST_MAX_SIZE/2);