From 74dd7128ac4f404e95c3bcf73c0fe14450c93d2b Mon Sep 17 00:00:00 2001
From: justinca <justinca@noemail>
Date: Sun, 31 Jan 2010 18:07:29 +0000
Subject: [PATCH] Fixup cholesky decomposition - rename to LLt and RtR to make
 convention clear

---
 cpp/Matrix.cpp     | 11 ++++++++---
 cpp/Matrix.h       |  5 ++++-
 cpp/testMatrix.cpp |  4 ++--
 3 files changed, 14 insertions(+), 6 deletions(-)

diff --git a/cpp/Matrix.cpp b/cpp/Matrix.cpp
index b1dffcf9a..562b8c385 100644
--- a/cpp/Matrix.cpp
+++ b/cpp/Matrix.cpp
@@ -787,7 +787,7 @@ double** createNRC(Matrix& A) {
 namespace BNU = boost::numeric::ublas;
 
 
-Matrix cholesky(const Matrix& A)
+Matrix LLt(const Matrix& A)
 {
 	assert(A.size1() == A.size2());
         Matrix L = zeros(A.size1(), A.size1());
@@ -808,12 +808,17 @@ Matrix cholesky(const Matrix& A)
         return L;
 }
 
+Matrix RtR(const Matrix &A)
+{
+	return trans(LLt(A));
+}
+
 /*
  * This is not ultra efficient, but not terrible, either.
  */
 Matrix cholesky_inverse(const Matrix &A)
 {
-        Matrix L = cholesky(A);
+        Matrix L = LLt(A);
         Matrix inv(boost::numeric::ublas::identity_matrix<double>(A.size1()));
         inplace_solve (L, inv, BNU::lower_tag ());
         return BNU::prod(trans(inv), inv);
@@ -883,7 +888,7 @@ Matrix inverse_square_root(const Matrix& A) {
 // inv(B) * inv(B)' == A
 // inv(B' * B) == 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()));
         inplace_solve (L, inv, BNU::lower_tag ());
 	return inv;
diff --git a/cpp/Matrix.h b/cpp/Matrix.h
index 53cefc220..bc75eb498 100644
--- a/cpp/Matrix.h
+++ b/cpp/Matrix.h
@@ -321,7 +321,10 @@ Matrix inverse_square_root(const Matrix& A);
 Matrix square_root_positive(const Matrix& A);
 
 /** 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. */
 Matrix cholesky_inverse(const Matrix &A);
diff --git a/cpp/testMatrix.cpp b/cpp/testMatrix.cpp
index c6b00094b..703282b4a 100644
--- a/cpp/testMatrix.cpp
+++ b/cpp/testMatrix.cpp
@@ -791,7 +791,7 @@ TEST( matrix, inverse_square_root )
 /* *********************************************************************** */
 // 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
-TEST( matrix, cholesky ) 
+TEST( matrix, LLt ) 
 {
 	Matrix M = Matrix_(5, 5,
 			   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.016485031422565, -0.012072546984405, -0.006621889326331, 0.014405837566082, 0.300462176944247);
 
-	EQUALITY(expected, cholesky(M));
+	EQUALITY(expected, LLt(M));
 }
 
 /* ************************************************************************* */