diff --git a/src/math/cholesky_helpers.h b/src/math/cholesky_helpers.h index 4bbf09e2..24142c3a 100644 --- a/src/math/cholesky_helpers.h +++ b/src/math/cholesky_helpers.h @@ -10,7 +10,9 @@ */ #include +#include #include +#include "math/explog_macros.h" namespace cholesky_helpers { @@ -42,4 +44,90 @@ inline double compute_inv_submatrix_i(const arma::mat& A, size_t i, return A(ii, jj) - A(ii, i) * A(jj, i) / A(i, i); } +/** + * Rank-2 matrix-determinant-lemma log-ratio log|K'| − log|K| for a precision + * change confined to entries (i,j), (j,i), (j,j). cc11/cc12/cc22 are the 2×2 + * update-Gram-matrix entries I + Vᵀ Σ U. + * + * Shared by GGMModel::log_det_ratio_edge and MixedMRFModel::log_det_ratio_yy_edge, + * which differ only in where the current precision entries come from (GGM stores + * K directly; Mixed stores −½K), so they pass the resolved precision scalars in. + * + * @param covariance Σ = K⁻¹. + * @param i, j Off-diagonal/diagonal indices (i < j). + * @param precision_curr_ij, precision_curr_jj Current K entries. + * @param precision_prop_ij, precision_prop_jj Proposed K entries. + */ +inline double log_det_ratio_edge_kernel( + const arma::mat& covariance, size_t i, size_t j, + double precision_curr_ij, double precision_prop_ij, + double precision_curr_jj, double precision_prop_jj) { + double Ui2 = precision_curr_ij - precision_prop_ij; + double Uj2 = (precision_curr_jj - precision_prop_jj) / 2; + + double cc11 = covariance(j, j); + double cc12 = 1.0 - (covariance(i, j) * Ui2 + covariance(j, j) * Uj2); + double cc22 = Ui2 * Ui2 * covariance(i, i) + + 2.0 * Ui2 * Uj2 * covariance(i, j) + + Uj2 * Uj2 * covariance(j, j); + + return MY_LOG(std::abs(cc11 * cc22 - cc12 * cc12)); +} + +/** + * Rank-1 specialisation of log_det_ratio_edge_kernel (Ui2 = 0): a diagonal-only + * precision change at (j,j). Shared by GGMModel::log_det_ratio_diag and + * MixedMRFModel::log_det_ratio_yy_diag. + */ +inline double log_det_ratio_diag_kernel( + const arma::mat& covariance, size_t j, + double precision_curr_jj, double precision_prop_jj) { + double Uj2 = (precision_curr_jj - precision_prop_jj) / 2; + + double cc11 = covariance(j, j); + double cc12 = 1.0 - covariance(j, j) * Uj2; + double cc22 = Uj2 * Uj2 * covariance(j, j); + + return MY_LOG(std::abs(cc11 * cc22 - cc12 * cc12)); +} + +/** + * Extract the six reparameterization constants for a rank-2 off-diagonal MH + * proposal in precision space (Roverato move). Shared by GGMModel::get_constants + * and MixedMRFModel::get_precision_constants, which differ only in where the + * current precision entries come from (K vs −½K), passed in as scalars. + * + * @param cholesky Upper-triangular Cholesky factor of K. + * @param covariance Σ = K⁻¹. + * @param i, j Off-diagonal/diagonal indices (i < j). + * @param precision_ij, precision_jj Current K entries K(i,j), K(j,j). + * @return {Phi_q1q, Phi_q1q1, c2, Phi_q1q1, c4, c5}. + */ +inline std::array precision_proposal_constants( + const arma::mat& cholesky, const arma::mat& covariance, + size_t i, size_t j, double precision_ij, double precision_jj) { + double logdet = get_log_det(cholesky); + + double log_adj_ii = logdet + MY_LOG(std::abs(covariance(i, i))); + double log_adj_ij = logdet + MY_LOG(std::abs(covariance(i, j))); + double log_adj_jj = logdet + MY_LOG(std::abs(covariance(j, j))); + + double inv_sub_jj = compute_inv_submatrix_i(covariance, i, j, j); + double log_abs_inv_sub_jj = log_adj_ii + MY_LOG(std::abs(inv_sub_jj)); + + double Phi_q1q = (2 * std::signbit(covariance(i, j)) - 1) * MY_EXP( + (log_adj_ij - (log_adj_jj + log_abs_inv_sub_jj) / 2) + ); + double Phi_q1q1 = MY_EXP((log_adj_jj - log_abs_inv_sub_jj) / 2); + + std::array c{}; + c[0] = Phi_q1q; + c[1] = Phi_q1q1; + c[2] = precision_ij - Phi_q1q * Phi_q1q1; + c[3] = Phi_q1q1; + c[4] = precision_jj - Phi_q1q * Phi_q1q; + c[5] = c[4] + c[2] * c[2] / (c[3] * c[3]); + return c; +} + } // namespace cholesky_helpers diff --git a/src/models/ggm/ggm_model.cpp b/src/models/ggm/ggm_model.cpp index e42f0244..767853ce 100644 --- a/src/models/ggm/ggm_model.cpp +++ b/src/models/ggm/ggm_model.cpp @@ -162,27 +162,10 @@ arma::vec GGMModel::get_active_inv_mass() const { void GGMModel::get_constants(size_t i, size_t j) { - - double logdet_omega = cholesky_helpers::get_log_det(cholesky_of_precision_); - - double log_adj_omega_ii = logdet_omega + MY_LOG(std::abs(covariance_matrix_(i, i))); - double log_adj_omega_ij = logdet_omega + MY_LOG(std::abs(covariance_matrix_(i, j))); - double log_adj_omega_jj = logdet_omega + MY_LOG(std::abs(covariance_matrix_(j, j))); - - double inv_omega_sub_j1j1 = cholesky_helpers::compute_inv_submatrix_i(covariance_matrix_, i, j, j); - double log_abs_inv_omega_sub_jj = log_adj_omega_ii + MY_LOG(std::abs(inv_omega_sub_j1j1)); - double Phi_q1q = (2 * std::signbit(covariance_matrix_(i, j)) - 1) * MY_EXP( - (log_adj_omega_ij - (log_adj_omega_jj + log_abs_inv_omega_sub_jj) / 2) - ); - double Phi_q1q1 = MY_EXP((log_adj_omega_jj - log_abs_inv_omega_sub_jj) / 2); - - constants_[0] = Phi_q1q; - constants_[1] = Phi_q1q1; - constants_[2] = precision_matrix_(i, j) - Phi_q1q * Phi_q1q1; - constants_[3] = Phi_q1q1; - constants_[4] = precision_matrix_(j, j) - Phi_q1q * Phi_q1q; - constants_[5] = constants_[4] + constants_[2] * constants_[2] / (constants_[3] * constants_[3]); - + // GGM stores K directly, so the precision entries are precision_matrix_. + constants_ = cholesky_helpers::precision_proposal_constants( + cholesky_of_precision_, covariance_matrix_, i, j, + precision_matrix_(i, j), precision_matrix_(j, j)); } double GGMModel::constrained_diagonal(const double x) const { @@ -205,31 +188,18 @@ double GGMModel::log_density_impl(const arma::mat& omega, const arma::mat& phi) double GGMModel::log_det_ratio_edge(size_t i, size_t j) const { // Rank-2 matrix-determinant lemma: log|K_prop| - log|K_curr| where K_prop - // differs from K_curr at entries (i,j), (j,i), and (j,j). cc11, cc12, cc22 - // are the entries of the 2x2 update Gram matrix I + V^T Sigma U used by - // the lemma; |det(...)| is its determinant. - double Ui2 = precision_matrix_(i, j) - precision_proposal_(i, j); - double Uj2 = (precision_matrix_(j, j) - precision_proposal_(j, j)) / 2; - - double cc11 = covariance_matrix_(j, j); - double cc12 = 1 - (covariance_matrix_(i, j) * Ui2 - + covariance_matrix_(j, j) * Uj2); - double cc22 = Ui2 * Ui2 * covariance_matrix_(i, i) - + 2 * Ui2 * Uj2 * covariance_matrix_(i, j) - + Uj2 * Uj2 * covariance_matrix_(j, j); - - return MY_LOG(std::abs(cc11 * cc22 - cc12 * cc12)); + // differs from K_curr at entries (i,j), (j,i), and (j,j). GGM stores K + // directly, so the precision scalars are precision_matrix_ entries. + return cholesky_helpers::log_det_ratio_edge_kernel( + covariance_matrix_, i, j, + precision_matrix_(i, j), precision_proposal_(i, j), + precision_matrix_(j, j), precision_proposal_(j, j)); } double GGMModel::log_det_ratio_diag(size_t j) const { - // Rank-1 specialisation of log_det_ratio_edge (Ui2 = 0). - double Uj2 = (precision_matrix_(j, j) - precision_proposal_(j, j)) / 2; - - double cc11 = covariance_matrix_(j, j); - double cc12 = 1 - covariance_matrix_(j, j) * Uj2; - double cc22 = Uj2 * Uj2 * covariance_matrix_(j, j); - - return MY_LOG(std::abs(cc11 * cc22 - cc12 * cc12)); + return cholesky_helpers::log_det_ratio_diag_kernel( + covariance_matrix_, j, + precision_matrix_(j, j), precision_proposal_(j, j)); } double GGMModel::log_density_impl_edge(size_t i, size_t j) const { diff --git a/src/models/mixed/mixed_mrf_metropolis.cpp b/src/models/mixed/mixed_mrf_metropolis.cpp index fa5cc300..d8093fe2 100644 --- a/src/models/mixed/mixed_mrf_metropolis.cpp +++ b/src/models/mixed/mixed_mrf_metropolis.cpp @@ -182,31 +182,14 @@ double MixedMRFModel::update_pairwise_discrete(int i, int j, std::optional(i); + size_t uj = static_cast(j); + cont_constants_ = cholesky_helpers::precision_proposal_constants( + cholesky_of_precision_, covariance_continuous_, ui, uj, + -2.0 * pairwise_effects_continuous_(ui, uj), + -2.0 * pairwise_effects_continuous_(uj, uj)); } double MixedMRFModel::precision_constrained_diagonal(double x) const { @@ -231,37 +214,20 @@ double MixedMRFModel::precision_constrained_diagonal(double x) const { // ============================================================================= double MixedMRFModel::log_det_ratio_yy_edge(int i, int j) const { - // Rank-2 matrix-determinant lemma applied to the Kyy block. Mirrors - // GGMModel::log_det_ratio_edge but uses pairwise_effects_continuous_ - // (Kyy = -2 * pairwise_effects_continuous_) and covariance_continuous_. + // Kyy = -2 * pairwise_effects_continuous_; delegate to the shared kernel. size_t ui = static_cast(i); size_t uj = static_cast(j); - double precision_ij_curr = -2.0 * pairwise_effects_continuous_(ui, uj); - double precision_jj_curr = -2.0 * pairwise_effects_continuous_(uj, uj); - double Ui = precision_ij_curr - precision_proposal_(ui, uj); - double Uj = (precision_jj_curr - precision_proposal_(uj, uj)) / 2.0; - - double cc11 = covariance_continuous_(uj, uj); - double cc12 = 1.0 - (covariance_continuous_(ui, uj) * Ui + - covariance_continuous_(uj, uj) * Uj); - double cc22 = Ui * Ui * covariance_continuous_(ui, ui) + - 2.0 * Ui * Uj * covariance_continuous_(ui, uj) + - Uj * Uj * covariance_continuous_(uj, uj); - - return MY_LOG(std::abs(cc11 * cc22 - cc12 * cc12)); + return cholesky_helpers::log_det_ratio_edge_kernel( + covariance_continuous_, ui, uj, + -2.0 * pairwise_effects_continuous_(ui, uj), precision_proposal_(ui, uj), + -2.0 * pairwise_effects_continuous_(uj, uj), precision_proposal_(uj, uj)); } double MixedMRFModel::log_det_ratio_yy_diag(int i) const { - // Rank-1 specialisation of log_det_ratio_yy_edge (Ui = 0). size_t ui = static_cast(i); - double precision_ii_curr = -2.0 * pairwise_effects_continuous_(ui, ui); - double Uj = (precision_ii_curr - precision_proposal_(ui, ui)) / 2.0; - - double cc11 = covariance_continuous_(ui, ui); - double cc12 = 1.0 - covariance_continuous_(ui, ui) * Uj; - double cc22 = Uj * Uj * covariance_continuous_(ui, ui); - - return MY_LOG(std::abs(cc11 * cc22 - cc12 * cc12)); + return cholesky_helpers::log_det_ratio_diag_kernel( + covariance_continuous_, ui, + -2.0 * pairwise_effects_continuous_(ui, ui), precision_proposal_(ui, ui)); }