@@ -139,8 +139,8 @@ inline void block_matrix_sqrt(WRootMat& W_root,
139139
140140template <typename WRootMat>
141141inline void block_matrix_chol_L (WRootMat& W_root,
142- const Eigen::SparseMatrix<double >& W,
143- const Eigen::Index block_size) {
142+ const Eigen::SparseMatrix<double >& W,
143+ const Eigen::Index block_size) {
144144 int n_block = W.cols () / block_size;
145145 Eigen::MatrixXd local_block (block_size, block_size);
146146 Eigen::MatrixXd local_block_sqrt (block_size, block_size);
@@ -161,13 +161,13 @@ inline void block_matrix_chol_L(WRootMat& W_root,
161161 // Compute square root of T
162162 Eigen::LLT<Eigen::Ref<Eigen::MatrixXd>> llt (local_block);
163163 if (llt.info () != Eigen::Success) {
164- throw std::runtime_error (" Cholesky failed on block " + std::to_string (i));
164+ throw std::runtime_error (" Cholesky failed on block "
165+ + std::to_string (i));
165166 }
166167 const auto Lb = llt.matrixL ();
167168 for (int k = 0 ; k < block_size; k++) {
168169 for (int j = k; j < block_size; j++) {
169- W_root.coeffRef (i * block_size + j, i * block_size + k)
170- = Lb (j, k);
170+ W_root.coeffRef (i * block_size + j, i * block_size + k) = Lb (j, k);
171171 }
172172 }
173173 } catch (const std::exception& e) {
@@ -210,7 +210,6 @@ inline void block_matrix_chol_L(WRootMat& W_root,
210210 }
211211}
212212
213-
214213/* *
215214 * @brief Performs a simple line search
216215 *
@@ -494,16 +493,15 @@ inline auto laplace_marginal_density_est(LLFun&& ll_fun, LLTupleArgs&& ll_args,
494493 }
495494 }
496495 B.noalias () = MatrixXd::Identity (theta_size, theta_size)
497- + W_r.asDiagonal () * covariance
498- * W_r.asDiagonal ();
496+ + W_r.asDiagonal () * covariance * W_r.asDiagonal ();
499497 Eigen::LLT<Eigen::Ref<Eigen::MatrixXd>> llt_B (B);
500498 auto L = llt_B.matrixL ();
501499 auto LT = llt_B.matrixU ();
502500 b.noalias () = W.diagonal ().cwiseProduct (theta) + theta_grad;
503- a.noalias () = b
504- - W_r. asDiagonal ()
505- * LT. solve (L. solve (
506- W_r.cwiseProduct (covariance * b)));
501+ a.noalias ()
502+ = b
503+ - W_r. asDiagonal ()
504+ * LT. solve (L. solve ( W_r.cwiseProduct (covariance * b)));
507505 // Simple Newton step
508506 theta.noalias () = covariance * a;
509507 objective_old = objective_new;
@@ -1009,21 +1007,38 @@ inline auto laplace_marginal_density(const LLFun& ll_fun, LLTupleArgs&& ll_args,
10091007 partial_parm, ll_args_filter);
10101008 if (options.solver == 1 ) {
10111009 if (options.hessian_block_size == 1 ) {
1012- // TODO(Steve): Solve without casting from sparse to dense
1013- Eigen::MatrixXd tmp
1014- = md_est.L .template triangularView <Eigen::Lower>().solve (
1015- md_est.W_r .toDense ());
1016- R = tmp.transpose () * tmp;
1017- arena_t <Eigen::MatrixXd> C
1018- = md_est.L .template triangularView <Eigen::Lower>().solve (
1019- md_est.W_r * md_est.covariance );
1020- if constexpr (!ll_args_contain_var) {
1021- s2.deep_copy (
1022- (0.5
1023- * (md_est.covariance .diagonal () - (C.transpose () * C).diagonal ())
1024- .cwiseProduct (laplace_likelihood::third_diff (
1025- ll_fun, md_est.theta , value_of (ll_args_copy), msgs))));
1010+ // TODO(Steve): Solve without casting from sparse to dense
1011+ Eigen::MatrixXd tmp
1012+ = md_est.L .template triangularView <Eigen::Lower>().solve (
1013+ md_est.W_r .toDense ());
1014+ R = tmp.transpose () * tmp;
1015+ arena_t <Eigen::MatrixXd> C
1016+ = md_est.L .template triangularView <Eigen::Lower>().solve (
1017+ md_est.W_r * md_est.covariance );
1018+ if constexpr (!ll_args_contain_var) {
1019+ s2.deep_copy (
1020+ (0.5
1021+ * (md_est.covariance .diagonal () - (C.transpose () * C).diagonal ())
1022+ .cwiseProduct (laplace_likelihood::third_diff (
1023+ ll_fun, md_est.theta , value_of (ll_args_copy), msgs))));
1024+ } else {
1025+ arena_t <Eigen::MatrixXd> A = md_est.covariance - C.transpose () * C;
1026+ auto s2_tmp = laplace_likelihood::compute_s2 (
1027+ ll_fun, md_est.theta , A, options.hessian_block_size , ll_args_copy,
1028+ msgs);
1029+ s2.deep_copy (s2_tmp);
1030+ copy_compute_s2 (partial_parm, ll_args_filter);
1031+ set_zero_adjoint (ll_args_filter);
1032+ }
1033+
10261034 } else {
1035+ Eigen::MatrixXd tmp
1036+ = md_est.L .template triangularView <Eigen::Lower>().solve (
1037+ md_est.W_r .toDense ());
1038+ R = tmp.transpose () * tmp;
1039+ arena_t <Eigen::MatrixXd> C
1040+ = md_est.L .template triangularView <Eigen::Lower>().solve (
1041+ md_est.W_r * md_est.covariance );
10271042 arena_t <Eigen::MatrixXd> A = md_est.covariance - C.transpose () * C;
10281043 auto s2_tmp = laplace_likelihood::compute_s2 (ll_fun, md_est.theta , A,
10291044 options.hessian_block_size ,
@@ -1032,23 +1047,6 @@ inline auto laplace_marginal_density(const LLFun& ll_fun, LLTupleArgs&& ll_args,
10321047 copy_compute_s2 (partial_parm, ll_args_filter);
10331048 set_zero_adjoint (ll_args_filter);
10341049 }
1035-
1036- } else {
1037- Eigen::MatrixXd tmp
1038- = md_est.L .template triangularView <Eigen::Lower>().solve (
1039- md_est.W_r .toDense ());
1040- R = tmp.transpose () * tmp;
1041- arena_t <Eigen::MatrixXd> C
1042- = md_est.L .template triangularView <Eigen::Lower>().solve (
1043- md_est.W_r * md_est.covariance );
1044- arena_t <Eigen::MatrixXd> A = md_est.covariance - C.transpose () * C;
1045- auto s2_tmp = laplace_likelihood::compute_s2 (ll_fun, md_est.theta , A,
1046- options.hessian_block_size ,
1047- ll_args_copy, msgs);
1048- s2.deep_copy (s2_tmp);
1049- copy_compute_s2 (partial_parm, ll_args_filter);
1050- set_zero_adjoint (ll_args_filter);
1051- }
10521050 } else if (options.solver == 2 ) {
10531051 R = md_est.W_r
10541052 - md_est.W_r * md_est.K_root
0 commit comments