From 0434a7bdd950534b0e05da46a220123398f2b192 Mon Sep 17 00:00:00 2001 From: Sean Pinkney Date: Wed, 6 Dec 2023 12:20:55 -0500 Subject: [PATCH 1/4] use less memory --- .../prim/prob/multi_normal_cholesky_lpdf.hpp | 62 ++++++++----------- 1 file changed, 27 insertions(+), 35 deletions(-) diff --git a/stan/math/prim/prob/multi_normal_cholesky_lpdf.hpp b/stan/math/prim/prob/multi_normal_cholesky_lpdf.hpp index bda51f702b6..ad216aae113 100644 --- a/stan/math/prim/prob/multi_normal_cholesky_lpdf.hpp +++ b/stan/math/prim/prob/multi_normal_cholesky_lpdf.hpp @@ -53,6 +53,9 @@ return_type_t multi_normal_cholesky_lpdf( using T_partials_return = partials_return_t; using matrix_partials_t = Eigen::Matrix; + using vector_partials_t = Eigen::Matrix; + using row_vector_partials_t + = Eigen::Matrix; using T_y_ref = ref_type_t; using T_mu_ref = ref_type_t; using T_L_ref = ref_type_t; @@ -119,59 +122,48 @@ return_type_t multi_normal_cholesky_lpdf( } if (include_summand::value) { - Eigen::Matrix - y_val_minus_mu_val(size_y, size_vec); + row_vector_partials_t half(size_vec); + vector_partials_t y_val_minus_mu_val(size_vec); + vector_partials_t scaled_diff(size_vec); + matrix_partials_t L_val = value_of(L_ref); + + T_partials_return sum_lp_vec(0.0); for (size_t i = 0; i < size_vec; i++) { decltype(auto) y_val = as_value_column_vector_or_scalar(y_vec[i]); decltype(auto) mu_val = as_value_column_vector_or_scalar(mu_vec[i]); - y_val_minus_mu_val.col(i) = y_val - mu_val; + y_val_minus_mu_val = eval(y_val - mu_val); + half = mdivide_left_tri(L_val, y_val_minus_mu_val).transpose(); + scaled_diff = mdivide_right_tri(half, L_val).transpose(); + + sum_lp_vec += dot_self(half); + + if (!is_constant_all::value) { + partials_vec<0>(ops_partials)[i] += -scaled_diff; + } + if (!is_constant_all::value) { + partials_vec<1>(ops_partials)[i] += scaled_diff; + } + if (!is_constant::value) { + partials_vec<2>(ops_partials)[i] += scaled_diff * half;; + } } - matrix_partials_t half; - matrix_partials_t scaled_diff; + logp += -0.5 * sum_lp_vec; // If the covariance is not autodiff, we can avoid computing a matrix // inverse if (is_constant::value) { - matrix_partials_t L_val = value_of(L_ref); - - half = mdivide_left_tri(L_val, y_val_minus_mu_val) - .transpose(); - - scaled_diff = mdivide_right_tri(half, L_val).transpose(); - if (include_summand::value) { logp -= sum(log(L_val.diagonal())) * size_vec; } } else { - matrix_partials_t inv_L_val + matrix_partials_t inv_L_val = mdivide_left_tri(value_of(L_ref)); - half = (inv_L_val.template triangularView() - * y_val_minus_mu_val) - .transpose(); - - scaled_diff = (half * inv_L_val.template triangularView()) - .transpose(); - logp += sum(log(inv_L_val.diagonal())) * size_vec; - partials<2>(ops_partials) -= size_vec * inv_L_val.transpose(); - for (size_t i = 0; i < size_vec; i++) { - partials_vec<2>(ops_partials)[i] += scaled_diff.col(i) * half.row(i); - } - } - - logp -= 0.5 * sum(columns_dot_self(half)); - - for (size_t i = 0; i < size_vec; i++) { - if (!is_constant_all::value) { - partials_vec<0>(ops_partials)[i] -= scaled_diff.col(i); - } - if (!is_constant_all::value) { - partials_vec<1>(ops_partials)[i] += scaled_diff.col(i); - } + partials<2>(ops_partials) -= size_vec * inv_L_val.transpose(); } } From ccf57838f808fd1905323a38d2063074258e16e3 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Wed, 6 Dec 2023 12:25:33 -0500 Subject: [PATCH 2/4] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/prim/prob/multi_normal_cholesky_lpdf.hpp | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/stan/math/prim/prob/multi_normal_cholesky_lpdf.hpp b/stan/math/prim/prob/multi_normal_cholesky_lpdf.hpp index ad216aae113..f0332d6e30b 100644 --- a/stan/math/prim/prob/multi_normal_cholesky_lpdf.hpp +++ b/stan/math/prim/prob/multi_normal_cholesky_lpdf.hpp @@ -54,7 +54,7 @@ return_type_t multi_normal_cholesky_lpdf( using matrix_partials_t = Eigen::Matrix; using vector_partials_t = Eigen::Matrix; - using row_vector_partials_t + using row_vector_partials_t = Eigen::Matrix; using T_y_ref = ref_type_t; using T_mu_ref = ref_type_t; @@ -132,9 +132,10 @@ return_type_t multi_normal_cholesky_lpdf( for (size_t i = 0; i < size_vec; i++) { decltype(auto) y_val = as_value_column_vector_or_scalar(y_vec[i]); decltype(auto) mu_val = as_value_column_vector_or_scalar(mu_vec[i]); - y_val_minus_mu_val = eval(y_val - mu_val); - half = mdivide_left_tri(L_val, y_val_minus_mu_val).transpose(); - scaled_diff = mdivide_right_tri(half, L_val).transpose(); + y_val_minus_mu_val = eval(y_val - mu_val); + half = mdivide_left_tri(L_val, y_val_minus_mu_val) + .transpose(); + scaled_diff = mdivide_right_tri(half, L_val).transpose(); sum_lp_vec += dot_self(half); @@ -145,7 +146,8 @@ return_type_t multi_normal_cholesky_lpdf( partials_vec<1>(ops_partials)[i] += scaled_diff; } if (!is_constant::value) { - partials_vec<2>(ops_partials)[i] += scaled_diff * half;; + partials_vec<2>(ops_partials)[i] += scaled_diff * half; + ; } } @@ -158,7 +160,7 @@ return_type_t multi_normal_cholesky_lpdf( logp -= sum(log(L_val.diagonal())) * size_vec; } } else { - matrix_partials_t inv_L_val + matrix_partials_t inv_L_val = mdivide_left_tri(value_of(L_ref)); logp += sum(log(inv_L_val.diagonal())) * size_vec; From 1284e64fd8d84d890526d13612f4041e60a17c08 Mon Sep 17 00:00:00 2001 From: Sean Pinkney Date: Wed, 6 Dec 2023 12:30:22 -0500 Subject: [PATCH 3/4] fix double semi-colon --- stan/math/prim/prob/multi_normal_cholesky_lpdf.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stan/math/prim/prob/multi_normal_cholesky_lpdf.hpp b/stan/math/prim/prob/multi_normal_cholesky_lpdf.hpp index ad216aae113..c83879583b9 100644 --- a/stan/math/prim/prob/multi_normal_cholesky_lpdf.hpp +++ b/stan/math/prim/prob/multi_normal_cholesky_lpdf.hpp @@ -145,7 +145,7 @@ return_type_t multi_normal_cholesky_lpdf( partials_vec<1>(ops_partials)[i] += scaled_diff; } if (!is_constant::value) { - partials_vec<2>(ops_partials)[i] += scaled_diff * half;; + partials_vec<2>(ops_partials)[i] += scaled_diff * half; } } From a08db5ea10a9d1d2cf8d8f7b0a1f63fdbf375c91 Mon Sep 17 00:00:00 2001 From: Sean Pinkney Date: Tue, 19 Dec 2023 22:22:38 -0500 Subject: [PATCH 4/4] Update stan/math/prim/prob/multi_normal_cholesky_lpdf.hpp Co-authored-by: Steve Bronder --- stan/math/prim/prob/multi_normal_cholesky_lpdf.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stan/math/prim/prob/multi_normal_cholesky_lpdf.hpp b/stan/math/prim/prob/multi_normal_cholesky_lpdf.hpp index b9e98ffa60b..47fecc0682e 100644 --- a/stan/math/prim/prob/multi_normal_cholesky_lpdf.hpp +++ b/stan/math/prim/prob/multi_normal_cholesky_lpdf.hpp @@ -132,7 +132,7 @@ return_type_t multi_normal_cholesky_lpdf( for (size_t i = 0; i < size_vec; i++) { decltype(auto) y_val = as_value_column_vector_or_scalar(y_vec[i]); decltype(auto) mu_val = as_value_column_vector_or_scalar(mu_vec[i]); - y_val_minus_mu_val = eval(y_val - mu_val); + y_val_minus_mu_val = y_val - mu_val; half = mdivide_left_tri(L_val, y_val_minus_mu_val) .transpose(); scaled_diff = mdivide_right_tri(half, L_val).transpose();