mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-01-24 14:45:14 +08:00
Code Issues Packages Projects Releases Wiki Activity

SVD:

Browse Source 
* fix unit test for rectangular matrices.
 * enforce that rows >= cols since various places in the code assume that.
This commit is contained in:
Benoit Jacob 2010-09-23 09:51:08 -04:00
parent 62bf04b339
commit c253cc3d53
2 changed files with 19 additions and 17 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
Eigen/src/SVD/SVD.h
View File

@ -38,6 +38,8 @@ template<typename MatrixType, typename Rhs> struct ei_svd_solve_impl;
* *
* This class performs a standard SVD decomposition of a real matrix A of size \c M x \c N. * This class performs a standard SVD decomposition of a real matrix A of size \c M x \c N.
* *
* Requires M >= N, in other words, at least as many rows as columns.
*
* \sa MatrixBase::SVD() * \sa MatrixBase::SVD()
*/ */
template<typename _MatrixType> class SVD template<typename _MatrixType> class SVD
@ -440,11 +442,8 @@ SVD<MatrixType>& SVD<MatrixType>::compute(const MatrixType& matrix)
} }
} }
} }
m_matU.setZero(); m_matU.leftCols(n) = A;
if (m>=n) m_matU.rightCols(m-n).setZero();
m_matU.block(0,0,m,n) = A;
else
m_matU = A.block(0,0,m,m);
m_isInitialized = true; m_isInitialized = true;
return *this; return *this;

27
test/svd.cpp
View File

@ -37,30 +37,30 @@ template<typename MatrixType> void svd(const MatrixType& m)
typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<Scalar>::Real RealScalar; typedef typename NumTraits<Scalar>::Real RealScalar;
MatrixType a = MatrixType::Random(rows,cols); MatrixType a = MatrixType::Random(rows,cols);
Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> b =
Matrix<Scalar, MatrixType::RowsAtCompileTime, 1>::Random(rows,1);
Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> x(cols,1), x2(cols,1); Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> x(cols,1), x2(cols,1);
{ {
SVD<MatrixType> svd(a); SVD<MatrixType> svd(a);
MatrixType sigma = MatrixType::Zero(rows,cols); MatrixType sigma = MatrixType::Zero(rows,cols);
MatrixType matU = MatrixType::Zero(rows,rows); MatrixType matU = MatrixType::Zero(rows,rows);
MatrixType matV = MatrixType::Zero(cols,cols);
sigma.diagonal() = svd.singularValues(); sigma.diagonal() = svd.singularValues();
matU = svd.matrixU(); matU = svd.matrixU();
VERIFY_IS_APPROX(a, matU * sigma * svd.matrixV().transpose()); //VERIFY_IS_UNITARY(matU);
matV = svd.matrixV();
//VERIFY_IS_UNITARY(matV);
VERIFY_IS_APPROX(a, matU * sigma * matV.transpose());
} }
if (rows>=cols) if (rows>=cols)
{ {
if (ei_is_same_type<RealScalar,float>::ret)
{
MatrixType a1 = MatrixType::Random(rows,cols);
a += a * a.adjoint() + a1 * a1.adjoint();
}
SVD<MatrixType> svd(a); SVD<MatrixType> svd(a);
x = svd.solve(b); Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> x = Matrix<Scalar, MatrixType::ColsAtCompileTime, 1>::Random(cols,1);
VERIFY_IS_APPROX(a * x,b); Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> b = a * x;
Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> result = svd.solve(b);
VERIFY_IS_APPROX(a * result, b);
} }
@ -102,8 +102,11 @@ void test_svd()
for(int i = 0; i < g_repeat; i++) { for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST_1( svd(Matrix3f()) ); CALL_SUBTEST_1( svd(Matrix3f()) );
CALL_SUBTEST_2( svd(Matrix4d()) ); CALL_SUBTEST_2( svd(Matrix4d()) );
CALL_SUBTEST_3( svd(MatrixXf(7,7)) ); int cols = ei_random<int>(2,50);
CALL_SUBTEST_4( svd(MatrixXd(14,7)) ); int rows = cols + ei_random<int>(0,50);
CALL_SUBTEST_3( svd(MatrixXf(rows,cols)) );
CALL_SUBTEST_4( svd(MatrixXd(rows,cols)) );
// complex are not implemented yet // complex are not implemented yet
// CALL_SUBTEST(svd(MatrixXcd(6,6)) ); // CALL_SUBTEST(svd(MatrixXcd(6,6)) );
// CALL_SUBTEST(svd(MatrixXcf(3,3)) ); // CALL_SUBTEST(svd(MatrixXcf(3,3)) );