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

JacobiSVD: move from Lapack to Matlab strategy for the default threshold

Browse Source
This commit is contained in:
Gael Guennebaud 2013-11-03 13:18:56 +01:00
parent 19521c83b8
commit 019dcfc21d
2 changed files with 59 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -530,7 +530,7 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
m_isAllocated(false),
m_usePrescribedThreshold(false),
m_computationOptions(0),
m_rows(-1), m_cols(-1)
m_rows(-1), m_cols(-1), m_diagSize(0)
{}
@ -726,7 +726,7 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
{
eigen_assert(m_isInitialized || m_usePrescribedThreshold);
return m_usePrescribedThreshold ? m_prescribedThreshold
: NumTraits<Scalar>::epsilon();
: (std::max<Index>)(1,m_diagSize)*NumTraits<Scalar>::epsilon();
}
inline Index rows() const { return m_rows; }
@ -918,11 +918,11 @@ struct solve_retval<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs>
// So A^{-1} = V S^{-1} U^*
Matrix<Scalar, Dynamic, Rhs::ColsAtCompileTime, 0, _MatrixType::MaxRowsAtCompileTime, Rhs::MaxColsAtCompileTime> tmp;
Index nonzeroSingVals = dec().rank();
Index rank = dec().rank();
tmp.noalias() = dec().matrixU().leftCols(nonzeroSingVals).adjoint() * rhs();
tmp = dec().singularValues().head(nonzeroSingVals).asDiagonal().inverse() * tmp;
dst = dec().matrixV().leftCols(nonzeroSingVals) * tmp;
tmp.noalias() = dec().matrixU().leftCols(rank).adjoint() * rhs();
tmp = dec().singularValues().head(rank).asDiagonal().inverse() * tmp;
dst = dec().matrixV().leftCols(rank) * tmp;
}
};
} // end namespace internal

53
test/jacobisvd.cpp
View File

@ -84,6 +84,59 @@ void jacobisvd_solve(const MatrixType& m, unsigned int computationOptions)
SolutionType x = svd.solve(rhs);
// evaluate normal equation which works also for least-squares solutions
VERIFY_IS_APPROX(m.adjoint()*m*x,m.adjoint()*rhs);
// check minimal norm solutions
{
// generate a full-rank m x n problem with m<n
enum {
RankAtCompileTime2 = ColsAtCompileTime==Dynamic ? Dynamic : (ColsAtCompileTime)/2+1,
RowsAtCompileTime3 = ColsAtCompileTime==Dynamic ? Dynamic : ColsAtCompileTime+1
};
typedef Matrix<Scalar, RankAtCompileTime2, ColsAtCompileTime> MatrixType2;
typedef Matrix<Scalar, RankAtCompileTime2, 1> RhsType2;
typedef Matrix<Scalar, ColsAtCompileTime, RankAtCompileTime2> MatrixType2T;
Index rank = RankAtCompileTime2==Dynamic ? internal::random<Index>(1,cols) : Index(RankAtCompileTime2);
MatrixType2 m2(rank,cols);
int guard = 0;
do {
m2.setRandom();
} while(m2.jacobiSvd().setThreshold(test_precision<Scalar>()).rank()!=rank && (++guard)<10);
VERIFY(guard<10);
RhsType2 rhs2 = RhsType2::Random(rank);
// use QR to find a reference minimal norm solution
HouseholderQR<MatrixType2T> qr(m2.adjoint());
Matrix<Scalar,Dynamic,1> tmp = qr.matrixQR().topLeftCorner(rank,rank).template triangularView<Upper>().adjoint().solve(rhs2);
tmp.conservativeResize(cols);
tmp.tail(cols-rank).setZero();
SolutionType x21 = qr.householderQ() * tmp;
// now check with SVD
JacobiSVD<MatrixType2, ColPivHouseholderQRPreconditioner> svd2(m2, computationOptions);
SolutionType x22 = svd2.solve(rhs2);
VERIFY_IS_APPROX(m2*x21, rhs2);
VERIFY_IS_APPROX(m2*x22, rhs2);
VERIFY_IS_APPROX(x21, x22);
// Now check with a rank deficient matrix
typedef Matrix<Scalar, RowsAtCompileTime3, ColsAtCompileTime> MatrixType3;
typedef Matrix<Scalar, RowsAtCompileTime3, 1> RhsType3;
Index rows3 = RowsAtCompileTime3==Dynamic ? internal::random<Index>(rank+1,2*cols) : Index(RowsAtCompileTime3);
Matrix<Scalar,RowsAtCompileTime3,Dynamic> C = Matrix<Scalar,RowsAtCompileTime3,Dynamic>::Random(rows3,rank);
MatrixType3 m3 = C * m2;
RhsType3 rhs3 = C * rhs2;
JacobiSVD<MatrixType3, ColPivHouseholderQRPreconditioner> svd3(m3, computationOptions);
SolutionType x3 = svd3.solve(rhs3);
if(svd3.rank()!=rank) {
std::cout << m3 << "\n\n";
std::cout << svd3.singularValues().transpose() << "\n";
std::cout << svd3.rank() << " == " << rank << "\n";
std::cout << x21.norm() << " == " << x3.norm() << "\n";
}
// VERIFY_IS_APPROX(m3*x3, rhs3);
VERIFY_IS_APPROX(m3*x21, rhs3);
VERIFY_IS_APPROX(m2*x3, rhs2);
VERIFY_IS_APPROX(x21, x3);
}
}
template<typename MatrixType, int QRPreconditioner>