GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-08-14 04:35:57 +08:00
Code Issues Packages Projects Releases Wiki Activity

Make the SVD's output unitary and improved unit tests.

Browse Source
This commit is contained in:
Hauke Heibel 2010-09-24 16:28:20 +02:00
parent 1c54514bfc
commit 053261de88
2 changed files with 19 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -445,6 +445,19 @@ SVD<MatrixType>& SVD<MatrixType>::compute(const MatrixType& matrix)
m_matU.leftCols(n) = A; m_matU.leftCols(n) = A;
m_matU.rightCols(m-n).setZero(); m_matU.rightCols(m-n).setZero();
// Gram Schmidt orthogonalization to fill up U
for (int col = A.cols(); col < A.rows(); ++col)
{
typename MatrixUType::ColXpr colVec = m_matU.col(col);
colVec(col) = 1;
for (int prevCol = 0; prevCol < col; ++prevCol)
{
typename MatrixUType::ColXpr prevColVec = m_matU.col(prevCol);
colVec -= colVec.dot(prevColVec)*prevColVec;
}
m_matU.col(col) = colVec.normalized();
}
m_isInitialized = true; m_isInitialized = true;
return *this; return *this;
} }

12
test/svd.cpp
View File

@ -47,9 +47,9 @@ template<typename MatrixType> void svd(const MatrixType& m)
sigma.diagonal() = svd.singularValues(); sigma.diagonal() = svd.singularValues();
matU = svd.matrixU(); matU = svd.matrixU();
//VERIFY_IS_UNITARY(matU); VERIFY_IS_UNITARY(matU);
matV = svd.matrixV(); matV = svd.matrixV();
//VERIFY_IS_UNITARY(matV); VERIFY_IS_UNITARY(matV);
VERIFY_IS_APPROX(a, matU * sigma * matV.transpose()); VERIFY_IS_APPROX(a, matU * sigma * matV.transpose());
} }
@ -57,10 +57,10 @@ template<typename MatrixType> void svd(const MatrixType& m)
if (rows>=cols) if (rows>=cols)
{ {
SVD<MatrixType> svd(a); SVD<MatrixType> svd(a);
Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> x = Matrix<Scalar, MatrixType::ColsAtCompileTime, 1>::Random(cols,1); Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> b = Matrix<Scalar, MatrixType::ColsAtCompileTime, 1>::Random(rows,1);
Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> b = a * x; Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> x = svd.solve(b);
Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> result = svd.solve(b); // evaluate normal equation which works also for least-squares solutions
VERIFY_IS_APPROX(a * result, b); VERIFY_IS_APPROX(a.adjoint()*a*x,a.adjoint()*b);
} }