Merge pull request #108 from stan-dev/feature/issue-48-cholesky_decom…

…pose

Fixes #48. Feature/issue 48 cholesky decompose

stan/math/prim/mat/err/check_pos_definite.hpp

@@ -12,11 +12,11 @@
 #include <stan/math/prim/scal/err/check_positive_size.hpp>
 #include <stan/math/prim/mat/fun/Eigen.hpp>
 #include <stan/math/prim/mat/fun/value_of_rec.hpp>
-
 namespace stan {
 
   namespace math {
     using Eigen::Dynamic;
+
     /**
      * Return <code>true</code> if the specified square, symmetric
      * matrix is positive definite.
@@ -51,12 +51,63 @@ namespace stan {
         = value_of_rec(y).ldlt();
       if (cholesky.info() != Eigen::Success
           || !cholesky.isPositive()
-          || (cholesky.vectorD().array() <= CONSTRAINT_TOLERANCE).any())
+          || (cholesky.vectorD().array() <= 0.0).any())
         domain_error(function, name, y, "is not positive definite:\n");
       check_not_nan(function, name, y);
       return true;
     }
 
+    /**
+     * Return <code>true</code> if the specified LDLT transform of a matrix
+     * is positive definite.
+     *
+     * @tparam Derived Derived type of the Eigen::LDLT transform.
+     *
+     * @param function Function name (for error messages)
+     * @param name Variable name (for error messages)
+     * @param cholesky Eigen::LDLT to test, whose progenitor 
+     * must not have any NaN elements
+     *
+     * @return <code>true</code> if the matrix is positive definite
+     * @throw <code>std::domain_error</code> if the matrix is not positive definite.
+     */
+    template <typename Derived>
+    inline bool
+    check_pos_definite(const char* function,
+                       const char* name,
+                       const Eigen::LDLT<Derived>& cholesky) {
+      if (cholesky.info() != Eigen::Success
+          || !cholesky.isPositive()
+          || !(cholesky.vectorD().array() > 0.0).all())
+        domain_error(function, "LDLT decomposition of", " failed", name);
+      return true;
+    }
+
+    /**
+     * Return <code>true</code> if the specified LLT transform of a matrix
+     * is positive definite.
+     *
+     * @tparam Derived Derived type of the Eigen::LLT transform.
+     *
+     * @param function Function name (for error messages)
+     * @param name Variable name (for error messages)
+     * @param cholesky Eigen::LLT to test, whose progenitor 
+     * must not have any NaN elements
+     *
+     * @return <code>true</code> if the matrix is positive definite
+     * @throw <code>std::domain_error</code> if the diagonal of the L matrix is not positive.
+     */
+    template <typename Derived>
+    inline bool
+    check_pos_definite(const char* function,
+                       const char* name,
+                       const Eigen::LLT<Derived>& cholesky) {
+      if (cholesky.info() != Eigen::Success
+          || !(cholesky.matrixLLT().diagonal().array() > 0.0).all())
+        domain_error(function, "Cholesky decomposition of", " failed", name);
+      return true;
+    }
+
   }
 }
 #endif

stan/math/prim/mat/fun/cholesky_decompose.hpp

@@ -2,6 +2,7 @@
 #define STAN_MATH_PRIM_MAT_FUN_CHOLESKY_DECOMPOSE_HPP
 
 #include <stan/math/prim/mat/fun/Eigen.hpp>
+#include <stan/math/prim/mat/err/check_pos_definite.hpp>
 #include <stan/math/prim/mat/err/check_square.hpp>
 #include <stan/math/prim/mat/err/check_symmetric.hpp>
 
@@ -16,7 +17,8 @@ namespace stan {
      * <p>\f$A = L \times L^T\f$.
      * @param m Symmetrix matrix.
      * @return Square root of matrix.
-     * @throw std::domain_error if m is not a symmetric matrix.
+     * @throw std::domain_error if m is not a symmetric matrix or
+     *   if m is not positive definite (if m has more than 0 elements)
      */
     template <typename T>
     Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>
@@ -27,6 +29,7 @@ namespace stan {
       Eigen::LLT<Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> >
         llt(m.rows());
       llt.compute(m);
+      stan::math::check_pos_definite("cholesky_decompose", "m", llt);
       return llt.matrixL();
     }
 

stan/math/prim/mat/fun/cov_matrix_constrain.hpp

@@ -73,6 +73,7 @@ namespace stan {
                                                  Eigen::Dynamic> >::type K,
          T& lp) {
       using std::exp;
+      using std::log;
 
       using Eigen::Dynamic;
       using Eigen::Matrix;

test/unit/math/fwd/mat/fun/cholesky_decompose_test.cpp

@@ -2,17 +2,49 @@
 #include <stan/math/fwd/mat/fun/typedefs.hpp>
 #include <gtest/gtest.h>
 #include <stan/math/fwd/scal/fun/value_of.hpp>
+#include <stan/math/prim/mat/fun/value_of_rec.hpp>
+#include <stan/math/fwd/scal/fun/value_of_rec.hpp>
 #include <stan/math/fwd/core.hpp>
-#include <stan/math/fwd/scal/fun/sqrt.hpp>
 #include <stan/math/fwd/scal/fun/fabs.hpp>
+#include <stan/math/fwd/scal/fun/sqrt.hpp>
+#include <cmath>
+
+template <typename T>
+void deriv_chol_fwd(const Eigen::Matrix<T, -1, -1>& parent_mat,
+                    Eigen::Matrix<double,-1, -1>& vals,
+                    Eigen::Matrix<double, -1, -1>& gradients) {
+  using stan::math::value_of_rec;
+  Eigen::Matrix<double,2,2> parent_mat_d = value_of_rec(parent_mat);
+  vals(0,0) = sqrt(parent_mat_d(0,0));
+  vals(1,0) = parent_mat_d(1,0) / sqrt(parent_mat_d(0,0));
+  vals(0,1) = 0.0;
+  vals(1,1) = sqrt(parent_mat_d(1,1) - pow(vals(1,0), 2));
+
+  double pow_neg_half_00 = pow(parent_mat_d(0,0), -0.5);
+  double pow_neg_half_comb = pow(parent_mat_d(1,1) 
+                                 - pow(parent_mat_d(1,0),2.0) 
+                                 / parent_mat_d(0,0), -0.5);
+  gradients(0,0) = pow_neg_half_00 / 2 * value_of_rec(parent_mat(0,0).d_);
+  gradients(1,0) = -0.5 * value_of_rec(parent_mat(0,0).d_) 
+    * pow(parent_mat_d(0,0), -1.5)
+    + value_of_rec(parent_mat(1,0).d_) * pow_neg_half_00 
+    * parent_mat_d(1,0);
+  gradients(0,1) = 0.0;
+  gradients(1,1) = 0.5 * pow_neg_half_comb 
+    * (value_of_rec(parent_mat(1,1).d_) 
+       * - 2 * parent_mat_d(1,0) / parent_mat_d(0,0) 
+       * value_of_rec(parent_mat(1,0).d_) 
+       + pow(parent_mat_d(1,0),2.0) / pow(parent_mat_d(0,0),2.0) 
+       * value_of_rec(parent_mat(0,0).d_));
+}
 
 TEST(AgradFwdMatrixCholeskyDecompose, exception_mat_fd) {
   stan::math::matrix_fd m;
 
   m.resize(2,2);
   m << 1.0, 2.0, 
     2.0, 3.0;
-  EXPECT_NO_THROW(stan::math::cholesky_decompose(m));
+  EXPECT_THROW(stan::math::cholesky_decompose(m),std::domain_error);
 
   m.resize(0, 0);
   EXPECT_NO_THROW(stan::math::cholesky_decompose(m));
@@ -32,7 +64,7 @@ TEST(AgradFwdMatrixCholeskyDecompose, exception_mat_ffd) {
   m.resize(2,2);
   m << 1.0, 2.0, 
     2.0, 3.0;
-  EXPECT_NO_THROW(stan::math::cholesky_decompose(m));
+  EXPECT_THROW(stan::math::cholesky_decompose(m),std::domain_error);
 
   m.resize(0, 0);
   EXPECT_NO_THROW(stan::math::cholesky_decompose(m));
@@ -48,7 +80,7 @@ TEST(AgradFwdMatrixCholeskyDecompose, exception_mat_ffd) {
 }
 TEST(AgradFwdMatrixCholeskyDecompose, mat_fd) {
   stan::math::matrix_fd m0(2,2);
-  m0 << 1, 2, 2, 4;
+  m0 << 2, 1, 1, 2;
   m0(0,0).d_ = 1.0;
   m0(0,1).d_ = 1.0;
   m0(1,0).d_ = 1.0;
@@ -57,19 +89,20 @@ TEST(AgradFwdMatrixCholeskyDecompose, mat_fd) {
   using stan::math::cholesky_decompose;
 
   stan::math::matrix_fd res = cholesky_decompose(m0);
-
-  EXPECT_FLOAT_EQ(1, res(0,0).val_);
-  EXPECT_FLOAT_EQ(0, res(0,1).val_);
-  EXPECT_FLOAT_EQ(2, res(1,0).val_);
-  EXPECT_FLOAT_EQ(4, res(1,1).val_);
-  EXPECT_FLOAT_EQ(0.5, res(0,0).d_);
-  EXPECT_FLOAT_EQ(0, res(0,1).d_);
-  EXPECT_FLOAT_EQ(0, res(1,0).d_);
-  EXPECT_FLOAT_EQ(1, res(1,1).d_);
+  Eigen::Matrix<double,-1,-1> res_mat(2,2);
+  Eigen::Matrix<double,-1,-1> d_mat(2,2);
+  deriv_chol_fwd(m0, res_mat, d_mat);
+  for (int i = 0; i < 2; ++i)
+    for (int j = 0; j < i; ++j) {
+     EXPECT_FLOAT_EQ(res_mat(i,j), res(i,j).val_);
+     EXPECT_FLOAT_EQ(d_mat(i,j),res(i,j).d_) << "Row: " << i
+       << "Col: " << j;
+    }
 }
+
 TEST(AgradFwdMatrixCholeskyDecompose, mat_ffd) {
   stan::math::matrix_ffd m0(2,2);
-  m0 << 1, 2, 2, 4;
+  m0 << 4, 1, 1, 4;
   m0(0,0).d_ = 1.0;
   m0(0,1).d_ = 1.0;
   m0(1,0).d_ = 1.0;
@@ -78,13 +111,13 @@ TEST(AgradFwdMatrixCholeskyDecompose, mat_ffd) {
   using stan::math::cholesky_decompose;
 
   stan::math::matrix_ffd res = cholesky_decompose(m0);
-
-  EXPECT_FLOAT_EQ(1, res(0,0).val_.val_);
-  EXPECT_FLOAT_EQ(0, res(0,1).val_.val_);
-  EXPECT_FLOAT_EQ(2, res(1,0).val_.val_);
-  EXPECT_FLOAT_EQ(4, res(1,1).val_.val_);
-  EXPECT_FLOAT_EQ(0.5, res(0,0).d_.val_);
-  EXPECT_FLOAT_EQ(0, res(0,1).d_.val_);
-  EXPECT_FLOAT_EQ(0, res(1,0).d_.val_);
-  EXPECT_FLOAT_EQ(1, res(1,1).d_.val_);
+  Eigen::Matrix<double,-1,-1> res_mat(2,2);
+  Eigen::Matrix<double,-1,-1> d_mat(2,2);
+  deriv_chol_fwd(m0, res_mat, d_mat);
+  for (int i = 0; i < 2; ++i)
+    for (int j = 0; j < i; ++j) {
+     EXPECT_FLOAT_EQ(res_mat(i,j), res(i,j).val_.val_);
+     EXPECT_FLOAT_EQ(d_mat(i,j),res(i,j).d_.val_) << "Row: " << i
+       << " Col: " << j;
+    }
 }

test/unit/math/fwd/mat/fun/mdivide_left_spd_test.cpp

@@ -187,6 +187,18 @@ TEST(AgradFwdMatrixMdivideLeftSPD,fd_exceptions) {
   matrix_d fd1(3,3), fd2(4,4);
   row_vector_d rvd1(3), rvd2(4);
   vector_d vd1(3), vd2(4);
+  fv1.setZero();
+  fv2.setZero();
+  rvf1.setZero();
+  rvf2.setZero();
+  vf1.setZero();
+  vf2.setZero();
+  fd1.setZero();
+  fd2.setZero();
+  rvd1.setZero();
+  rvd2.setZero();
+  vd1.setZero();
+  vd2.setZero();
 
   EXPECT_THROW(mdivide_left_spd(fv1, fd2), std::domain_error);
   EXPECT_THROW(mdivide_left_spd(fd1, fv2), std::domain_error);
@@ -387,7 +399,18 @@ TEST(AgradFwdMatrixMdivideLeftSPD,ffd_exceptions) {
   matrix_d fd1(3,3), fd2(4,4);
   row_vector_d rvd1(3), rvd2(4);
   vector_d vd1(3), vd2(4);
-
+  fv1.setZero();
+  fv2.setZero();
+  rvf1.setZero();
+  rvf2.setZero();
+  vf1.setZero();
+  vf2.setZero();
+  fd1.setZero();
+  fd2.setZero();
+  rvd1.setZero();
+  rvd2.setZero();
+  vd1.setZero();
+  vd2.setZero();
   EXPECT_THROW(mdivide_left_spd(fv1, fd2), std::domain_error);
   EXPECT_THROW(mdivide_left_spd(fd1, fv2), std::domain_error);
   EXPECT_THROW(mdivide_left_spd(fv1, fv2), std::domain_error);

test/unit/math/fwd/mat/fun/mdivide_right_spd_test.cpp

@@ -190,7 +190,18 @@ TEST(AgradFwdMatrixMdivideRightSPD,fd_exceptions) {
   matrix_d fd1(3,3), fd2(4,4);
   row_vector_d rvd1(3), rvd2(4);
   row_vector_d vd1(3), vd2(4);
-
+  fv1.setZero();
+  fv2.setZero();
+  rvf1.setZero();
+  rvf2.setZero();
+  vf1.setZero();
+  vf2.setZero();
+  fd1.setZero();
+  fd2.setZero();
+  rvd1.setZero();
+  rvd2.setZero();
+  vd1.setZero();
+  vd2.setZero();
   EXPECT_THROW(mdivide_right_spd(fd2, fv1), std::invalid_argument);
   EXPECT_THROW(mdivide_right_spd(fv2, fd1), std::invalid_argument);
   EXPECT_THROW(mdivide_right_spd(fv2, fv1), std::invalid_argument);