diff options
author | Gael Guennebaud <g.gael@free.fr> | 2013-11-03 13:18:56 +0100 |
---|---|---|
committer | Gael Guennebaud <g.gael@free.fr> | 2013-11-03 13:18:56 +0100 |
commit | 019dcfc21dc74f0c735d72918731b5349a62a26a (patch) | |
tree | c0eec562f65934030e25453427e6c32b21cb7bac /test | |
parent | 19521c83b8d5c9b0567ceadebe7ad4b1119dd7cc (diff) |
JacobiSVD: move from Lapack to Matlab strategy for the default threshold
Diffstat (limited to 'test')
-rw-r--r-- | test/jacobisvd.cpp | 53 |
1 files changed, 53 insertions, 0 deletions
diff --git a/test/jacobisvd.cpp b/test/jacobisvd.cpp index 76157c30f..e378de477 100644 --- a/test/jacobisvd.cpp +++ b/test/jacobisvd.cpp @@ -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> |