mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2024-12-21 07:19:46 +08:00
Code Issues Packages Projects Releases Wiki Activity

bug #843: fix jacobisvd for complexes and extend respective unit test to chack with random tricky matrices,

Browse Source 
and backport other JacobiSVD fixes
This commit is contained in:
Gael Guennebaud 2014-07-17 17:09:15 +02:00
parent e215740e8e
commit 9b00035438
2 changed files with 86 additions and 28 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -375,17 +375,19 @@ struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, true>
Scalar z;
JacobiRotation<Scalar> rot;
RealScalar n = sqrt(numext::abs2(work_matrix.coeff(p,p)) + numext::abs2(work_matrix.coeff(q,p)));
if(n==0)
{
z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q);
work_matrix.row(p) *= z;
if(svd.computeU()) svd.m_matrixU.col(p) *= conj(z);
if(work_matrix.coeff(q,q)!=Scalar(0))
{
z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q);
else
z = Scalar(0);
work_matrix.row(q) *= z;
if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z);
work_matrix.row(q) *= z;
if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z);
}
// otherwise the second row is already zero, so we have nothing to do.
}
else
{
@ -415,6 +417,7 @@ void real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q,
JacobiRotation<RealScalar> *j_right)
{
using std::sqrt;
using std::abs;
Matrix<RealScalar,2,2> m;
m << numext::real(matrix.coeff(p,p)), numext::real(matrix.coeff(p,q)),
numext::real(matrix.coeff(q,p)), numext::real(matrix.coeff(q,q));
@ -428,9 +431,11 @@ void real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q,
}
else
{
RealScalar u = d / t;
rot1.c() = RealScalar(1) / sqrt(RealScalar(1) + numext::abs2(u));
rot1.s() = rot1.c() * u;
RealScalar t2d2 = numext::hypot(t,d);
rot1.c() = abs(t)/t2d2;
rot1.s() = d/t2d2;
if(t<RealScalar(0))
rot1.s() = -rot1.s();
}
m.applyOnTheLeft(0,1,rot1);
j_right->makeJacobi(m,0,1);
@ -831,6 +836,11 @@ JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsig
if(m_computeFullV) m_matrixV.setIdentity(m_cols,m_cols);
if(m_computeThinV) m_matrixV.setIdentity(m_cols, m_diagSize);
}
// Scaling factor to reduce over/under-flows
RealScalar scale = m_workMatrix.cwiseAbs().maxCoeff();
if(scale==RealScalar(0)) scale = RealScalar(1);
m_workMatrix /= scale;
/*** step 2. The main Jacobi SVD iteration. ***/
@ -900,6 +910,8 @@ JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsig
if(computeV()) m_matrixV.col(pos).swap(m_matrixV.col(i));
}
}
m_singularValues *= scale;
m_isInitialized = true;
return *this;

88
test/jacobisvd.cpp
View File

@ -67,6 +67,7 @@ template<typename MatrixType, int QRPreconditioner>
void jacobisvd_solve(const MatrixType& m, unsigned int computationOptions)
{
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
typedef typename MatrixType::Index Index;
Index rows = m.rows();
Index cols = m.cols();
@ -81,9 +82,37 @@ void jacobisvd_solve(const MatrixType& m, unsigned int computationOptions)
RhsType rhs = RhsType::Random(rows, internal::random<Index>(1, cols));
JacobiSVD<MatrixType, QRPreconditioner> svd(m, computationOptions);
if(internal::is_same<RealScalar,double>::value) svd.setThreshold(1e-8);
else if(internal::is_same<RealScalar,float>::value) svd.setThreshold(1e-4);
SolutionType x = svd.solve(rhs);
RealScalar residual = (m*x-rhs).norm();
// Check that there is no significantly better solution in the neighborhood of x
if(!test_isMuchSmallerThan(residual,rhs.norm()))
{
// If the residual is very small, then we have an exact solution, so we are already good.
for(int k=0;k<x.rows();++k)
{
SolutionType y(x);
y.row(k).array() += 2*NumTraits<RealScalar>::epsilon();
RealScalar residual_y = (m*y-rhs).norm();
VERIFY( test_isApprox(residual_y,residual) || residual < residual_y );
y.row(k) = x.row(k).array() - 2*NumTraits<RealScalar>::epsilon();
residual_y = (m*y-rhs).norm();
VERIFY( test_isApprox(residual_y,residual) || residual < residual_y );
}
}
// evaluate normal equation which works also for least-squares solutions
VERIFY_IS_APPROX(m.adjoint()*m*x,m.adjoint()*rhs);
if(internal::is_same<RealScalar,double>::value)
{
// This test is not stable with single precision.
// This is probably because squaring m signicantly affects the precision.
VERIFY_IS_APPROX(m.adjoint()*m*x,m.adjoint()*rhs);
}
// check minimal norm solutions
{
@ -145,10 +174,9 @@ void jacobisvd_test_all_computation_options(const MatrixType& m)
if (QRPreconditioner == NoQRPreconditioner && m.rows() != m.cols())
return;
JacobiSVD<MatrixType, QRPreconditioner> fullSvd(m, ComputeFullU|ComputeFullV);
jacobisvd_check_full(m, fullSvd);
jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeFullU | ComputeFullV);
CALL_SUBTEST(( jacobisvd_check_full(m, fullSvd) ));
CALL_SUBTEST(( jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeFullU | ComputeFullV) ));
#if defined __INTEL_COMPILER
// remark #111: statement is unreachable
#pragma warning disable 111
@ -156,20 +184,20 @@ void jacobisvd_test_all_computation_options(const MatrixType& m)
if(QRPreconditioner == FullPivHouseholderQRPreconditioner)
return;
jacobisvd_compare_to_full(m, ComputeFullU, fullSvd);
jacobisvd_compare_to_full(m, ComputeFullV, fullSvd);
jacobisvd_compare_to_full(m, 0, fullSvd);
CALL_SUBTEST(( jacobisvd_compare_to_full(m, ComputeFullU, fullSvd) ));
CALL_SUBTEST(( jacobisvd_compare_to_full(m, ComputeFullV, fullSvd) ));
CALL_SUBTEST(( jacobisvd_compare_to_full(m, 0, fullSvd) ));
if (MatrixType::ColsAtCompileTime == Dynamic) {
// thin U/V are only available with dynamic number of columns
jacobisvd_compare_to_full(m, ComputeFullU|ComputeThinV, fullSvd);
jacobisvd_compare_to_full(m, ComputeThinV, fullSvd);
jacobisvd_compare_to_full(m, ComputeThinU|ComputeFullV, fullSvd);
jacobisvd_compare_to_full(m, ComputeThinU , fullSvd);
jacobisvd_compare_to_full(m, ComputeThinU|ComputeThinV, fullSvd);
jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeFullU | ComputeThinV);
jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeThinU | ComputeFullV);
jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeThinU | ComputeThinV);
CALL_SUBTEST(( jacobisvd_compare_to_full(m, ComputeFullU|ComputeThinV, fullSvd) ));
CALL_SUBTEST(( jacobisvd_compare_to_full(m, ComputeThinV, fullSvd) ));
CALL_SUBTEST(( jacobisvd_compare_to_full(m, ComputeThinU|ComputeFullV, fullSvd) ));
CALL_SUBTEST(( jacobisvd_compare_to_full(m, ComputeThinU , fullSvd) ));
CALL_SUBTEST(( jacobisvd_compare_to_full(m, ComputeThinU|ComputeThinV, fullSvd) ));
CALL_SUBTEST(( jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeFullU | ComputeThinV) ));
CALL_SUBTEST(( jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeThinU | ComputeFullV) ));
CALL_SUBTEST(( jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeThinU | ComputeThinV) ));
// test reconstruction
typedef typename MatrixType::Index Index;
@ -182,12 +210,29 @@ void jacobisvd_test_all_computation_options(const MatrixType& m)
template<typename MatrixType>
void jacobisvd(const MatrixType& a = MatrixType(), bool pickrandom = true)
{
MatrixType m = pickrandom ? MatrixType::Random(a.rows(), a.cols()) : a;
MatrixType m = a;
if(pickrandom)
{
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
typedef typename MatrixType::Index Index;
Index diagSize = (std::min)(a.rows(), a.cols());
RealScalar s = std::numeric_limits<RealScalar>::max_exponent10/4;
s = internal::random<RealScalar>(1,s);
Matrix<RealScalar,Dynamic,1> d = Matrix<RealScalar,Dynamic,1>::Random(diagSize);
for(Index k=0; k<diagSize; ++k)
d(k) = d(k)*std::pow(RealScalar(10),internal::random<RealScalar>(-s,s));
m = Matrix<Scalar,Dynamic,Dynamic>::Random(a.rows(),diagSize) * d.asDiagonal() * Matrix<Scalar,Dynamic,Dynamic>::Random(diagSize,a.cols());
// cancel some coeffs
Index n = internal::random<Index>(0,m.size()-1);
for(Index i=0; i<n; ++i)
m(internal::random<Index>(0,m.rows()-1), internal::random<Index>(0,m.cols()-1)) = Scalar(0);
}
jacobisvd_test_all_computation_options<MatrixType, FullPivHouseholderQRPreconditioner>(m);
jacobisvd_test_all_computation_options<MatrixType, ColPivHouseholderQRPreconditioner>(m);
jacobisvd_test_all_computation_options<MatrixType, HouseholderQRPreconditioner>(m);
jacobisvd_test_all_computation_options<MatrixType, NoQRPreconditioner>(m);
CALL_SUBTEST(( jacobisvd_test_all_computation_options<MatrixType, FullPivHouseholderQRPreconditioner>(m) ));
CALL_SUBTEST(( jacobisvd_test_all_computation_options<MatrixType, ColPivHouseholderQRPreconditioner>(m) ));
CALL_SUBTEST(( jacobisvd_test_all_computation_options<MatrixType, HouseholderQRPreconditioner>(m) ));
CALL_SUBTEST(( jacobisvd_test_all_computation_options<MatrixType, NoQRPreconditioner>(m) ));
}
template<typename MatrixType> void jacobisvd_verify_assert(const MatrixType& m)
@ -381,6 +426,7 @@ void test_jacobisvd()
TEST_SET_BUT_UNUSED_VARIABLE(r)
TEST_SET_BUT_UNUSED_VARIABLE(c)
CALL_SUBTEST_10(( jacobisvd<MatrixXd>(MatrixXd(r,c)) ));
CALL_SUBTEST_7(( jacobisvd<MatrixXf>(MatrixXf(r,c)) ));
CALL_SUBTEST_8(( jacobisvd<MatrixXcd>(MatrixXcd(r,c)) ));
(void) r;