bugfix in CholeskyWithoutSquareRoot::solve found by Timothy Hunter

This commit is contained in:
Gael Guennebaud 2008-08-16 07:09:39 +00:00
parent d9d69de382
commit 5546e963c9
3 changed files with 11 additions and 7 deletions

View File

@ -25,6 +25,8 @@ namespace Eigen {
* \endcode * \endcode
*/ */
#include "src/Array/CwiseOperators.h"
#include "src/Array/Functors.h"
#include "src/Cholesky/Cholesky.h" #include "src/Cholesky/Cholesky.h"
#include "src/Cholesky/CholeskyWithoutSquareRoot.h" #include "src/Cholesky/CholeskyWithoutSquareRoot.h"

View File

@ -149,10 +149,9 @@ typename Derived::Eval CholeskyWithoutSquareRoot<MatrixType>::solve(const Matrix
return m_matrix.adjoint().template part<UnitUpper>() return m_matrix.adjoint().template part<UnitUpper>()
.solveTriangular( .solveTriangular(
(matrixL() ( m_matrix.cwise().inverse().diagonal().asDiagonal()
.solveTriangular(b)) * matrixL().solveTriangular(b))
.cwise()/m_matrix.diagonal() );
);
} }
/** \cholesky_module /** \cholesky_module

View File

@ -39,16 +39,19 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
MatrixType a = MatrixType::Random(rows,cols); MatrixType a = MatrixType::Random(rows,cols);
VectorType b = VectorType::Random(rows); VectorType vecB = VectorType::Random(rows);
MatrixType matB = MatrixType::Random(rows,cols);
SquareMatrixType covMat = a * a.adjoint(); SquareMatrixType covMat = a * a.adjoint();
CholeskyWithoutSquareRoot<SquareMatrixType> cholnosqrt(covMat); CholeskyWithoutSquareRoot<SquareMatrixType> cholnosqrt(covMat);
VERIFY_IS_APPROX(covMat, cholnosqrt.matrixL() * cholnosqrt.vectorD().asDiagonal() * cholnosqrt.matrixL().adjoint()); VERIFY_IS_APPROX(covMat, cholnosqrt.matrixL() * cholnosqrt.vectorD().asDiagonal() * cholnosqrt.matrixL().adjoint());
VERIFY_IS_APPROX(covMat * cholnosqrt.solve(b), b); VERIFY_IS_APPROX(covMat * cholnosqrt.solve(vecB), vecB);
VERIFY_IS_APPROX(covMat * cholnosqrt.solve(matB), matB);
Cholesky<SquareMatrixType> chol(covMat); Cholesky<SquareMatrixType> chol(covMat);
VERIFY_IS_APPROX(covMat, chol.matrixL() * chol.matrixL().adjoint()); VERIFY_IS_APPROX(covMat, chol.matrixL() * chol.matrixL().adjoint());
VERIFY_IS_APPROX(covMat * chol.solve(b), b); VERIFY_IS_APPROX(covMat * chol.solve(vecB), vecB);
VERIFY_IS_APPROX(covMat * chol.solve(matB), matB);
} }
void test_cholesky() void test_cholesky()