12345qiupeng
/
gtsam
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Fixup cholesky decomposition - rename to LLt and RtR to make convention clear

release/4.3a0
justinca 2010-01-31 18:07:29 +00:00
parent 87cc573d34
commit 74dd7128ac
3 changed files with 14 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
cpp/Matrix.cpp
View File

@ -787,7 +787,7 @@ double** createNRC(Matrix& A) {
namespace BNU = boost::numeric::ublas; namespace BNU = boost::numeric::ublas;
Matrix cholesky(const Matrix& A) Matrix LLt(const Matrix& A)
{ {
assert(A.size1() == A.size2()); assert(A.size1() == A.size2());
Matrix L = zeros(A.size1(), A.size1()); Matrix L = zeros(A.size1(), A.size1());
@ -808,12 +808,17 @@ Matrix cholesky(const Matrix& A)
return L; return L;
} }
Matrix RtR(const Matrix &A)
{
return trans(LLt(A));
}
/* /*
* This is not ultra efficient, but not terrible, either. * This is not ultra efficient, but not terrible, either.
*/ */
Matrix cholesky_inverse(const Matrix &A) Matrix cholesky_inverse(const Matrix &A)
{ {
Matrix L = cholesky(A); Matrix L = LLt(A);
Matrix inv(boost::numeric::ublas::identity_matrix<double>(A.size1())); Matrix inv(boost::numeric::ublas::identity_matrix<double>(A.size1()));
inplace_solve (L, inv, BNU::lower_tag ()); inplace_solve (L, inv, BNU::lower_tag ());
return BNU::prod(trans(inv), inv); return BNU::prod(trans(inv), inv);
@ -883,7 +888,7 @@ Matrix inverse_square_root(const Matrix& A) {
// inv(B) * inv(B)' == A // inv(B) * inv(B)' == A
// inv(B' * B) == A // inv(B' * B) == A
Matrix inverse_square_root(const Matrix& A) { Matrix inverse_square_root(const Matrix& A) {
Matrix L = cholesky(A); Matrix L = LLt(A);
Matrix inv(boost::numeric::ublas::identity_matrix<double>(A.size1())); Matrix inv(boost::numeric::ublas::identity_matrix<double>(A.size1()));
inplace_solve (L, inv, BNU::lower_tag ()); inplace_solve (L, inv, BNU::lower_tag ());
return inv; return inv;

5
cpp/Matrix.h
View File

@ -321,7 +321,10 @@ Matrix inverse_square_root(const Matrix& A);
Matrix square_root_positive(const Matrix& A); Matrix square_root_positive(const Matrix& A);
/** Calculate the LL^t decomposition of a S.P.D matrix */ /** Calculate the LL^t decomposition of a S.P.D matrix */
Matrix cholesky(const Matrix& A); Matrix LLt(const Matrix& A);
/** Calculate the R^tR decomposition of a S.P.D matrix */
Matrix RtR(const Matrix& A);
/** Return the inverse of a S.P.D. matrix. Inversion is done via Cholesky decomposition. */ /** Return the inverse of a S.P.D. matrix. Inversion is done via Cholesky decomposition. */
Matrix cholesky_inverse(const Matrix &A); Matrix cholesky_inverse(const Matrix &A);

4
cpp/testMatrix.cpp
View File

@ -791,7 +791,7 @@ TEST( matrix, inverse_square_root )
/* *********************************************************************** */ /* *********************************************************************** */
// M was generated as the covariance of a set of random numbers. L that // M was generated as the covariance of a set of random numbers. L that
// we are checking against was generated via chol(M)' on octave // we are checking against was generated via chol(M)' on octave
TEST( matrix, cholesky ) TEST( matrix, LLt )
{ {
Matrix M = Matrix_(5, 5, Matrix M = Matrix_(5, 5,
0.0874197, -0.0030860, 0.0116969, 0.0081463, 0.0048741, 0.0874197, -0.0030860, 0.0116969, 0.0081463, 0.0048741,
@ -807,7 +807,7 @@ TEST( matrix, cholesky )
0.027552165831157, 0.043423266737274, 0.021695600982708, 0.267613525371710, 0.000000000000000, 0.027552165831157, 0.043423266737274, 0.021695600982708, 0.267613525371710, 0.000000000000000,
0.016485031422565, -0.012072546984405, -0.006621889326331, 0.014405837566082, 0.300462176944247); 0.016485031422565, -0.012072546984405, -0.006621889326331, 0.014405837566082, 0.300462176944247);
EQUALITY(expected, cholesky(M)); EQUALITY(expected, LLt(M));
} }
/* ************************************************************************* */ /* ************************************************************************* */