adding checks for cholesky factors; we found problems with user-facin…

…g functions
stan-dev · syclik · Jan 13, 2024 · Jan 17, 2024 · Jan 17, 2024 · Jan 19, 2024
commit 66fff6d2b6f2d55eff11d9633dc1a4a9cd50ca60
diff --git a/stan/math/prim/prob/inv_wishart_cholesky_lpdf.hpp b/stan/math/prim/prob/inv_wishart_cholesky_lpdf.hpp
@@ -49,24 +49,16 @@ return_type_t<T_y, T_dof, T_scale> inv_wishart_cholesky_lpdf(
   using T_return = return_type_t<T_y, T_dof, T_scale>;
   static const char* function = "inv_wishart_cholesky_lpdf";
   Eigen::Index k = L_Y.rows();
-  check_size_match(function, "Rows of Random variable", L_Y.rows(),
-                   "columns of scale parameter", L_S.rows());
+  check_cholesky_factor(function, "Cholesky random variable", L_Y);
+  check_greater(function, "Degrees of freedom parameter", nu, k - 1);
+  check_cholesky_factor(function, "Cholesky Scale matrix", L_S);
   check_size_match(function, "Rows of random variable", L_Y.rows(),
-                   "columns of Random variable", L_Y.cols());
-  check_size_match(function, "Rows of scale parameter", L_S.rows(),
-                   "columns of scale parameter", L_S.cols());
+                   "columns of scale parameter", L_S.rows());
 
   T_L_Y_ref L_Y_ref = L_Y;
   T_nu_ref nu_ref = nu;
   T_L_S_ref L_S_ref = L_S;
 
-  check_greater(function, "Degrees of freedom parameter", nu_ref, k - 1);
-  check_positive(function, "Cholesky Random variable", L_Y_ref.diagonal());
-  check_positive(function, "columns of Cholesky Random variable",
-                 L_Y_ref.cols());
-  check_positive(function, "Cholesky Scale matrix", L_S_ref.diagonal());
-  check_positive(function, "columns of Cholesky Scale matrix", L_S_ref.cols());
-
   T_return lp(0.0);
 
   if (include_summand<propto, T_dof>::value) {

diff --git a/stan/math/prim/prob/inv_wishart_cholesky_rng.hpp b/stan/math/prim/prob/inv_wishart_cholesky_rng.hpp
@@ -3,7 +3,7 @@
 
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err.hpp>
-#include <stan/math/prim/fun/mdivide_left_tri_low.hpp>
+#include <stan/math/prim/fun/mdivide_right_tri_low.hpp>
 
 namespace stan {
 namespace math {
@@ -34,10 +34,8 @@ inline Eigen::MatrixXd inv_wishart_cholesky_rng(double nu,
   using Eigen::MatrixXd;
   static const char* function = "inv_wishart_cholesky_rng";
   index_type_t<MatrixXd> k = L_S.rows();
-  check_square(function, "Cholesky Scale matrix", L_S);
   check_greater(function, "degrees of freedom > dims - 1", nu, k - 1);
-  check_positive(function, "Cholesky Scale matrix", L_S.diagonal());
-  check_positive(function, "columns of Cholesky Scale matrix", L_S.cols());
+  check_cholesky_factor(function, "Cholesky Scale matrix", L_S);
 
   MatrixXd B = MatrixXd::Zero(k, k);
   for (int j = 0; j < k; ++j) {
@@ -46,8 +44,8 @@ inline Eigen::MatrixXd inv_wishart_cholesky_rng(double nu,
     }
     B(j, j) = std::sqrt(chi_square_rng(nu - k + j + 1, rng));
   }
-
-  return mdivide_left_tri_low(B, L_S);
+  
+  return mdivide_right_tri_low(L_S, B);
 }
 
 }  // namespace math

diff --git a/stan/math/prim/prob/lkj_corr_cholesky_lpdf.hpp b/stan/math/prim/prob/lkj_corr_cholesky_lpdf.hpp
@@ -19,10 +19,10 @@ return_type_t<T_covar, T_shape> lkj_corr_cholesky_lpdf(const T_covar& L,
                                                        const T_shape& eta) {
   using lp_ret = return_type_t<T_covar, T_shape>;
   static const char* function = "lkj_corr_cholesky_lpdf";
-  const auto& L_ref = to_ref(L);
   check_positive(function, "Shape parameter", eta);
-  check_lower_triangular(function, "Random variable", L_ref);
+  check_cholesky_factor(function, "Random variable", L);
 
+  const auto& L_ref = to_ref(L);
   const unsigned int K = L.rows();
   if (K == 0) {
     return 0.0;

diff --git a/stan/math/prim/prob/multi_gp_cholesky_lpdf.hpp b/stan/math/prim/prob/multi_gp_cholesky_lpdf.hpp
@@ -47,11 +47,13 @@ return_type_t<T_y, T_covar, T_w> multi_gp_cholesky_lpdf(const T_y& y,
                    "Size of kernel scales (w)", w.size());
   check_size_match(function, "Size of random variable", y.cols(),
                    "rows of covariance parameter", L.rows());
+  check_positive_finite(function, "Kernel scales", w);
+  check_finite(function, "Random variable", y);
+  check_cholesky_factor(function, "Cholesky decomposition of kernel matrix", L);
+
   const auto& y_ref = to_ref(y);
   const auto& L_ref = to_ref(L);
   const auto& w_ref = to_ref(w);
-  check_positive_finite(function, "Kernel scales", w_ref);
-  check_finite(function, "Random variable", y_ref);
 
   if (y.rows() == 0) {
     return 0;

diff --git a/stan/math/prim/prob/multi_normal_cholesky_lpdf.hpp b/stan/math/prim/prob/multi_normal_cholesky_lpdf.hpp
@@ -99,13 +99,13 @@ return_type_t<T_y, T_loc, T_covar> multi_normal_cholesky_lpdf(
                    "size of location parameter", size_mu);
   check_size_match(function, "Size of random variable", size_y,
                    "rows of covariance parameter", L.rows());
-  check_size_match(function, "Size of random variable", size_y,
-                   "columns of covariance parameter", L.cols());
 
   for (size_t i = 0; i < size_vec; i++) {
     check_finite(function, "Location parameter", mu_vec[i]);
     check_not_nan(function, "Random variable", y_vec[i]);
   }
+  check_cholesky_factor(function, "Cholesky decomposition of a variance matrix",
+                        L);
 
   if (unlikely(size_y == 0)) {
     return T_return(0);

diff --git a/stan/math/prim/prob/multi_normal_cholesky_rng.hpp b/stan/math/prim/prob/multi_normal_cholesky_rng.hpp
@@ -55,6 +55,7 @@ multi_normal_cholesky_rng(
   for (size_t i = 0; i < N; i++) {
     check_finite(function, "Location parameter", mu_vec[i]);
   }
+  check_cholesky_factor(function, "Cholesky factor of covariance matrix", L);
 
   const auto& L_ref = to_ref(L);
 

diff --git a/stan/math/prim/prob/wishart_cholesky_lpdf.hpp b/stan/math/prim/prob/wishart_cholesky_lpdf.hpp
@@ -49,23 +49,17 @@ return_type_t<T_y, T_dof, T_scale> wishart_cholesky_lpdf(const T_y& L_Y,
   using T_return = return_type_t<T_y, T_dof, T_scale>;
   static const char* function = "wishart_cholesky_lpdf";
   Eigen::Index k = L_Y.rows();
-  check_size_match(function, "Rows of RSandom variable", L_Y.rows(),
-                   "columns of scale parameter", L_S.rows());
+
+  check_cholesky_factor(function, "Cholesky random variable", L_Y);
+  check_greater(function, "Degrees of freedom parameter", nu, k - 1);
+  check_cholesky_factor(function, "Cholesky scale matrix", L_S);
   check_size_match(function, "Rows of random variable", L_Y.rows(),
-                   "columns of random variable", L_Y.cols());
-  check_size_match(function, "Rows of scale parameter", L_S.rows(),
-                   "columns of scale parameter", L_S.cols());
+                   "columns of scale parameter", L_S.rows());
+
   T_L_Y_ref L_Y_ref = L_Y;
   T_nu_ref nu_ref = nu;
   T_L_S_ref L_S_ref = L_S;
 
-  check_greater(function, "Degrees of freedom parameter", nu_ref, k - 1);
-  check_positive(function, "Cholesky Random variable", L_Y_ref.diagonal());
-  check_positive(function, "columns of Cholesky Random variable",
-                 L_Y_ref.cols());
-  check_positive(function, "Cholesky scale matrix", L_S_ref.diagonal());
-  check_positive(function, "columns of Cholesky scale matrix", L_S_ref.cols());
-
   T_return lp(0.0);
 
   if (include_summand<propto, T_dof>::value) {

diff --git a/stan/math/prim/prob/wishart_cholesky_rng.hpp b/stan/math/prim/prob/wishart_cholesky_rng.hpp
@@ -34,10 +34,8 @@ inline Eigen::MatrixXd wishart_cholesky_rng(double nu,
   using Eigen::MatrixXd;
   static const char* function = "wishart_cholesky_rng";
   index_type_t<MatrixXd> k = L_S.rows();
-  check_square(function, "Cholesky Scale matrix", L_S);
   check_greater(function, "degrees of freedom > dims - 1", nu, k - 1);
-  check_positive(function, "Cholesky Scale matrix", L_S.diagonal());
-  check_positive(function, "columns of Cholesky Scale matrix", L_S.cols());
+  check_cholesky_factor(function, "Cholesky Scale matrix", L_S);
 
   MatrixXd B = MatrixXd::Zero(k, k);
   for (int j = 0; j < k; ++j) {

diff --git a/test/unit/math/prim/prob/inv_wishart_cholesky_rng_test.cpp b/test/unit/math/prim/prob/inv_wishart_cholesky_rng_test.cpp
@@ -17,7 +17,7 @@ TEST(ProbDistributionsInvWishartCholesky, rng) {
 
   MatrixXd omega(3, 4);
   EXPECT_THROW(inv_wishart_cholesky_rng(3.0, omega, rng),
-               std::invalid_argument);
+               std::domain_error);
 
   MatrixXd sigma(3, 3);
   sigma << 9.0, -3.0, 0.0, -3.0, 4.0, 1.0, 0.0, 1.0, 3.0;
@@ -124,19 +124,22 @@ TEST(ProbDistributionsInvWishartCholesky, compareToInvWishart) {
   using stan::math::qr_thin_Q;
 
   boost::random::mt19937 rng(92343U);
-  int N = 1e5;
+  int N = 1e4;
   double tol = 0.05;
-  for (int k = 1; k < 5; k++) {
+  for (int k = 1; k < 4; k++) {
     MatrixXd L = qr_thin_Q(MatrixXd::Random(k, k)).transpose();
     L.diagonal() = stan::math::abs(L.diagonal());
     MatrixXd sigma = multiply_lower_tri_self_transpose(L);
-    MatrixXd Z_mean = sigma / (k + 3);
+    MatrixXd L_S = stan::math::cholesky_decompose(sigma);
+    MatrixXd Z_mean = MatrixXd::Zero(k, k);
     MatrixXd Z_est = MatrixXd::Zero(k, k);
     for (int i = 0; i < N; i++) {
       Z_est += multiply_lower_tri_self_transpose(
-          inv_wishart_cholesky_rng(k + 4, L, rng));
+          inv_wishart_cholesky_rng(k + 4, L_S, rng));
+      Z_mean += inv_wishart_rng(k + 4, sigma, rng);
     }
     Z_est /= N;
+    Z_mean /= N;
     for (int j = 0; j < k; j++) {
       for (int i = 0; i < j; i++) {
         EXPECT_NEAR(Z_est(i, j), Z_mean(i, j), tol);

diff --git a/test/unit/math/prim/prob/inv_wishart_cholesky_test.cpp b/test/unit/math/prim/prob/inv_wishart_cholesky_test.cpp
@@ -163,12 +163,12 @@ TEST(ProbDistributionsInvWishartCholesky, Error) {
   nu = 5;
   MatrixXd Sigma(2, 1);
   EXPECT_THROW(inv_wishart_cholesky_lpdf(MatrixXd::Identity(1, 1), nu, Sigma),
-               std::invalid_argument);
+               std::domain_error);
 
   nu = 5;
   MatrixXd Y(2, 1);
   EXPECT_THROW(inv_wishart_cholesky_lpdf(Y, nu, MatrixXd::Identity(2, 2)),
-               std::invalid_argument);
+               std::domain_error);
 
   nu = 5;
   EXPECT_THROW(inv_wishart_cholesky_lpdf(MatrixXd::Identity(3, 3), nu,

diff --git a/test/unit/math/prim/prob/wishart_cholesky_test.cpp b/test/unit/math/prim/prob/wishart_cholesky_test.cpp
@@ -211,12 +211,12 @@ TEST(ProbDistributionsWishartCholesky, error) {
   nu = 5;
   MatrixXd Sigma(2, 1);
   EXPECT_THROW(wishart_cholesky_lpdf(MatrixXd::Identity(1, 1), nu, Sigma),
-               std::invalid_argument);
+               std::domain_error);
 
   nu = 5;
   MatrixXd Y(2, 1);
   EXPECT_THROW(wishart_cholesky_lpdf(Y, nu, MatrixXd::Identity(2, 2)),
-               std::invalid_argument);
+               std::domain_error);
 
   nu = 5;
   EXPECT_THROW(wishart_cholesky_lpdf(MatrixXd::Identity(3, 3), nu,