From 2f32e485174c2049a7c364745b2aad091449602b Mon Sep 17 00:00:00 2001
From: Gael Guennebaud <g.gael@free.fr>
Date: Mon, 20 Jun 2011 15:05:50 +0200
Subject: [PATCH] New feature: add rank one update in Cholesky decomposition

---
 Eigen/src/Cholesky/LLT.h | 52 ++++++++++++++++++++++++++++++++++++++++
 test/cholesky.cpp        | 23 +++++++++++++++++-
 2 files changed, 74 insertions(+), 1 deletion(-)

diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h
index a8fc525e8..33a5c4165 100644
--- a/Eigen/src/Cholesky/LLT.h
+++ b/Eigen/src/Cholesky/LLT.h
@@ -178,6 +178,9 @@ template<typename _MatrixType, int _UpLo> class LLT
     inline Index rows() const { return m_matrix.rows(); }
     inline Index cols() const { return m_matrix.cols(); }
 
+    template<typename VectorType>
+    void rankUpdate(const VectorType& vec);
+
   protected:
     /** \internal
       * Used to compute and store L
@@ -254,6 +257,34 @@ template<> struct llt_inplace<Lower>
     }
     return -1;
   }
+
+  template<typename MatrixType, typename VectorType>
+  static void rankUpdate(MatrixType& mat, const VectorType& vec)
+  {
+    typedef typename MatrixType::ColXpr ColXpr;
+    typedef typename internal::remove_all<ColXpr>::type ColXprCleaned;
+    typedef typename ColXprCleaned::SegmentReturnType ColXprSegment;
+    typedef typename VectorType::SegmentReturnType VecSegment;
+
+    int n = mat.cols();
+    eigen_assert(mat.rows()==n && vec.size()==n);
+    typedef typename MatrixType::Scalar Scalar;
+    Matrix<Scalar,Dynamic,1> temp(vec);
+
+    for(int i=0; i<n; ++i)
+    {
+      JacobiRotation<Scalar> g;
+      g.makeGivens(mat(i,i), -temp(i), &mat(i,i));
+
+      int rs = n-i-1;
+      if(rs>0)
+      {
+        ColXprSegment x(mat.col(i).tail(rs));
+        VecSegment y(temp.tail(rs));
+        apply_rotation_in_the_plane(x, y, g);
+      }
+    }
+  }
 };
 
 template<> struct llt_inplace<Upper>
@@ -270,6 +301,12 @@ template<> struct llt_inplace<Upper>
     Transpose<MatrixType> matt(mat);
     return llt_inplace<Lower>::blocked(matt);
   }
+  template<typename MatrixType, typename VectorType>
+  static void rankUpdate(MatrixType& mat, const VectorType& vec)
+  {
+    Transpose<MatrixType> matt(mat);
+    return llt_inplace<Lower>::rankUpdate(matt, vec);
+  }
 };
 
 template<typename MatrixType> struct LLT_Traits<MatrixType,Lower>
@@ -314,6 +351,20 @@ LLT<MatrixType,_UpLo>& LLT<MatrixType,_UpLo>::compute(const MatrixType& a)
   return *this;
 }
 
+/** Performs a rank one update of the current decomposition.
+  * If A = LL^* before the rank one update,
+  * then after it we have LL^* = A + vv^* where \a v must be a vector
+  * of same dimension.
+  *
+  */
+template<typename MatrixType, int _UpLo>
+template<typename VectorType>
+void LLT<MatrixType,_UpLo>::rankUpdate(const VectorType& v)
+{
+  EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorType);
+  internal::llt_inplace<UpLo>::rankUpdate(m_matrix,v);
+}
+    
 namespace internal {
 template<typename _MatrixType, int UpLo, typename Rhs>
 struct solve_retval<LLT<_MatrixType, UpLo>, Rhs>
@@ -384,3 +435,4 @@ SelfAdjointView<MatrixType, UpLo>::llt() const
 }
 
 #endif // EIGEN_LLT_H
+
diff --git a/test/cholesky.cpp b/test/cholesky.cpp
index 293a375ad..35fa5953d 100644
--- a/test/cholesky.cpp
+++ b/test/cholesky.cpp
@@ -166,6 +166,10 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
       VERIFY_EVALUATION_COUNT(matX = ldltup.solve(matX), 0);
       VERIFY_IS_APPROX(matX, ldltup.solve(matB).eval());
     }
+
+    // restore
+    if(sign == -1)
+      symm = -symm;
   }
 
   // test some special use cases of SelfCwiseBinaryOp:
@@ -182,6 +186,24 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
   m2 = m1;
   m2.noalias() -= symmLo.template selfadjointView<Lower>().llt().solve(matB);
   VERIFY_IS_APPROX(m2, m1 - symmLo.template selfadjointView<Lower>().llt().solve(matB));
+
+  // Cholesky update/downdate
+  {
+    MatrixType symmLo = symm.template triangularView<Lower>();
+    MatrixType symmUp = symm.template triangularView<Upper>();
+    
+    VectorType vec = VectorType::Random(rows);
+
+    MatrixType symmCpy = symm + vec * vec.adjoint();
+
+    LLT<MatrixType,Lower> chollo(symmLo);
+    chollo.rankUpdate(vec);
+    VERIFY_IS_APPROX(symmCpy,  chollo.reconstructedMatrix());
+
+    LLT<MatrixType,Upper> cholup(symmUp);
+    cholup.rankUpdate(vec);
+    VERIFY_IS_APPROX(symmCpy,  cholup.reconstructedMatrix());
+  }
   
 }
 
@@ -242,7 +264,6 @@ template<typename MatrixType> void cholesky_cplx(const MatrixType& m)
 //     matX = ldltlo.solve(matB);
 //     VERIFY_IS_APPROX(symm * matX, matB);
   }
-
 }
 
 template<typename MatrixType> void cholesky_verify_assert()