devtools/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2026-04-10 11:34:33 +08:00
Code Issues Packages Projects Releases Wiki Activity

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

Browse Source 
(grafted from 019dcfc21d
)
This commit is contained in:
Gael Guennebaud
2013-11-03 13:18:56 +01:00
parent e16e52d493
commit 0c4fc69d62
2 changed files with 59 additions and 6 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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>