From 6f7d80100e2dfdac0b47144b0c259f69a3358753 Mon Sep 17 00:00:00 2001 From: Maarten Marsman Date: Thu, 4 Jun 2026 16:27:04 +0200 Subject: [PATCH 1/2] refactor(ggm): route gradient-cache invalidation through a named mutator The MH within-K edge/diag tuning sweep and the two edge-selection accept branches each flipped constraint_dirty_ = true; theta_valid_ = false; inline -- three copies of the same 'precision/graph changed, NUTS caches are stale' rule. Extract invalidate_gradient_cache() and call it at all three sites. set_determinant_tilt is intentionally left as-is: it only sets constraint_dirty_ (gradient-engine rebuild) and must NOT also clear theta_valid_, so it does not route through the mutator. Behavior-preserving: identical flag writes. Verified bitwise -- adaptive- Metropolis GGM raw draws (incl. the edge-selection path) digest-identical. --- src/models/ggm/ggm_model.cpp | 9 +++------ src/models/ggm/ggm_model.h | 15 +++++++++++++++ 2 files changed, 18 insertions(+), 6 deletions(-) diff --git a/src/models/ggm/ggm_model.cpp b/src/models/ggm/ggm_model.cpp index 03fa5029..501d7564 100644 --- a/src/models/ggm/ggm_model.cpp +++ b/src/models/ggm/ggm_model.cpp @@ -482,8 +482,7 @@ void GGMModel::update_edge_indicator_parameter_pair(size_t i, size_t j) { cholesky_update_after_edge(omega_ij_old, omega_jj_old, i, j); - constraint_dirty_ = true; - theta_valid_ = false; + invalidate_gradient_cache(); } } else { @@ -539,8 +538,7 @@ void GGMModel::update_edge_indicator_parameter_pair(size_t i, size_t j) { cholesky_update_after_edge(omega_ij_old, omega_jj_old, i, j); - constraint_dirty_ = true; - theta_valid_ = false; + invalidate_gradient_cache(); } } } @@ -641,8 +639,7 @@ void GGMModel::tune_proposal_sd(int iteration, const WarmupSchedule& schedule) { } // Invalidate gradient cache after MH updates - constraint_dirty_ = true; - theta_valid_ = false; + invalidate_gradient_cache(); } void GGMModel::refresh_cholesky() { diff --git a/src/models/ggm/ggm_model.h b/src/models/ggm/ggm_model.h index 2cbd3417..9886064d 100644 --- a/src/models/ggm/ggm_model.h +++ b/src/models/ggm/ggm_model.h @@ -680,6 +680,21 @@ class GGMModel : public BaseModel { mutable arma::vec theta_; public: + /** + * Mark the NUTS gradient caches stale after a change to the precision + * matrix or the edge set. Both the constraint structure (active dimension + * depends on the edge indicators) and the cached theta parameterization + * become invalid, so the MH within-K, edge-selection, and proposal-SD + * tuning paths all funnel their invalidation through here rather than + * repeating the two flag writes. (set_determinant_tilt is intentionally + * NOT routed through this: it only needs a gradient-engine rebuild, not a + * theta invalidation.) + */ + void invalidate_gradient_cache() { + constraint_dirty_ = true; + theta_valid_ = false; + } + /** * Rebuild the constraint structure and gradient engine from current * edge indicators. Called lazily before gradient evaluation. From 75f9480a209eb6ed78b8317bbc825a177a11eea1 Mon Sep 17 00:00:00 2001 From: Maarten Marsman Date: Thu, 4 Jun 2026 16:31:00 +0200 Subject: [PATCH 2/2] refactor(ggm): delegate the raw-data constructor to the suf-stat constructor The two GGMModel constructors shared ~18 lines of identical member initialization (priors, precision/Cholesky/covariance identity warm-start, edge indicators, proposal SDs, sparse-graph detection, MLE warm-start). The raw-data ctor differs only in deriving n (n-1 effective df) and S = X'X from the observations and in retaining observations_ for imputation. Delegate raw-data -> suf-stat; the only raw-specific work (storing observations_ when na_impute) moves to the body. initialize_precision_from_mle reads suf_stat_/ n_/edge_indicators_, never observations_, so running it in the delegated ctor before the body sets observations_ is safe. has_missing_/missing_index_ are set by a separate setter, not either ctor, so they are unaffected. Verified bitwise: na_action=impute GGM fits (edge-sel + no-edge, exercising the raw-data ctor and the in-chain imputation RNG) plus a complete-data fit are digest-identical pre/post. --- src/models/ggm/ggm_model.h | 38 +++++++++++++++----------------------- 1 file changed, 15 insertions(+), 23 deletions(-) diff --git a/src/models/ggm/ggm_model.h b/src/models/ggm/ggm_model.h index 9886064d..2e96f034 100644 --- a/src/models/ggm/ggm_model.h +++ b/src/models/ggm/ggm_model.h @@ -47,30 +47,22 @@ class GGMModel : public BaseModel { std::unique_ptr interaction_prior, std::unique_ptr diagonal_prior, const bool na_impute = false - ) : n_(observations.n_rows - 1), // centered data has n-1 effective df - p_(observations.n_cols), - dim_((p_ * (p_ + 1)) / 2), - suf_stat_(observations.t() * observations), - inclusion_probability_(inclusion_probability), - edge_selection_(edge_selection), - interaction_prior_(std::move(interaction_prior)), - diagonal_prior_(std::move(diagonal_prior)), - precision_matrix_(arma::eye(p_, p_)), - cholesky_of_precision_(arma::eye(p_, p_)), - inv_cholesky_of_precision_(arma::eye(p_, p_)), - covariance_matrix_(arma::eye(p_, p_)), - edge_indicators_(initial_edge_indicators), - vectorized_parameters_(dim_), - vectorized_indicator_parameters_(edge_selection_ ? dim_ : 0), - proposal_sds_(arma::mat(dim_, 1, arma::fill::ones) * 0.25), - num_pairwise_(p_ * (p_ - 1) / 2), - observations_(na_impute ? observations : arma::mat()), - precision_proposal_(arma::mat(p_, p_, arma::fill::none)) + ) // Delegate to the sufficient-statistics constructor: the raw-data path + // differs only in deriving n (centered data has n-1 effective df) and + // S = X'X from the observations, and in retaining the observations for + // missing-data imputation. All member init and the MLE warm-start live + // in the delegated constructor. + : GGMModel(observations.n_rows - 1, + observations.t() * observations, + inclusion_probability, + initial_edge_indicators, + edge_selection, + std::move(interaction_prior), + std::move(diagonal_prior)) { - int num_edges = arma::accu(edge_indicators_) / 2; - int max_edges = static_cast(p_ * (p_ - 1) / 2); - has_sparse_graph_ = !edge_selection_ && (num_edges < max_edges); - initialize_precision_from_mle(); + if (na_impute) { + observations_ = observations; + } } /**