diff --git a/src/mcmc/algorithms/metropolis.cpp b/src/mcmc/algorithms/metropolis.cpp
index e57099e2..995a35f7 100644
--- a/src/mcmc/algorithms/metropolis.cpp
+++ b/src/mcmc/algorithms/metropolis.cpp
@@ -23,3 +23,23 @@ StepResult metropolis_step(
 
   return {State, accept_prob};
 }
+
+
+StepResult metropolis_step_cached(
+    double current_state,
+    double step_size,
+    double log_post_current,
+    const std::function<double(double)>& log_post_proposed,
+    SafeRNG& rng
+) {
+  double proposed_state = rnorm(rng, current_state, step_size);
+  double log_accept = log_post_proposed(proposed_state) - log_post_current;
+  double accept_prob = std::min(1.0, MY_EXP(log_accept));
+
+  double state = (runif(rng) < accept_prob) ? proposed_state : current_state;
+
+  arma::vec State(1);
+  State[0] = state;
+
+  return {State, accept_prob};
+}
diff --git a/src/mcmc/algorithms/metropolis.h b/src/mcmc/algorithms/metropolis.h
index 6d5e1252..ec75f0af 100644
--- a/src/mcmc/algorithms/metropolis.h
+++ b/src/mcmc/algorithms/metropolis.h
@@ -24,3 +24,26 @@ StepResult metropolis_step(
     const std::function<double(double)>& log_post,
     SafeRNG& rng
 );
+
+/**
+ * Performs one Random Walk Metropolis step with a cached current-state value
+ *
+ * Variant of metropolis_step() for callers that maintain the current-state
+ * log-posterior across proposals: only the proposed state is evaluated. The
+ * random-number stream matches metropolis_step() (one normal draw, one
+ * uniform draw).
+ *
+ * @param current_state     Current scalar parameter value
+ * @param step_size         Standard deviation of the Gaussian proposal
+ * @param log_post_current  Log-posterior value at current_state
+ * @param log_post_proposed Log-posterior function evaluated at the proposal
+ * @param rng               Thread-safe random number generator
+ * @return StepResult with accepted state (1-element vector) and acceptance probability
+ */
+StepResult metropolis_step_cached(
+    double current_state,
+    double step_size,
+    double log_post_current,
+    const std::function<double(double)>& log_post_proposed,
+    SafeRNG& rng
+);
diff --git a/src/models/bgmCompare/bgmCompare_helper.cpp b/src/models/bgmCompare/bgmCompare_helper.cpp
index b0b47860..be855d77 100644
--- a/src/models/bgmCompare/bgmCompare_helper.cpp
+++ b/src/models/bgmCompare/bgmCompare_helper.cpp
@@ -1,10 +1,74 @@
 #include <RcppArmadillo.h>
 #include <cmath>
 #include "models/bgmCompare/bgmCompare_helper.h"
+#include "models/bgmCompare/bgmCompare_state.h"
 #include "utils/common_helpers.h"
 
 
 
+// Converts the integer observations to double, whole-matrix and per group.
+// Weight and residual members are left untouched; call
+// rebuild_sweep_state_weights() afterwards to make the state consistent.
+void initialize_sweep_state_observations(
+    CompareSweepState& state,
+    const arma::imat& observations,
+    const arma::imat& group_indices,
+    const int num_groups
+) {
+  state.obs_double_all = arma::conv_to<arma::mat>::from(observations);
+  state.obs_double.resize(num_groups);
+  for (int g = 0; g < num_groups; g++) {
+    const int r0 = group_indices(g, 0);
+    const int r1 = group_indices(g, 1);
+    state.obs_double[g] = state.obs_double_all.rows(r0, r1);
+  }
+
+  const int num_variables = observations.n_cols;
+  state.log_normalizer.zeros(num_variables, num_groups);
+  state.normalizer_valid.zeros(num_variables);
+}
+
+
+
+// Recomputes the per-group effective pairwise weights from the current
+// pairwise effects and inclusion indicators, and the residual matrices as
+// one matrix product per group. Invalidates the normalizer cache.
+void rebuild_sweep_state_weights(
+    CompareSweepState& state,
+    const arma::mat& pairwise_effects,
+    const arma::imat& pairwise_effect_indices,
+    const arma::imat& inclusion_indicator,
+    const arma::mat& projection,
+    const int num_groups
+) {
+  const int num_variables = inclusion_indicator.n_rows;
+  state.pairwise_group.resize(num_groups);
+  state.residual.resize(num_groups);
+
+  for (int g = 0; g < num_groups; g++) {
+    const arma::vec proj_g = projection.row(g).t();
+
+    arma::mat& pairwise_g = state.pairwise_group[g];
+    pairwise_g.zeros(num_variables, num_variables);
+    for (int v = 0; v < num_variables - 1; v++) {
+      for (int u = v + 1; u < num_variables; u++) {
+        double w = compute_group_pairwise_effects(
+          v, u, num_groups, pairwise_effects, pairwise_effect_indices,
+          inclusion_indicator, proj_g
+        );
+        pairwise_g(v, u) = w;
+        pairwise_g(u, v) = w;
+      }
+    }
+
+    state.residual[g] = state.obs_double[g] * pairwise_g;
+  }
+
+  state.normalizer_valid.zeros();
+}
+
+
+
 // Computes group-specific main effects for a given variable (bgmCompare model).
 //
 // For a variable, the main-effect parameters are stored with:
diff --git a/src/models/bgmCompare/bgmCompare_logp_and_grad.cpp b/src/models/bgmCompare/bgmCompare_logp_and_grad.cpp
index 38b52ace..eb34cd41 100644
--- a/src/models/bgmCompare/bgmCompare_logp_and_grad.cpp
+++ b/src/models/bgmCompare/bgmCompare_logp_and_grad.cpp
@@ -771,21 +771,19 @@ std::pair<double, arma::vec> logp_and_gradient(
 //  - For each group:
 //    * Construct group-specific main effects for the selected variable
 //      with `compute_group_main_effects()`.
-//    * Construct group-specific pairwise effects for the variable.
 //    * Add linear contributions from sufficient statistics.
-//    * Subtract log normalizing constants from the group-specific likelihood.
+//    * Subtract log normalizing constants from the group-specific likelihood,
+//      reading the variable's rest scores from the maintained residual
+//      matrices (pairwise effects are fixed during main-effect updates).
 //  - Add prior contribution:
-//    * Logistic–Beta prior for baseline (h == 0).
+//    * Threshold prior for baseline (h == 0).
 //    * Difference prior for group differences (h > 0), if included.
 //
 // Inputs:
 //  - main_effects: Matrix of main-effect parameters (rows = categories, cols = groups).
-//  - pairwise_effects: Matrix of pairwise-effect parameters (rows = pairs, cols = groups).
 //  - main_effect_indices: Index ranges [row_start,row_end] for each variable.
-//  - pairwise_effect_indices: Lookup table mapping (var1,var2) → row in pairwise_effects.
 //  - projection: Group projection matrix (num_groups × (num_groups − 1)).
-//  - observations: Observation matrix (persons × variables).
-//  - group_indices: Row ranges [start,end] for each group in observations.
+//  - residual_groups: Per-group rest-score matrices (persons × variables).
 //  - num_categories: Number of categories per variable.
 //  - counts_per_category_group: Per-group category counts (for ordinal variables).
 //  - blume_capel_stats_group: Per-group sufficient statistics (for Blume–Capel variables).
@@ -793,12 +791,15 @@ std::pair<double, arma::vec> logp_and_gradient(
 //  - inclusion_indicator: Symmetric binary matrix of active variables (diag) and pairs (off-diag).
 //  - is_ordinal_variable: Indicator (1 = ordinal, 0 = Blume–Capel).
 //  - baseline_category: Reference categories for Blume–Capel variables.
-//  - main_alpha, main_beta: Hyperparameters for Beta priors on main effects.
-//  - difference_scale: Scale parameter of the difference prior.
 //  - variable: Index of the variable of interest.
 //  - category: Category index (only used if variable is ordinal).
 //  - par: Parameter index (0 = linear, 1 = quadratic; used for Blume–Capel).
 //  - h: Column index (0 = overall baseline, >0 = group difference).
+//  - difference_prior, threshold_prior: Parameter priors.
+//  - normalizers_in: If non-null, cached per-group log-normalizer sums for
+//    the variable (length G); the normalizer computation is skipped.
+//  - normalizers_out: If non-null, receives the computed per-group
+//    log-normalizer sums (length G). Ignored when normalizers_in is given.
 //
 // Returns:
 //  - The scalar log pseudoposterior contribution of the selected parameter.
@@ -810,12 +811,9 @@ std::pair<double, arma::vec> logp_and_gradient(
 //  - Consistent with the full `log_pseudoposterior()` for bgmCompare.
 double log_pseudoposterior_main_component(
     const arma::mat& main_effects,
-    const arma::mat& pairwise_effects,
     const arma::imat& main_effect_indices,
-    const arma::imat& pairwise_effect_indices,
     const arma::mat& projection,
-    const arma::imat& observations,
-    const arma::imat& group_indices,
+    const std::vector<arma::mat>& residual_groups,
     const arma::ivec& num_categories,
     const std::vector<arma::imat>& counts_per_category_group,
     const std::vector<arma::imat>& blume_capel_stats_group,
@@ -828,18 +826,16 @@ double log_pseudoposterior_main_component(
     int par, // for Blume-Capel variables only
     int h, // Overall = 0, differences are 1,2,...
     const BaseParameterPrior& difference_prior,
-    const BaseParameterPrior& threshold_prior
+    const BaseParameterPrior& threshold_prior,
+    const arma::vec* normalizers_in,
+    arma::vec* normalizers_out
 ) {
   if(h > 0 && inclusion_indicator(variable, variable) == 0) {
     return 0.0; // No contribution if differences not included
   }
 
-  const int num_variables = observations.n_cols;
   double log_pp = 0.0;
 
-  // group-specific pairwise weights of `variable` with all other variables
-  arma::vec pairwise_col(num_variables, arma::fill::zeros);
-
   // --- per group ---
   for (int group = 0; group < num_groups; ++group) {
     const arma::imat& counts_per_category = counts_per_category_group[group];
@@ -847,20 +843,11 @@ double log_pseudoposterior_main_component(
 
     const arma::vec proj_g = projection.row(group).t(); // length = num_groups-1
 
-    // ---- build group-specific main & pairwise effects ----
+    // ---- build group-specific main effects ----
     arma::vec me = compute_group_main_effects(
       variable, num_groups, main_effects, main_effect_indices, proj_g
     );
 
-    // pairwise weights with the other variables; entry `variable` stays zero
-    for (int u = 0; u < num_variables; u++) {
-      if(u == variable) continue;
-      pairwise_col(u) = compute_group_pairwise_effects(
-        variable, u, num_groups, pairwise_effects, pairwise_effect_indices,
-        inclusion_indicator, proj_g
-      );
-    }
-
     // ---- data contribution pseudolikelihood (linear terms) ----
     if (is_ordinal_variable(variable)) {
       log_pp += static_cast<double>(counts_per_category(category, variable)) *
@@ -869,13 +856,17 @@ double log_pseudoposterior_main_component(
       log_pp += static_cast<double>(blume_capel_stats(par, variable)) * me(par);
     }
 
-    // ---- data contribution pseudolikelihood (quadratic terms) ----
-    const int r0 = group_indices(group, 0);
-    const int r1 = group_indices(group, 1);
-    const arma::mat obs = arma::conv_to<arma::mat>::from(observations.rows(r0, r1));
-
     // ---- pseudolikelihood normalizing constants (per variable) ----
-    const arma::vec rest_score = obs * pairwise_col;
+    // The caller may supply the variable's cached per-group log-normalizer
+    // sums for the current state.
+    if (normalizers_in != nullptr) {
+      log_pp -= (*normalizers_in)(group);
+      continue;
+    }
+
+    // Rest scores come from the maintained residual matrix; the pairwise
+    // effects do not change during main-effect updates.
+    const arma::vec rest_score = residual_groups[group].col(variable);
     const int num_cats = num_categories(variable);
 
     // bound to stabilize exp; clamp at 0 (the reference-category exponent) so
@@ -900,7 +891,11 @@ double log_pseudoposterior_main_component(
     }
 
     // - sum_i [ bound_i + log denom_i ]
-    log_pp -= arma::accu(bound + ARMA_MY_LOG(denom));
+    const double normalizer = arma::accu(bound + ARMA_MY_LOG(denom));
+    log_pp -= normalizer;
+    if (normalizers_out != nullptr) {
+      (*normalizers_out)(group) = normalizer;
+    }
   }
 
   // ---- priors ----
@@ -963,6 +958,11 @@ double log_pseudoposterior_main_component(
 //  - variable1, variable2: Indices of the variable pair.
 //  - h: Column index (0 = baseline, > 0 = group difference).
 //  - delta: Parameter change (proposed - current).
+//  - normalizers_in: If non-null, cached per-group log-normalizer sums for
+//    the pair's endpoint variables (G x 2: column 0 = variable1, column 1 =
+//    variable2); the normalizer computation is skipped.
+//  - normalizers_out: If non-null, receives the computed per-group
+//    log-normalizer sums (G x 2). Ignored when normalizers_in is given.
 //
 // Returns:
 //  - The log pseudoposterior value at the proposed state.
@@ -990,7 +990,9 @@ double log_pseudoposterior_pair_component(
     int h,
     double delta,
     const BaseParameterPrior& interaction_prior,
-    const BaseParameterPrior& difference_prior
+    const BaseParameterPrior& difference_prior,
+    const arma::mat* normalizers_in,
+    arma::mat* normalizers_out
 ) {
   if(h > 0 && inclusion_indicator(variable1, variable2) == 0) {
     return 0.0;
@@ -1008,6 +1010,25 @@ double log_pseudoposterior_pair_component(
   for (int group = 0; group < num_groups; ++group) {
     const arma::vec proj_g = projection.row(group).t();
 
+    // ---- data contribution pseudolikelihood ----
+    const arma::mat& pairwise_stats = pairwise_stats_group[group];
+    const double suff_pair = pairwise_stats(variable1, variable2);
+
+    if(h == 0) {
+      log_pp += 2.0 * suff_pair * proposed_value;
+    } else {
+      log_pp += 2.0 * suff_pair * proj_g(h-1) * proposed_value;
+    }
+
+    // ---- pseudolikelihood normalizing constants ----
+    // The caller may supply cached per-group log-normalizer sums for the
+    // pair's endpoint variables at the current state.
+    if (normalizers_in != nullptr) {
+      log_pp -= (*normalizers_in)(group, 0);
+      log_pp -= (*normalizers_in)(group, 1);
+      continue;
+    }
+
     // Compute group-specific delta: how much pairwise_group(var1,var2) changes for this group
     double delta_g = (h == 0) ? delta : delta * proj_g(h - 1);
 
@@ -1020,19 +1041,9 @@ double log_pseudoposterior_pair_component(
       main_group(v, arma::span(0, me.n_elem - 1)) = me.t();
     }
 
-    // ---- data contribution pseudolikelihood ----
-    const arma::mat& pairwise_stats = pairwise_stats_group[group];
-    const double suff_pair = pairwise_stats(variable1, variable2);
-
-    if(h == 0) {
-      log_pp += 2.0 * suff_pair * proposed_value;
-    } else {
-      log_pp += 2.0 * suff_pair * proj_g(h-1) * proposed_value;
-    }
-
-    // ---- pseudolikelihood normalizing constants (using residual matrix + delta) ----
     const arma::mat& obs_g = obs_double_groups[group];
 
+    int slot = 0;
     for (int v : {variable1, variable2}) {
       const int num_cats = num_categories(v);
       const int other = (v == variable1) ? variable2 : variable1;
@@ -1055,7 +1066,12 @@ double log_pseudoposterior_pair_component(
         denom = compute_denom_blume_capel(rest_score, lin_effect, quad_effect, ref, num_cats, bound);
       }
 
-      log_pp -= arma::accu(bound + ARMA_MY_LOG(denom));
+      const double normalizer = arma::accu(bound + ARMA_MY_LOG(denom));
+      log_pp -= normalizer;
+      if (normalizers_out != nullptr) {
+        (*normalizers_out)(group, slot) = normalizer;
+      }
+      ++slot;
     }
   }
 
@@ -1070,147 +1086,6 @@ double log_pseudoposterior_pair_component(
 
 
 
-// Computes the log-ratio of pseudolikelihood normalizing constants
-// for a single variable under current vs. proposed parameters (bgmCompare model).
-//
-// This function is used in Metropolis–Hastings updates for main-effect parameters.
-// It evaluates how the normalizing constant (denominator of the pseudolikelihood)
-// changes when switching from the current to the proposed parameter values.
-//
-// Procedure:
-//  - For each group:
-//    * Construct group-specific main effects (current vs. proposed).
-//    * Construct group-specific pairwise weights for the variable.
-//    * Compute residual scores for observations under both models.
-//    * Calculate denominators with stability bounds (ordinal vs. Blume–Capel cases).
-//    * Accumulate the log-ratio contribution across all observations.
-//
-// Inputs:
-//  - current_main_effects, proposed_main_effects: Matrices of main-effect parameters
-//    (rows = categories, cols = groups).
-//  - current_pairwise_effects, proposed_pairwise_effects: Matrices of pairwise-effect parameters
-//    (rows = pairs, cols = groups).
-//  - main_effect_indices: Index ranges [row_start,row_end] for each variable.
-//  - pairwise_effect_indices: Lookup table mapping (var1,var2) → row in pairwise_effects.
-//  - projection: Group projection matrix (num_groups × (num_groups − 1)).
-//  - observations: Observation matrix (persons × variables).
-//  - group_indices: Row ranges [start,end] for each group in observations.
-//  - num_categories: Number of categories per variable.
-//  - num_groups: Number of groups.
-//  - inclusion_indicator: Symmetric binary matrix of active variables (diag) and pairs (off-diag).
-//  - is_ordinal_variable: Indicator (1 = ordinal, 0 = Blume–Capel).
-//  - baseline_category: Reference categories for Blume–Capel variables.
-//  - variable: Index of the variable being updated.
-//
-// Returns:
-//  - The scalar log-ratio of pseudolikelihood constants
-//    (current model vs. proposed model).
-//
-// Notes:
-//  - For ordinal variables, denominators include exp(-bound) and category terms.
-//  - For Blume–Capel variables, denominators use linear/quadratic scores
-//    with baseline centering.
-//  - Stability bounds (`bound_current`, `bound_proposed`) are applied to avoid overflow.
-double log_ratio_pseudolikelihood_constant_variable(
-    const arma::mat& current_main_effects,
-    const arma::mat& current_pairwise_effects,
-    const arma::mat& proposed_main_effects,
-    const arma::mat& proposed_pairwise_effects,
-    const arma::imat& main_effect_indices,
-    const arma::imat& pairwise_effect_indices,
-    const arma::mat& projection,
-    const arma::imat& observations,
-    const arma::imat& group_indices,
-    const arma::ivec& num_categories,
-    const int num_groups,
-    const arma::imat& inclusion_indicator,
-    const arma::uvec& is_ordinal_variable,
-    const arma::ivec& baseline_category,
-    const int variable
-) {
-  const int num_cats = num_categories(variable);
-  const int num_variables = observations.n_cols;
-
-  double log_ratio = 0.0;
-
-  // --- per group ---
-  for (int group = 0; group < num_groups; ++group) {
-    const arma::vec proj_g = projection.row(group).t();
-
-    // --- group-specific main effects (current/proposed) ---
-    const arma::vec main_current = compute_group_main_effects(
-      variable, num_groups, current_main_effects, main_effect_indices, proj_g
-    );
-    const arma::vec main_proposed = compute_group_main_effects(
-      variable, num_groups, proposed_main_effects, main_effect_indices, proj_g
-    );
-
-    // --- group-specific pairwise effects for this variable (column) ---
-    arma::vec weights_current(num_variables, arma::fill::zeros);
-    arma::vec weights_proposed(num_variables, arma::fill::zeros);
-    for (int u = 0; u < num_variables; ++u) {
-      if (u == variable) continue;
-      weights_current(u) = compute_group_pairwise_effects(
-        variable, u, num_groups, current_pairwise_effects,
-        pairwise_effect_indices, inclusion_indicator, proj_g
-      );
-      weights_proposed(u) = compute_group_pairwise_effects(
-        variable, u, num_groups, proposed_pairwise_effects,
-        pairwise_effect_indices, inclusion_indicator, proj_g
-      );
-    }
-
-    // --- group observations and rest scores ---
-    const int r0 = group_indices(group, 0);
-    const int r1 = group_indices(group, 1);
-    const arma::mat obs = arma::conv_to<arma::mat>::from(observations.rows(r0, r1));
-
-    const arma::vec rest_current = obs * weights_current;
-    const arma::vec rest_proposed = obs * weights_proposed;
-
-    // --- denominators with stability bounds ---
-    arma::vec bound_current;
-    arma::vec bound_proposed;
-    arma::vec denom_current(rest_current.n_elem, arma::fill::zeros);
-    arma::vec denom_proposed(rest_proposed.n_elem, arma::fill::zeros);
-
-    if (is_ordinal_variable (variable)) {
-      bound_current = arma::clamp(rest_current * num_cats, 0.0, arma::datum::inf);
-      bound_proposed = arma::clamp(rest_proposed * num_cats, 0.0, arma::datum::inf);
-
-      denom_current += compute_denom_ordinal(
-        rest_current, main_current, bound_current
-      );
-      denom_proposed += compute_denom_ordinal(
-        rest_proposed, main_proposed, bound_proposed
-      );
-    } else {
-      // Binary or categorical variable: linear + quadratic score
-      const int ref_cat = baseline_category (variable);
-      bound_current = rest_current * num_cats;
-      bound_proposed = rest_proposed * num_cats;
-
-      denom_current = compute_denom_blume_capel(
-        rest_current, main_current (0), main_current (1),
-        ref_cat, num_cats, /*Updated in place:*/bound_current
-      );
-
-      denom_proposed = compute_denom_blume_capel(
-        rest_proposed, main_proposed (0), main_proposed (1),
-        ref_cat, num_cats, /*Updated in place:*/bound_proposed
-      );
-    }
-
-    // --- accumulate contribution ---
-    log_ratio += arma::accu((bound_current - bound_proposed) +
-      ARMA_MY_LOG(denom_current) - ARMA_MY_LOG(denom_proposed));
-  }
-
-  return log_ratio;
-}
-
-
-
 // Computes the log pseudolikelihood ratio for updating a single main-effect parameter (bgmCompare model).
 //
 // This function is used in Metropolis–Hastings updates for main effects.
@@ -1223,23 +1098,20 @@ double log_ratio_pseudolikelihood_constant_variable(
 //    * Compute group-specific main effects for the variable (current vs. proposed).
 //    * Add contributions from observed sufficient statistics
 //      (category counts or Blume–Capel stats).
-//  - Add the ratio of pseudolikelihood normalizing constants by calling
-//    `log_ratio_pseudolikelihood_constant_variable()`.
+//    * Add the ratio of pseudolikelihood normalizing constants. The pairwise
+//      effects are identical in both states, so both sides share the rest
+//      scores held in the maintained residual matrices.
 //
 // Inputs:
 //  - current_main_effects: Matrix of main-effect parameters (current state).
 //  - proposed_main_effects: Matrix of main-effect parameters (candidate state).
-//  - current_pairwise_effects: Matrix of pairwise-effect parameters (fixed at current state).
 //  - main_effect_indices: Index ranges [row_start,row_end] for each variable.
-//  - pairwise_effect_indices: Lookup table mapping (var1,var2) → row in pairwise_effects.
 //  - projection: Group projection matrix (num_groups × (num_groups − 1)).
-//  - observations: Observation matrix (persons × variables).
-//  - group_indices: Row ranges [start,end] for each group in observations.
+//  - residual_groups: Per-group rest-score matrices (persons × variables).
 //  - num_categories: Number of categories per variable.
 //  - counts_per_category_group: Per-group category counts (for ordinal variables).
 //  - blume_capel_stats_group: Per-group sufficient statistics (for Blume–Capel variables).
 //  - num_groups: Number of groups.
-//  - inclusion_indicator: Symmetric binary matrix of active variables (diag) and pairs (off-diag).
 //  - is_ordinal_variable: Indicator (1 = ordinal, 0 = Blume–Capel).
 //  - baseline_category: Reference categories for Blume–Capel variables.
 //  - variable: Index of the variable being updated.
@@ -1248,34 +1120,26 @@ double log_ratio_pseudolikelihood_constant_variable(
 //  - The scalar log pseudolikelihood ratio (proposed vs. current).
 //
 // Notes:
-//  - A temporary copy of `inclusion_indicator` is made to ensure the
-//    variable’s self-term (diagonal entry) is included.
 //  - Only the variable under update changes between current and proposed states;
 //    all other variables and pairwise effects remain fixed.
 //  - This function does not add prior contributions — only pseudolikelihood terms.
 double log_pseudolikelihood_ratio_main(
     const arma::mat& current_main_effects,
     const arma::mat& proposed_main_effects,
-    const arma::mat& current_pairwise_effects,
     const arma::imat& main_effect_indices,
-    const arma::imat& pairwise_effect_indices,
-    const arma::mat&  projection,
-    const arma::imat& observations,
-    const arma::imat& group_indices,
+    const arma::mat& projection,
+    const std::vector<arma::mat>& residual_groups,
     const arma::ivec& num_categories,
     const std::vector<arma::imat>& counts_per_category_group,
     const std::vector<arma::imat>& blume_capel_stats_group,
     const int num_groups,
-    const arma::imat& inclusion_indicator,
     const arma::uvec& is_ordinal_variable,
     const arma::ivec& baseline_category,
     const int variable
 ) {
+  const int num_cats = num_categories(variable);
   double lr = 0.0;
-  arma::imat tmp_ind = inclusion_indicator;
-  tmp_ind(variable, variable) = 1; // Ensure self-interaction is included
 
-  // Add data contribution (group-specific parameters via projection)
   for (int g = 0; g < num_groups; ++g) {
     const arma::vec proj_g = projection.row(g).t();
 
@@ -1286,9 +1150,9 @@ double log_pseudolikelihood_ratio_main(
       variable, num_groups, proposed_main_effects, main_effect_indices, proj_g
     );
 
+    // Add data contribution (group-specific parameters via projection)
     if (is_ordinal_variable(variable)) {
       const arma::imat& num_obs = counts_per_category_group[g];
-      const int num_cats = num_categories(variable);
       for (int c = 0; c < num_cats; ++c) {
         lr += (main_prop(c) - main_cur(c)) * static_cast<double>(num_obs(c, variable));
       }
@@ -1297,16 +1161,44 @@ double log_pseudolikelihood_ratio_main(
       lr += (main_prop(0) - main_cur(0)) * static_cast<double>(suff(0, variable));
       lr += (main_prop(1) - main_cur(1)) * static_cast<double>(suff(1, variable));
     }
-  }
 
-  // Add ratio of normalizing constants
-  lr += log_ratio_pseudolikelihood_constant_variable(
-    current_main_effects, current_pairwise_effects, proposed_main_effects,
-    /* same */ current_pairwise_effects, main_effect_indices,
-    pairwise_effect_indices, projection, observations, group_indices,
-    num_categories, num_groups, tmp_ind, is_ordinal_variable,
-    baseline_category, variable
-  );
+    // Add ratio of normalizing constants; rest scores are shared between
+    // the two states because the pairwise effects are unchanged.
+    const arma::vec rest_score = residual_groups[g].col(variable);
+
+    arma::vec bound_current;
+    arma::vec bound_proposed;
+    arma::vec denom_current(rest_score.n_elem, arma::fill::zeros);
+    arma::vec denom_proposed(rest_score.n_elem, arma::fill::zeros);
+
+    if (is_ordinal_variable(variable)) {
+      bound_current = arma::clamp(rest_score * num_cats, 0.0, arma::datum::inf);
+      bound_proposed = bound_current;
+
+      denom_current += compute_denom_ordinal(
+        rest_score, main_cur, bound_current
+      );
+      denom_proposed += compute_denom_ordinal(
+        rest_score, main_prop, bound_proposed
+      );
+    } else {
+      const int ref_cat = baseline_category(variable);
+      bound_current = rest_score * num_cats;
+      bound_proposed = rest_score * num_cats;
+
+      denom_current = compute_denom_blume_capel(
+        rest_score, main_cur(0), main_cur(1),
+        ref_cat, num_cats, /*Updated in place:*/bound_current
+      );
+      denom_proposed = compute_denom_blume_capel(
+        rest_score, main_prop(0), main_prop(1),
+        ref_cat, num_cats, /*Updated in place:*/bound_proposed
+      );
+    }
+
+    lr += arma::accu((bound_current - bound_proposed) +
+      ARMA_MY_LOG(denom_current) - ARMA_MY_LOG(denom_proposed));
+  }
 
   return lr;
 }
@@ -1324,9 +1216,10 @@ double log_pseudolikelihood_ratio_main(
 //  - For each group:
 //    * Compute group-specific pairwise effect for (var1,var2), current vs. proposed.
 //    * Add linear contribution from the pairwise sufficient statistic.
-//  - Add the ratio of pseudolikelihood normalizing constants for both variables:
-//    * Call `log_ratio_pseudolikelihood_constant_variable()` separately for var1 and var2,
-//      comparing current vs. proposed pairwise weights.
+//    * Add the ratio of pseudolikelihood normalizing constants for both
+//      variables. The current-state rest scores come from the maintained
+//      residual matrices; the proposed-state rest scores adjust them by the
+//      group-specific change in the effective weight.
 //
 // Inputs:
 //  - main_effects: Matrix of main-effect parameters (fixed).
@@ -1335,8 +1228,8 @@ double log_pseudolikelihood_ratio_main(
 //  - main_effect_indices: Index ranges [row_start,row_end] for each variable.
 //  - pairwise_effect_indices: Lookup table mapping (var1,var2) → row in pairwise_effects.
 //  - projection: Group projection matrix (num_groups × (num_groups − 1)).
-//  - observations: Observation matrix (persons × variables).
-//  - group_indices: Row ranges [start,end] for each group in observations.
+//  - obs_double_groups: Per-group observation matrices converted to double.
+//  - residual_groups: Per-group rest-score matrices (persons × variables).
 //  - num_categories: Number of categories per variable.
 //  - pairwise_stats_group: Per-group pairwise sufficient statistics.
 //  - num_groups: Number of groups.
@@ -1360,8 +1253,8 @@ double log_pseudolikelihood_ratio_pairwise(
     const arma::imat& main_effect_indices,
     const arma::imat& pairwise_effect_indices,
     const arma::mat& projection,
-    const arma::imat& observations,
-    const arma::imat& group_indices,
+    const std::vector<arma::mat>& obs_double_groups,
+    const std::vector<arma::mat>& residual_groups,
     const arma::ivec& num_categories,
     const std::vector<arma::mat>& pairwise_stats_group,
     const int num_groups,
@@ -1377,7 +1270,6 @@ double log_pseudolikelihood_ratio_pairwise(
   tmp_ind(var1, var2) = 1;
   tmp_ind(var2, var1) = 1;
 
-  // Add data contribution
   for (int g = 0; g < num_groups; ++g) {
     const arma::vec proj_g = projection.row(g).t();
     const arma::mat& suff  = pairwise_stats_group[g];
@@ -1390,25 +1282,58 @@ double log_pseudolikelihood_ratio_pairwise(
       var1, var2, num_groups, proposed_pairwise_effects,
       pairwise_effect_indices, tmp_ind, proj_g
     );
+    const double delta_g = w_prop - w_cur;
 
+    // Add data contribution
     lr += 2.0 * (w_prop - w_cur) * suff(var1, var2);
-  }
 
-  // Add ratio of normalizing constant for `var1`
-  lr += log_ratio_pseudolikelihood_constant_variable(
-    main_effects, current_pairwise_effects, /* same */ main_effects,
-    proposed_pairwise_effects, main_effect_indices, pairwise_effect_indices,
-    projection, observations, group_indices, num_categories, num_groups,
-    tmp_ind, is_ordinal_variable, baseline_category, var1
-  );
+    // Add ratio of normalizing constants for both endpoint variables
+    const arma::mat& obs_g = obs_double_groups[g];
+    for (int v : {var1, var2}) {
+      const int other = (v == var1) ? var2 : var1;
+      const int num_cats = num_categories(v);
 
-  // Add ratio of normalizing constant for `var2`
-  lr += log_ratio_pseudolikelihood_constant_variable(
-    main_effects, current_pairwise_effects, /* same */ main_effects,
-    proposed_pairwise_effects, main_effect_indices, pairwise_effect_indices,
-    projection, observations, group_indices, num_categories, num_groups,
-    tmp_ind, is_ordinal_variable, baseline_category, var2
-  );
+      const arma::vec me = compute_group_main_effects(
+        v, num_groups, main_effects, main_effect_indices, proj_g
+      );
+
+      const arma::vec rest_current = residual_groups[g].col(v);
+      const arma::vec rest_proposed = rest_current + obs_g.col(other) * delta_g;
+
+      arma::vec bound_current;
+      arma::vec bound_proposed;
+      arma::vec denom_current(rest_current.n_elem, arma::fill::zeros);
+      arma::vec denom_proposed(rest_proposed.n_elem, arma::fill::zeros);
+
+      if (is_ordinal_variable(v)) {
+        bound_current = arma::clamp(rest_current * num_cats, 0.0, arma::datum::inf);
+        bound_proposed = arma::clamp(rest_proposed * num_cats, 0.0, arma::datum::inf);
+
+        denom_current += compute_denom_ordinal(
+          rest_current, me, bound_current
+        );
+        denom_proposed += compute_denom_ordinal(
+          rest_proposed, me, bound_proposed
+        );
+      } else {
+        const int ref_cat = baseline_category(v);
+        bound_current = rest_current * num_cats;
+        bound_proposed = rest_proposed * num_cats;
+
+        denom_current = compute_denom_blume_capel(
+          rest_current, me(0), me(1),
+          ref_cat, num_cats, /*Updated in place:*/bound_current
+        );
+        denom_proposed = compute_denom_blume_capel(
+          rest_proposed, me(0), me(1),
+          ref_cat, num_cats, /*Updated in place:*/bound_proposed
+        );
+      }
+
+      lr += arma::accu((bound_current - bound_proposed) +
+        ARMA_MY_LOG(denom_current) - ARMA_MY_LOG(denom_proposed));
+    }
+  }
 
   return lr;
-}
\ No newline at end of file
+}
diff --git a/src/models/bgmCompare/bgmCompare_logp_and_grad.h b/src/models/bgmCompare/bgmCompare_logp_and_grad.h
index ca00e246..e205ea9d 100644
--- a/src/models/bgmCompare/bgmCompare_logp_and_grad.h
+++ b/src/models/bgmCompare/bgmCompare_logp_and_grad.h
@@ -174,22 +174,26 @@ std::pair<double, arma::vec> logp_and_gradient(
  * Log-pseudoposterior contribution of a single main-effect parameter.
  *
  * Used by element-wise Metropolis updates. Evaluates the pseudolikelihood
- * and prior for one (variable, category/par, column h) entry.
+ * and prior for one (variable, category/par, column h) entry. Rest scores
+ * are read from the maintained residual matrices; the pairwise effects do
+ * not change during main-effect updates.
  *
+ * @param residual_groups  Per-group rest-score matrices (n_g x V)
  * @param variable   Variable index
  * @param category   Category index (ordinal variables only)
  * @param par        Parameter index: 0 = linear, 1 = quadratic (Blume-Capel only)
  * @param h          Column index: 0 = overall baseline, >0 = group difference
+ * @param normalizers_in   Optional cached per-group log-normalizer sums for
+ *                         the variable (length G); skips their computation
+ * @param normalizers_out  Optional output for the computed per-group
+ *                         log-normalizer sums (length G)
  * @see gradient() for remaining parameter descriptions
  */
 double log_pseudoposterior_main_component(
     const arma::mat& main_effects,
-    const arma::mat& pairwise_effects,
     const arma::imat& main_effect_indices,
-    const arma::imat& pairwise_effect_indices,
     const arma::mat& projection,
-    const arma::imat& observations,
-    const arma::imat& group_indices,
+    const std::vector<arma::mat>& residual_groups,
     const arma::ivec& num_categories,
     const std::vector<arma::imat>& counts_per_category_group,
     const std::vector<arma::imat>& blume_capel_stats_group,
@@ -202,7 +206,9 @@ double log_pseudoposterior_main_component(
     int par,
     int h,
     const BaseParameterPrior& difference_prior,
-    const BaseParameterPrior& threshold_prior
+    const BaseParameterPrior& threshold_prior,
+    const arma::vec* normalizers_in = nullptr,
+    arma::vec* normalizers_out = nullptr
 );
 
 /**
@@ -217,6 +223,10 @@ double log_pseudoposterior_main_component(
  * @param variable2          Second variable index
  * @param h                  Column index: 0 = overall baseline, >0 = group difference
  * @param delta              Proposed change to pairwise_effects(idx, h)
+ * @param normalizers_in     Optional cached per-group log-normalizer sums for
+ *                           the endpoint variables (G x 2); skips their computation
+ * @param normalizers_out    Optional output for the computed per-group
+ *                           log-normalizer sums (G x 2)
  * @see gradient() for remaining parameter descriptions
  */
 double log_pseudoposterior_pair_component(
@@ -238,7 +248,9 @@ double log_pseudoposterior_pair_component(
     int h,
     double delta,
     const BaseParameterPrior& interaction_prior,
-    const BaseParameterPrior& difference_prior
+    const BaseParameterPrior& difference_prior,
+    const arma::mat* normalizers_in = nullptr,
+    arma::mat* normalizers_out = nullptr
 );
 
 
@@ -247,27 +259,25 @@ double log_pseudoposterior_pair_component(
  *
  * Compares proposed vs. current main-effect parameters across all groups,
  * combining sufficient-statistic differences with normalizing-constant ratios.
+ * Both states share the rest scores held in the maintained residual matrices.
  * Used by the Metropolis-Hastings indicator update for main effects.
  *
  * @param current_main_effects   Current main-effect matrix
  * @param proposed_main_effects  Proposed main-effect matrix
+ * @param residual_groups        Per-group rest-score matrices (n_g x V)
  * @param variable               Variable whose main effect is being toggled
  * @see gradient() for remaining parameter descriptions
  */
 double log_pseudolikelihood_ratio_main(
     const arma::mat& current_main_effects,
     const arma::mat& proposed_main_effects,
-    const arma::mat& current_pairwise_effects,
     const arma::imat& main_effect_indices,
-    const arma::imat& pairwise_effect_indices,
     const arma::mat&  projection,
-    const arma::imat& observations,
-    const arma::imat& group_indices,
+    const std::vector<arma::mat>& residual_groups,
     const arma::ivec& num_categories,
     const std::vector<arma::imat>& counts_per_category_group,
     const std::vector<arma::imat>& blume_capel_stats_group,
     const int num_groups,
-    const arma::imat& inclusion_indicator,
     const arma::uvec& is_ordinal_variable,
     const arma::ivec& baseline_category,
     const int variable
@@ -278,10 +288,15 @@ double log_pseudolikelihood_ratio_main(
  *
  * Compares proposed vs. current pairwise-effect parameters for a single edge,
  * summing the data contribution and normalizing-constant ratios for both
- * endpoint variables. Used by the Metropolis-Hastings indicator update.
+ * endpoint variables. Current-state rest scores come from the maintained
+ * residual matrices; proposed-state rest scores adjust them by the change in
+ * the group-specific effective weight. Used by the Metropolis-Hastings
+ * indicator update.
  *
  * @param current_pairwise_effects   Current pairwise-effect matrix
  * @param proposed_pairwise_effects  Proposed pairwise-effect matrix
+ * @param obs_double_groups          Per-group observation matrices as double
+ * @param residual_groups            Per-group rest-score matrices (n_g x V)
  * @param var1                       First variable index
  * @param var2                       Second variable index
  * @see gradient() for remaining parameter descriptions
@@ -293,8 +308,8 @@ double log_pseudolikelihood_ratio_pairwise(
     const arma::imat& main_effect_indices,
     const arma::imat& pairwise_effect_indices,
     const arma::mat& projection,
-    const arma::imat& observations,
-    const arma::imat& group_indices,
+    const std::vector<arma::mat>& obs_double_groups,
+    const std::vector<arma::mat>& residual_groups,
     const arma::ivec& num_categories,
     const std::vector<arma::mat>& pairwise_stats_group,
     const int num_groups,
diff --git a/src/models/bgmCompare/bgmCompare_sampler.cpp b/src/models/bgmCompare/bgmCompare_sampler.cpp
index f72b5220..93af5b6f 100644
--- a/src/models/bgmCompare/bgmCompare_sampler.cpp
+++ b/src/models/bgmCompare/bgmCompare_sampler.cpp
@@ -2,6 +2,7 @@
 #include "models/bgmCompare/bgmCompare_helper.h"
 #include "models/bgmCompare/bgmCompare_logp_and_grad.h"
 #include "models/bgmCompare/bgmCompare_sampler.h"
+#include "models/bgmCompare/bgmCompare_state.h"
 #include "models/bgmCompare/bgmCompare_output.h"
 #include "mcmc/samplers/metropolis_adaptation.h"
 #include "mcmc/samplers/nuts_adaptation.h"
@@ -198,6 +199,186 @@ void impute_missing_bgmcompare(
 
 
 
+// Performs one cached random-walk Metropolis update of a single main-effect
+// parameter (row, h) of `variable`.
+//
+// The current-state log-pseudoposterior is reconstructed from the sweep
+// state's per-variable normalizer cache, computing and caching the
+// normalizer sums on a miss; only the proposed state is evaluated in full.
+// On acceptance the parameter and the variable's cached normalizer sums are
+// updated in place. Shared by the Metropolis sweep and the proposal-sd tuner.
+static StepResult metropolis_update_main_effect_cached(
+    arma::mat& main_effects,
+    const arma::imat& main_effect_indices,
+    const arma::imat& inclusion_indicator,
+    const arma::mat& projection,
+    const arma::ivec& num_categories,
+    CompareSweepState& sweep_state,
+    const int num_groups,
+    const std::vector<arma::imat>& counts_per_category,
+    const std::vector<arma::imat>& blume_capel_stats,
+    const arma::uvec& is_ordinal_variable,
+    const arma::ivec& baseline_category,
+    const int row,
+    const int variable,
+    const int category,
+    const int par,
+    const int h,
+    const double proposal_sd,
+    SafeRNG& rng,
+    const BaseParameterPrior& difference_prior,
+    const BaseParameterPrior& threshold_prior
+) {
+  const double current = main_effects(row, h);
+
+  auto component = [&](const arma::vec* norm_in, arma::vec* norm_out) {
+    return log_pseudoposterior_main_component(
+      main_effects, main_effect_indices, projection,
+      sweep_state.residual, num_categories, counts_per_category,
+      blume_capel_stats, num_groups, inclusion_indicator,
+      is_ordinal_variable, baseline_category,
+      variable, category, par, h,
+      difference_prior, threshold_prior,
+      norm_in, norm_out
+    );
+  };
+
+  // Current-state value: reconstruct from the cached per-group normalizer
+  // sums, or compute and cache them on a miss.
+  double logp_current;
+  if (sweep_state.normalizer_valid(variable)) {
+    const arma::vec norm_cur = sweep_state.log_normalizer.row(variable).t();
+    logp_current = component(&norm_cur, nullptr);
+  } else {
+    arma::vec norm_cur(num_groups);
+    logp_current = component(nullptr, &norm_cur);
+    sweep_state.log_normalizer.row(variable) = norm_cur.t();
+    sweep_state.normalizer_valid(variable) = 1;
+  }
+
+  arma::vec norm_prop(num_groups);
+  auto log_post_proposed = [&](double theta) {
+    main_effects(row, h) = theta;
+    return component(nullptr, &norm_prop);
+  };
+
+  StepResult result = metropolis_step_cached(
+    current, proposal_sd, logp_current, log_post_proposed, rng
+  );
+  main_effects(row, h) = result.state[0];
+  if (result.state[0] != current) {
+    // Accepted: the proposed-state normalizer sums become current.
+    sweep_state.log_normalizer.row(variable) = norm_prop.t();
+  }
+
+  return result;
+}
+
+
+
+// Performs one cached random-walk Metropolis update of a single
+// pairwise-effect parameter (idx, h) of the pair (var1, var2).
+//
+// The current-state log-pseudoposterior is reconstructed from the sweep
+// state's normalizer cache of both endpoint variables, computing and caching
+// the normalizer sums on a miss; only the proposed state is evaluated in
+// full. On acceptance the parameter, the effective weights, the residual
+// columns, and both variables' cached normalizer sums are updated in place.
+// Shared by the Metropolis sweep and the proposal-sd tuner.
+static StepResult metropolis_update_pairwise_effect_cached(
+    arma::mat& main_effects,
+    arma::mat& pairwise_effects,
+    const arma::imat& main_effect_indices,
+    const arma::imat& pairwise_effect_indices,
+    const arma::imat& inclusion_indicator,
+    const arma::mat& projection,
+    const arma::ivec& num_categories,
+    CompareSweepState& sweep_state,
+    const int num_groups,
+    const std::vector<arma::mat>& pairwise_stats,
+    const arma::uvec& is_ordinal_variable,
+    const arma::ivec& baseline_category,
+    const int var1,
+    const int var2,
+    const int h,
+    const double proposal_sd,
+    SafeRNG& rng,
+    const BaseParameterPrior& interaction_prior,
+    const BaseParameterPrior& difference_prior
+) {
+  const int idx = pairwise_effect_indices(var1, var2);
+  const double current = pairwise_effects(idx, h);
+
+  auto component = [&](double delta, const arma::mat* norm_in,
+                       arma::mat* norm_out) {
+    return log_pseudoposterior_pair_component(
+      main_effects, pairwise_effects, main_effect_indices,
+      pairwise_effect_indices, projection, sweep_state.obs_double,
+      num_categories, pairwise_stats, sweep_state.residual, num_groups,
+      inclusion_indicator, is_ordinal_variable, baseline_category,
+      var1, var2, h, delta,
+      interaction_prior, difference_prior,
+      norm_in, norm_out
+    );
+  };
+
+  // Current-state value: reconstruct from the cached per-group normalizer
+  // sums of both endpoints, or compute and cache them on a miss.
+  double logp_current;
+  arma::mat norm_cur(num_groups, 2);
+  if (sweep_state.normalizer_valid(var1) &&
+      sweep_state.normalizer_valid(var2)) {
+    norm_cur.col(0) = sweep_state.log_normalizer.row(var1).t();
+    norm_cur.col(1) = sweep_state.log_normalizer.row(var2).t();
+    logp_current = component(0.0, &norm_cur, nullptr);
+  } else {
+    logp_current = component(0.0, nullptr, &norm_cur);
+    sweep_state.log_normalizer.row(var1) = norm_cur.col(0).t();
+    sweep_state.log_normalizer.row(var2) = norm_cur.col(1).t();
+    sweep_state.normalizer_valid(var1) = 1;
+    sweep_state.normalizer_valid(var2) = 1;
+  }
+
+  arma::mat norm_prop(num_groups, 2);
+  auto log_post_proposed = [&](double theta) {
+    return component(theta - current, nullptr, &norm_prop);
+  };
+
+  StepResult result = metropolis_step_cached(
+    current, proposal_sd, logp_current, log_post_proposed, rng
+  );
+  const double value = result.state[0];
+
+  // Update the parameter and sweep state if the move was accepted
+  if (current != value) {
+    const double delta = value - current;
+    pairwise_effects(idx, h) = value;
+
+    for (int g = 0; g < num_groups; ++g) {
+      const arma::vec proj_g = projection.row(g).t();
+      double delta_g = (h == 0) ? delta : delta * proj_g(h - 1);
+
+      // Update pairwise_group for this group
+      sweep_state.pairwise_group[g](var1, var2) += delta_g;
+      sweep_state.pairwise_group[g](var2, var1) += delta_g;
+
+      // Update residual matrix columns
+      sweep_state.residual[g].col(var1) +=
+        sweep_state.obs_double[g].col(var2) * delta_g;
+      sweep_state.residual[g].col(var2) +=
+        sweep_state.obs_double[g].col(var1) * delta_g;
+    }
+
+    // The proposed-state normalizer sums become current.
+    sweep_state.log_normalizer.row(var1) = norm_prop.col(0).t();
+    sweep_state.log_normalizer.row(var2) = norm_prop.col(1).t();
+  }
+
+  return result;
+}
+
+
+
 // Updates main effect parameters in bgmCompare using a random-walk Metropolis step.
 //
 // For each variable, the function proposes new parameter values for either:
@@ -222,15 +403,14 @@ void impute_missing_bgmcompare(
 // Inputs:
 //  - main_effects: Matrix of main effect parameters [rows = effects, cols = groups];
 //                  updated in place.
-//  - pairwise_effects: Current pairwise effects (passed through to log posterior).
 //  - main_effect_indices: Row index ranges for each variable’s main effects.
-//  - pairwise_effect_indices: Index map for pairwise effects.
 //  - inclusion_indicator: Indicator matrix; diagonal entries control group differences.
 //  - projection: Group projection matrix.
 //  - num_categories: Number of categories for each variable.
-//  - observations: Data matrix [persons × variables].
+//  - sweep_state: Maintained per-group sweep state; the residual matrices
+//                 supply the rest scores, and accepted moves update the
+//                 per-variable normalizer cache.
 //  - num_groups: Number of groups (G).
-//  - group_indices: Row ranges per group in `observations`.
 //  - counts_per_category, blume_capel_stats: Group-specific sufficient statistics.
 //  - is_ordinal_variable: Indicator for ordinal vs. Blume–Capel.
 //  - baseline_category: Reference categories (Blume–Capel only).
@@ -244,21 +424,18 @@ void impute_missing_bgmcompare(
 //
 // Notes:
 //  - Acceptance probabilities are stored per parameter and fed to `metropolis_adapt.update()`.
-//  - This function does not alter pairwise effects, but passes them into
-//    the posterior for likelihood consistency.
+//  - This function does not alter pairwise effects; the sweep state stays
+//    valid throughout the sweep.
 //  - The helper lambda `do_update` encapsulates the proposal/accept/revert loop
 //    for a single parameter, improving readability.
 void update_main_effects_metropolis_bgmcompare (
     arma::mat& main_effects,
-    arma::mat& pairwise_effects,
     const arma::imat& main_effect_indices,
-    const arma::imat& pairwise_effect_indices,
     const arma::imat& inclusion_indicator,
     const arma::mat& projection,
     const arma::ivec& num_categories,
-    const arma::imat& observations,
+    CompareSweepState& sweep_state,
     const int num_groups,
-    const arma::imat& group_indices,
     const std::vector<arma::imat>& counts_per_category,
     const std::vector<arma::imat>& blume_capel_stats,
     const arma::uvec& is_ordinal_variable,
@@ -270,7 +447,7 @@ void update_main_effects_metropolis_bgmcompare (
     const BaseParameterPrior& difference_prior,
     const BaseParameterPrior& threshold_prior
 ) {
-  const int num_vars = observations.n_cols;
+  const int num_vars = inclusion_indicator.n_rows;
   arma::umat index_mask_main = arma::zeros<arma::umat>(proposal_sd_main.n_rows,
                                                        proposal_sd_main.n_cols);
   arma::mat accept_prob_main = arma::zeros<arma::mat>(proposal_sd_main.n_rows,
@@ -279,24 +456,13 @@ void update_main_effects_metropolis_bgmcompare (
   // --- helper for one update ---
   auto do_update = [&](int row, int variable, int category, int par, int h) {
     index_mask_main(row, h) = 1;
-    double& current = main_effects(row, h);
-    double proposal_sd = proposal_sd_main(row, h);
-
-    auto log_post = [&](double theta) {
-      main_effects(row, h) = theta;
-      return log_pseudoposterior_main_component(
-        main_effects, pairwise_effects, main_effect_indices,
-        pairwise_effect_indices, projection, observations,
-        group_indices, num_categories, counts_per_category,
-        blume_capel_stats, num_groups, inclusion_indicator,
-        is_ordinal_variable, baseline_category,
-        variable, category, par, h,
-        difference_prior, threshold_prior
-      );
-    };
-
-    StepResult result = metropolis_step(current, proposal_sd, log_post, rng);
-    current = result.state[0];
+    StepResult result = metropolis_update_main_effect_cached(
+      main_effects, main_effect_indices, inclusion_indicator, projection,
+      num_categories, sweep_state, num_groups, counts_per_category,
+      blume_capel_stats, is_ordinal_variable, baseline_category,
+      row, variable, category, par, h, proposal_sd_main(row, h), rng,
+      difference_prior, threshold_prior
+    );
     accept_prob_main(row, h) = result.accept_prob;
   };
 
@@ -361,9 +527,9 @@ void update_main_effects_metropolis_bgmcompare (
 //  - inclusion_indicator: Indicator matrix; off-diagonal entries control group differences.
 //  - projection: Group projection matrix.
 //  - num_categories: Number of categories per variable.
-//  - observations: Data matrix [persons × variables].
+//  - sweep_state: Maintained per-group sweep state; accepted moves update
+//                 the effective weights and residual matrices in place.
 //  - num_groups: Number of groups (G).
-//  - group_indices: Row ranges per group in `observations`.
 //  - pairwise_stats: Group-specific sufficient statistics for pairwise effects.
 //  - is_ordinal_variable: Indicator for ordinal vs. Blume–Capel variables.
 //  - baseline_category: Reference categories (Blume–Capel only).
@@ -388,9 +554,8 @@ void update_pairwise_effects_metropolis_bgmcompare (
     const arma::imat& inclusion_indicator,
     const arma::mat& projection,
     const arma::ivec& num_categories,
-    const arma::imat& observations,
+    CompareSweepState& sweep_state,
     const int num_groups,
-    const arma::imat& group_indices,
     const std::vector<arma::mat>& pairwise_stats,
     const arma::uvec& is_ordinal_variable,
     const arma::ivec& baseline_category,
@@ -401,81 +566,23 @@ void update_pairwise_effects_metropolis_bgmcompare (
     const BaseParameterPrior& interaction_prior,
     const BaseParameterPrior& difference_prior
 ) {
-  int num_variables = observations.n_cols;
+  int num_variables = inclusion_indicator.n_rows;
   int num_pairs = num_variables * (num_variables - 1) / 2;
   arma::mat accept_prob_pair = arma::zeros<arma::mat>(num_pairs, num_groups);
   arma::umat index_mask_pair = arma::zeros<arma::umat>(num_pairs, num_groups);
 
-  // --- Build group-specific pairwise matrices and residual matrices ---
-  std::vector<arma::mat> pairwise_groups(num_groups);
-  std::vector<arma::mat> residual_matrices(num_groups);
-  std::vector<arma::mat> obs_double_groups(num_groups);
-
-  for (int g = 0; g < num_groups; ++g) {
-    const int r0 = group_indices(g, 0);
-    const int r1 = group_indices(g, 1);
-    const arma::vec proj_g = projection.row(g).t();
-
-    // Build group-specific pairwise effects matrix
-    arma::mat pairwise_group(num_variables, num_variables, arma::fill::zeros);
-    for (int v = 0; v < num_variables - 1; ++v) {
-      for (int u = v + 1; u < num_variables; ++u) {
-        double w = compute_group_pairwise_effects(
-          v, u, num_groups, pairwise_effects, pairwise_effect_indices,
-          inclusion_indicator, proj_g
-        );
-        pairwise_group(v, u) = w;
-        pairwise_group(u, v) = w;
-      }
-    }
-    pairwise_groups[g] = pairwise_group;
-
-    // Convert observations to double and compute residual matrix
-    obs_double_groups[g] = arma::conv_to<arma::mat>::from(observations.rows(r0, r1));
-    residual_matrices[g] = obs_double_groups[g] * pairwise_group;
-  }
-
-  // --- helper for one update using optimized residual-based function ---
+  // --- helper for one update using the cached pairwise Metropolis step ---
   auto do_update = [&](int var1, int var2, int h) {
     int idx = pairwise_effect_indices(var1, var2);
     index_mask_pair(idx, h) = 1;
-    double current = pairwise_effects(idx, h);
-    double proposal_sd = proposal_sd_pair(idx, h);
-
-    auto log_post = [&](double theta) {
-      double delta = theta - current;
-      return log_pseudoposterior_pair_component(
-        main_effects, pairwise_effects, main_effect_indices,
-        pairwise_effect_indices, projection, obs_double_groups,
-        num_categories, pairwise_stats, residual_matrices, num_groups,
-        inclusion_indicator, is_ordinal_variable, baseline_category,
-        var1, var2, h, delta,
-        interaction_prior, difference_prior
-      );
-    };
-
-    StepResult result = metropolis_step(current, proposal_sd, log_post, rng);
-    double value = result.state[0];
-
-    // Update residual matrices if move was accepted
-    if (current != value) {
-      double delta = value - current;
-      pairwise_effects(idx, h) = value;
-
-      for (int g = 0; g < num_groups; ++g) {
-        const arma::vec proj_g = projection.row(g).t();
-        double delta_g = (h == 0) ? delta : delta * proj_g(h - 1);
-
-        // Update pairwise_group for this group
-        pairwise_groups[g](var1, var2) += delta_g;
-        pairwise_groups[g](var2, var1) += delta_g;
-
-        // Update residual matrix columns
-        residual_matrices[g].col(var1) += obs_double_groups[g].col(var2) * delta_g;
-        residual_matrices[g].col(var2) += obs_double_groups[g].col(var1) * delta_g;
-      }
-    }
-
+    StepResult result = metropolis_update_pairwise_effect_cached(
+      main_effects, pairwise_effects, main_effect_indices,
+      pairwise_effect_indices, inclusion_indicator, projection,
+      num_categories, sweep_state, num_groups, pairwise_stats,
+      is_ordinal_variable, baseline_category,
+      var1, var2, h, proposal_sd_pair(idx, h), rng,
+      interaction_prior, difference_prior
+    );
     accept_prob_pair(idx, h) = result.accept_prob;
   };
 
@@ -549,6 +656,7 @@ double find_initial_stepsize_bgmcompare(
     const arma::mat& projection,
     const arma::ivec& num_categories,
     const arma::imat& observations,
+    const arma::mat& obs_double,
     const int num_groups,
     const arma::imat& group_indices,
     const std::vector<arma::imat>& counts_per_category,
@@ -569,9 +677,6 @@ double find_initial_stepsize_bgmcompare(
   arma::mat current_main = main_effects;
   arma::mat current_pair = pairwise_effects;
 
-  // Pre-convert observations to double once (avoids repeated conversion in gradient evaluations)
-  const arma::mat obs_double = arma::conv_to<arma::mat>::from(observations);
-
   auto index_maps = build_index_maps(
     main_effects, pairwise_effects,
     inclusion_indicator,
@@ -647,6 +752,7 @@ double find_initial_stepsize_bgmcompare(
 //  - projection: Group projection matrix for contrasts.
 //  - num_categories: Number of categories per variable [V].
 //  - observations: Data matrix [N × V].
+//  - obs_double: Data matrix pre-converted to double [N × V].
 //  - num_groups: Number of groups.
 //  - group_indices: Row ranges for each group in `observations`.
 //  - counts_per_category, blume_capel_stats: Per-group sufficient statistics.
@@ -683,6 +789,7 @@ StepResult update_nuts_bgmcompare(
     const arma::mat& projection,
     const arma::ivec& num_categories,
     const arma::imat& observations,
+    const arma::mat& obs_double,
     const int num_groups,
     const arma::imat& group_indices,
     const std::vector<arma::imat>& counts_per_category,
@@ -700,9 +807,6 @@ StepResult update_nuts_bgmcompare(
     const BaseParameterPrior& difference_prior,
     const BaseParameterPrior& threshold_prior
 ) {
-  // Pre-convert observations to double once (avoids repeated conversion in gradient evaluations)
-  const arma::mat obs_double = arma::conv_to<arma::mat>::from(observations);
-
   arma::vec current_state = vectorize_model_parameters_bgmcompare(
     main_effects, pairwise_effects, inclusion_indicator,
     main_effect_indices, pairwise_effect_indices, num_categories,
@@ -833,9 +937,9 @@ StepResult update_nuts_bgmcompare(
 //  - inclusion_indicator: Marks which group differences are active.
 //  - projection: Group projection matrix.
 //  - num_categories: Categories per variable.
-//  - observations: Data matrix [N × V].
+//  - sweep_state: Maintained per-group sweep state; accepted pairwise moves
+//                 update the effective weights and residual matrices in place.
 //  - num_groups: Number of groups.
-//  - group_indices: Row ranges per group in `observations`.
 //  - counts_per_category, blume_capel_stats: Per-group sufficient statistics for main effects.
 //  - pairwise_stats: Per-group sufficient statistics for pairwise effects.
 //  - is_ordinal_variable: Marks ordinal vs. Blume–Capel variables.
@@ -866,9 +970,8 @@ void tune_proposal_sd_bgmcompare(
     const arma::imat& inclusion_indicator,
     const arma::mat& projection,
     const arma::ivec& num_categories,
-    const arma::imat& observations,
+    CompareSweepState& sweep_state,
     int num_groups,
-    const arma::imat& group_indices,
     const std::vector<arma::imat>& counts_per_category,
     const std::vector<arma::imat>& blume_capel_stats,
     const std::vector<arma::mat>& pairwise_stats,
@@ -889,7 +992,7 @@ void tune_proposal_sd_bgmcompare(
   double t = iteration - sched.stage3b_start + 1;
   double rm_weight = std::pow(t, -rm_decay);
 
-  const int V = observations.n_cols;
+  const int V = inclusion_indicator.n_rows;
 
   // --- MAIN EFFECTS ---
   for (int var = 0; var < V; ++var) {
@@ -904,25 +1007,16 @@ void tune_proposal_sd_bgmcompare(
       for (int c = 0; c < ncat; ++c) {
         int row = start + c;
         for (int h = 0; h < hmax; ++h) {
-          double& current = main_effects(row, h);
           double& prop_sd = proposal_sd_main_effects(row, h);
 
-          auto log_post = [&](double theta) {
-            main_effects(row, h) = theta;
-            return log_pseudoposterior_main_component(
-              main_effects, pairwise_effects,
-              main_effect_indices, pairwise_effect_indices,
-              projection, observations, group_indices,
-              num_categories, counts_per_category, blume_capel_stats,
-              num_groups, inclusion_indicator, is_ordinal_variable,
-              baseline_category,
-              var, c, -1, h,
-              difference_prior, threshold_prior
-            );
-          };
-
-          StepResult result = metropolis_step(current, prop_sd, log_post, rng);
-          current = result.state[0];
+          StepResult result = metropolis_update_main_effect_cached(
+            main_effects, main_effect_indices, inclusion_indicator,
+            projection, num_categories, sweep_state, num_groups,
+            counts_per_category, blume_capel_stats, is_ordinal_variable,
+            baseline_category,
+            row, var, c, -1, h, prop_sd, rng,
+            difference_prior, threshold_prior
+          );
           prop_sd = update_proposal_sd_with_robbins_monro(
             prop_sd, MY_LOG(result.accept_prob), rm_weight, target_accept
           );
@@ -933,25 +1027,16 @@ void tune_proposal_sd_bgmcompare(
       for (int par = 0; par < 2; ++par) {
         int row = start + par;
         for (int h = 0; h < hmax; ++h) {
-          double& current = main_effects(row, h);
           double& prop_sd = proposal_sd_main_effects(row, h);
 
-          auto log_post = [&](double theta) {
-            main_effects(row, h) = theta;
-            return log_pseudoposterior_main_component(
-              main_effects, pairwise_effects,
-              main_effect_indices, pairwise_effect_indices,
-              projection, observations, group_indices,
-              num_categories, counts_per_category, blume_capel_stats,
-              num_groups, inclusion_indicator, is_ordinal_variable,
-              baseline_category,
-              var, -1, par, h,
-              difference_prior, threshold_prior
-            );
-          };
-
-          StepResult result = metropolis_step(current, prop_sd, log_post, rng);
-          current = result.state[0];
+          StepResult result = metropolis_update_main_effect_cached(
+            main_effects, main_effect_indices, inclusion_indicator,
+            projection, num_categories, sweep_state, num_groups,
+            counts_per_category, blume_capel_stats, is_ordinal_variable,
+            baseline_category,
+            row, var, -1, par, h, prop_sd, rng,
+            difference_prior, threshold_prior
+          );
           prop_sd = update_proposal_sd_with_robbins_monro(
             prop_sd, MY_LOG(result.accept_prob), rm_weight, target_accept
           );
@@ -961,32 +1046,6 @@ void tune_proposal_sd_bgmcompare(
   }
 
   // --- PAIRWISE EFFECTS ---
-  // Build group-specific pairwise matrices and residual matrices
-  std::vector<arma::mat> pairwise_groups(num_groups);
-  std::vector<arma::mat> residual_matrices(num_groups);
-  std::vector<arma::mat> obs_double_groups(num_groups);
-
-  for (int g = 0; g < num_groups; ++g) {
-    const int r0 = group_indices(g, 0);
-    const int r1 = group_indices(g, 1);
-    const arma::vec proj_g = projection.row(g).t();
-
-    arma::mat pairwise_group(V, V, arma::fill::zeros);
-    for (int v = 0; v < V - 1; ++v) {
-      for (int u = v + 1; u < V; ++u) {
-        double w = compute_group_pairwise_effects(
-          v, u, num_groups, pairwise_effects, pairwise_effect_indices,
-          inclusion_indicator, proj_g
-        );
-        pairwise_group(v, u) = w;
-        pairwise_group(u, v) = w;
-      }
-    }
-    pairwise_groups[g] = pairwise_group;
-    obs_double_groups[g] = arma::conv_to<arma::mat>::from(observations.rows(r0, r1));
-    residual_matrices[g] = obs_double_groups[g] * pairwise_group;
-  }
-
   for (int v1 = 0; v1 < V - 1; ++v1) {
     for (int v2 = v1 + 1; v2 < V; ++v2) {
       int idx = pairwise_effect_indices(v1, v2);
@@ -994,40 +1053,16 @@ void tune_proposal_sd_bgmcompare(
       int hmax = group_differences ? num_groups : 1;
 
       for (int h = 0; h < hmax; ++h) {
-        double& current = pairwise_effects(idx, h);
         double& prop_sd = proposal_sd_pairwise_effects(idx, h);
 
-        auto log_post = [&](double theta) {
-          double delta = theta - current;
-          return log_pseudoposterior_pair_component(
-            main_effects, pairwise_effects,
-            main_effect_indices, pairwise_effect_indices,
-            projection, obs_double_groups,
-            num_categories, pairwise_stats, residual_matrices, num_groups,
-            inclusion_indicator, is_ordinal_variable, baseline_category,
-            v1, v2, h, delta,
-            interaction_prior, difference_prior
-          );
-        };
-
-        StepResult result = metropolis_step(current, prop_sd, log_post, rng);
-        double value = result.state[0];
-
-        if (current != value) {
-          double delta = value - current;
-          current = value;
-
-          for (int g = 0; g < num_groups; ++g) {
-            const arma::vec proj_g = projection.row(g).t();
-            double delta_g = (h == 0) ? delta : delta * proj_g(h - 1);
-
-            pairwise_groups[g](v1, v2) += delta_g;
-            pairwise_groups[g](v2, v1) += delta_g;
-
-            residual_matrices[g].col(v1) += obs_double_groups[g].col(v2) * delta_g;
-            residual_matrices[g].col(v2) += obs_double_groups[g].col(v1) * delta_g;
-          }
-        }
+        StepResult result = metropolis_update_pairwise_effect_cached(
+          main_effects, pairwise_effects, main_effect_indices,
+          pairwise_effect_indices, inclusion_indicator, projection,
+          num_categories, sweep_state, num_groups, pairwise_stats,
+          is_ordinal_variable, baseline_category,
+          v1, v2, h, prop_sd, rng,
+          interaction_prior, difference_prior
+        );
 
         prop_sd = update_proposal_sd_with_robbins_monro(
           prop_sd, MY_LOG(result.accept_prob), rm_weight, target_accept
@@ -1061,9 +1096,9 @@ void tune_proposal_sd_bgmcompare(
 //  - main_effects, pairwise_effects: Parameter matrices, updated in place.
 //  - main_effect_indices, pairwise_effect_indices: Index maps for parameters.
 //  - projection: Group projection matrix.
-//  - observations: Data matrix [N × V].
+//  - sweep_state: Maintained per-group sweep state; accepted pairwise flips
+//    update the effective weights and residual matrices in place.
 //  - num_groups: Number of groups.
-//  - group_indices: Row ranges per group in `observations`.
 //  - num_categories: Categories per variable [V × G].
 //  - inclusion_indicator: Indicator matrix for differences, updated in place.
 //  - is_ordinal_variable: Marks ordinal vs. Blume–Capel variables [V].
@@ -1094,9 +1129,8 @@ void update_indicator_differences_metropolis_bgmcompare (
     const arma::imat& main_effect_indices,
     const arma::imat& pairwise_effect_indices,
     const arma::mat& projection,
-    const arma::imat& observations,
+    CompareSweepState& sweep_state,
     const int num_groups,
-    const arma::imat& group_indices,
     const arma::imat& num_categories,
     arma::imat& inclusion_indicator,
     const arma::uvec& is_ordinal_variable,
@@ -1110,7 +1144,7 @@ void update_indicator_differences_metropolis_bgmcompare (
     SafeRNG& rng,
     const BaseParameterPrior& difference_prior
 ) {
-  const int num_variables = observations.n_cols;
+  const int num_variables = inclusion_indicator.n_rows;
 
   // --- main effects ---
   // Skip main effect indicator updates if main_difference_selection is disabled
@@ -1140,10 +1174,10 @@ void update_indicator_differences_metropolis_bgmcompare (
 
     // Calculate log acceptance probability
     double log_accept = log_pseudolikelihood_ratio_main(
-      current_main_effects, proposed_main_effects, pairwise_effects,
-      main_effect_indices, pairwise_effect_indices, projection,
-      observations, group_indices, num_categories, counts_per_category,
-      blume_capel_stats, num_groups, inclusion_indicator,
+      current_main_effects, proposed_main_effects,
+      main_effect_indices, projection, sweep_state.residual,
+      num_categories, counts_per_category,
+      blume_capel_stats, num_groups,
       is_ordinal_variable, baseline_category, var
     );
 
@@ -1185,6 +1219,10 @@ void update_indicator_differences_metropolis_bgmcompare (
       inclusion_indicator(var, var) = proposed_ind;
       main_effects.rows(start, stop).cols(1, num_groups - 1) =
         proposed_main_effects.rows(start, stop).cols(1, num_groups - 1);
+
+      // The variable's main effects changed; its cached normalizer sums
+      // are stale.
+      sweep_state.normalizer_valid(var) = 0;
     }
   }
   } // end if (main_difference_selection)
@@ -1216,8 +1254,9 @@ void update_indicator_differences_metropolis_bgmcompare (
     // Calculate log acceptance probability
     double log_accept = log_pseudolikelihood_ratio_pairwise(
       main_effects, current_pairwise_effects, proposed_pairwise_effects,
-      main_effect_indices, pairwise_effect_indices, projection, observations,
-      group_indices, num_categories, pairwise_stats, num_groups,
+      main_effect_indices, pairwise_effect_indices, projection,
+      sweep_state.obs_double, sweep_state.residual,
+      num_categories, pairwise_stats, num_groups,
       inclusion_indicator, is_ordinal_variable, baseline_category, var1, var2
     );
 
@@ -1269,6 +1308,31 @@ void update_indicator_differences_metropolis_bgmcompare (
       for (int h = 1; h < num_groups; h++) {
         pairwise_effects(int_index, h) = proposed_pairwise_effects(int_index, h);
       }
+
+      // Maintain the sweep state: the flip changes the pair's effective
+      // weight per group, which shifts two residual columns per group.
+      for (int g = 0; g < num_groups; g++) {
+        const arma::vec proj_g = projection.row(g).t();
+        const double w_old = sweep_state.pairwise_group[g](var1, var2);
+        const double w_new = compute_group_pairwise_effects(
+          var1, var2, num_groups, pairwise_effects, pairwise_effect_indices,
+          inclusion_indicator, proj_g
+        );
+        const double delta_g = w_new - w_old;
+
+        sweep_state.pairwise_group[g](var1, var2) = w_new;
+        sweep_state.pairwise_group[g](var2, var1) = w_new;
+
+        sweep_state.residual[g].col(var1) +=
+          sweep_state.obs_double[g].col(var2) * delta_g;
+        sweep_state.residual[g].col(var2) +=
+          sweep_state.obs_double[g].col(var1) * delta_g;
+      }
+
+      // Both endpoints' rest scores changed; their cached normalizer sums
+      // are stale.
+      sweep_state.normalizer_valid(var1) = 0;
+      sweep_state.normalizer_valid(var2) = 0;
     }
   }
 }
@@ -1296,6 +1360,8 @@ void update_indicator_differences_metropolis_bgmcompare (
 //
 // Inputs:
 //  - observations: Data matrix [N × V].
+//  - sweep_state: Maintained per-group sweep state (double observations,
+//    effective weights, residual matrices), kept consistent across updates.
 //  - num_categories: Number of categories per variable [V].
 //  - pairwise_scale, difference_scale: Prior scale parameters.
 //  - counts_per_category, blume_capel_stats: Sufficient statistics per group.
@@ -1332,6 +1398,7 @@ void update_indicator_differences_metropolis_bgmcompare (
 //    global (NUTS).
 void gibbs_update_step_bgmcompare (
     const arma::imat& observations,
+    CompareSweepState& sweep_state,
     const arma::ivec& num_categories,
     const std::vector<arma::imat>& counts_per_category,
     const std::vector<arma::imat>& blume_capel_stats,
@@ -1375,14 +1442,20 @@ void gibbs_update_step_bgmcompare (
       inclusion_indicator, main_effects, pairwise_effects, main_effect_indices,
       pairwise_effect_indices, inclusion_probability, main_difference_selection, rng
     );
+    // The excluded pairs' difference columns were zeroed; refresh the
+    // effective weights and residual matrices.
+    rebuild_sweep_state_weights(
+      sweep_state, pairwise_effects, pairwise_effect_indices,
+      inclusion_indicator, projection, num_groups
+    );
   }
 
   // Step 1: Difference selection via MH indicator updates (if enabled)
   if (schedule.selection_enabled(iteration)) {
     update_indicator_differences_metropolis_bgmcompare (
         inclusion_probability, index, main_effects, pairwise_effects,
-        main_effect_indices, pairwise_effect_indices, projection, observations,
-        num_groups, group_indices, num_categories, inclusion_indicator,
+        main_effect_indices, pairwise_effect_indices, projection, sweep_state,
+        num_groups, num_categories, inclusion_indicator,
         is_ordinal_variable, baseline_category, proposal_sd_main,
         proposal_sd_pair, counts_per_category,
         blume_capel_stats, pairwise_stats, main_difference_selection, rng,
@@ -1393,9 +1466,9 @@ void gibbs_update_step_bgmcompare (
   // Step 2: Update parameters
   if(update_method == adaptive_metropolis) {
     update_main_effects_metropolis_bgmcompare (
-        main_effects, pairwise_effects, main_effect_indices,
-        pairwise_effect_indices, inclusion_indicator, projection,
-        num_categories, observations, num_groups, group_indices,
+        main_effects, main_effect_indices,
+        inclusion_indicator, projection,
+        num_categories, sweep_state, num_groups,
         counts_per_category, blume_capel_stats, is_ordinal_variable,
         baseline_category, iteration,
         metropolis_adapt_main, rng, proposal_sd_main,
@@ -1405,7 +1478,7 @@ void gibbs_update_step_bgmcompare (
     update_pairwise_effects_metropolis_bgmcompare (
         main_effects, pairwise_effects, main_effect_indices,
         pairwise_effect_indices, inclusion_indicator, projection,
-        num_categories, observations, num_groups, group_indices,
+        num_categories, sweep_state, num_groups,
         pairwise_stats, is_ordinal_variable, baseline_category,
         iteration, metropolis_adapt_pair, rng,
         proposal_sd_pair,
@@ -1415,7 +1488,8 @@ void gibbs_update_step_bgmcompare (
     StepResult result = update_nuts_bgmcompare(
       main_effects, pairwise_effects, main_effect_indices,
       pairwise_effect_indices, inclusion_indicator, projection, num_categories,
-      observations, num_groups, group_indices, counts_per_category,
+      observations, sweep_state.obs_double_all, num_groups, group_indices,
+      counts_per_category,
       blume_capel_stats, pairwise_stats, is_ordinal_variable,
       baseline_category,
       nuts_max_depth, iteration, nuts_adapt, learn_mass_matrix,
@@ -1423,6 +1497,13 @@ void gibbs_update_step_bgmcompare (
       interaction_prior, difference_prior, threshold_prior
     );
 
+    // The NUTS update changes the parameters wholesale; rebuild the
+    // effective weights and residual matrices.
+    rebuild_sweep_state_weights(
+      sweep_state, pairwise_effects, pairwise_effect_indices,
+      inclusion_indicator, projection, num_groups
+    );
+
     if (iteration >= schedule.total_warmup) {
       int sample_index = iteration - schedule.total_warmup;
       if (auto diag = std::dynamic_pointer_cast<NUTSDiagnostics>(result.diagnostics)) {
@@ -1438,8 +1519,8 @@ void gibbs_update_step_bgmcompare (
   tune_proposal_sd_bgmcompare(
     proposal_sd_main, proposal_sd_pair, main_effects,
     pairwise_effects, main_effect_indices, pairwise_effect_indices,
-    inclusion_indicator, projection, num_categories, observations, num_groups,
-    group_indices, counts_per_category, blume_capel_stats,
+    inclusion_indicator, projection, num_categories, sweep_state, num_groups,
+    counts_per_category, blume_capel_stats,
     pairwise_stats, is_ordinal_variable, baseline_category,
     iteration, rng, schedule,
     interaction_prior, difference_prior, threshold_prior
@@ -1599,13 +1680,25 @@ bgmCompareOutput run_gibbs_sampler_bgmCompare(
   arma::mat proposal_sd_main(num_main, num_groups, arma::fill::ones);
   arma::mat proposal_sd_pair(num_pair, num_groups, arma::fill::ones);
 
+  // --- Persistent per-group sweep state (double observations, effective
+  //     pairwise weights, residual matrices), maintained across iterations.
+  CompareSweepState sweep_state;
+  initialize_sweep_state_observations(
+    sweep_state, observations, group_indices, num_groups
+  );
+  rebuild_sweep_state_weights(
+    sweep_state, pairwise_effects, pairwise_effect_indices,
+    inclusion_indicator, projection, num_groups
+  );
+
   // --- Optional NUTS warmup stage
   double initial_step_size = 1.0;
   if (update_method == nuts) {
     initial_step_size = find_initial_stepsize_bgmcompare(
       main_effects, pairwise_effects, main_effect_indices,
       pairwise_effect_indices, inclusion_indicator, projection, num_categories,
-      observations, num_groups, group_indices, counts_per_category,
+      observations, sweep_state.obs_double_all, num_groups, group_indices,
+      counts_per_category,
       blume_capel_stats, pairwise_stats, is_ordinal_variable,
       baseline_category,
       target_accept, rng,
@@ -1656,11 +1749,20 @@ bgmCompareOutput run_gibbs_sampler_bgmCompare(
           num_categories, missing_data_indices, is_ordinal_variable,
           baseline_category, rng
       );
+
+      // Imputation may change observations; refresh the sweep state.
+      initialize_sweep_state_observations(
+        sweep_state, observations, group_indices, num_groups
+      );
+      rebuild_sweep_state_weights(
+        sweep_state, pairwise_effects, pairwise_effect_indices,
+        inclusion_indicator, projection, num_groups
+      );
     }
 
     // Main Gibbs update step for parameters
     gibbs_update_step_bgmcompare (
-        observations, num_categories, counts_per_category,
+        observations, sweep_state, num_categories, counts_per_category,
         blume_capel_stats, inclusion_indicator,
         pairwise_effects, main_effects, is_ordinal_variable, baseline_category,
         iteration, pairwise_effect_indices, pairwise_stats, nuts_max_depth,
diff --git a/src/models/bgmCompare/bgmCompare_state.h b/src/models/bgmCompare/bgmCompare_state.h
new file mode 100644
index 00000000..41271897
--- /dev/null
+++ b/src/models/bgmCompare/bgmCompare_state.h
@@ -0,0 +1,85 @@
+#pragma once
+
+/**
+ * @file bgmCompare_state.h
+ * @brief Persistent per-group sweep state for the bgmCompare sampler.
+ */
+
+#include <RcppArmadillo.h>
+#include <vector>
+
+/**
+ * Persistent per-group state shared by the bgmCompare Gibbs sweeps.
+ *
+ * Holds the double-precision observation matrices, the group-specific
+ * effective pairwise weights, and the residual (rest-score) matrices
+ * `residual[g] = obs_double[g] * pairwise_group[g]`. The state is owned by
+ * run_gibbs_sampler_bgmCompare() and threaded by reference through the
+ * sweeps, which keep it consistent with the model parameters:
+ *  - accepted pairwise-effect moves and accepted difference-indicator flips
+ *    update `pairwise_group` and two residual columns per group;
+ *  - accepted main-effect moves leave the state untouched;
+ *  - imputation refills the observation matrices and triggers a full weight
+ *    rebuild once per pass;
+ *  - NUTS updates and graph re-initialisation trigger a full weight rebuild
+ *    via rebuild_sweep_state_weights().
+ *
+ * The state also caches the per-variable, per-group log-normalizer sums of
+ * the pseudolikelihood, `log_normalizer(v, g) = sum_i [bound_i + log
+ * denom_i]` at the current parameters and rest scores. A variable's row is
+ * valid while its main effects and its residual columns are unchanged:
+ * accepted Metropolis moves store the proposed-state normalizers, accepted
+ * indicator flips invalidate the affected variables, and weight rebuilds
+ * invalidate every variable.
+ */
+struct CompareSweepState {
+  arma::mat obs_double_all;              ///< Observations as double (n x V)
+  std::vector<arma::mat> obs_double;     ///< Per-group observations as double (n_g x V)
+  std::vector<arma::mat> pairwise_group; ///< Per-group effective pairwise weights (V x V)
+  std::vector<arma::mat> residual;       ///< Per-group rest scores obs_double[g] * pairwise_group[g] (n_g x V)
+  arma::mat log_normalizer;              ///< Cached pseudolikelihood log-normalizer sums (V x G)
+  arma::uvec normalizer_valid;           ///< 1 if a variable's log_normalizer row is current (V)
+};
+
+/**
+ * Fill the observation members of a CompareSweepState.
+ *
+ * Converts the integer observation matrix to double, both as a whole and
+ * per group. Does not touch the weight or residual members.
+ *
+ * @param[out] state          Sweep state to fill
+ * @param observations        Integer observation matrix (n x V)
+ * @param group_indices       Group start/end row indices per group (G x 2)
+ * @param num_groups          Number of groups (G)
+ */
+void initialize_sweep_state_observations(
+    CompareSweepState& state,
+    const arma::imat& observations,
+    const arma::imat& group_indices,
+    const int num_groups
+);
+
+/**
+ * Rebuild the effective pairwise weights and residual matrices.
+ *
+ * Recomputes `pairwise_group[g]` from the pairwise-effect matrix and the
+ * inclusion indicators, and `residual[g]` as one matrix product per group.
+ * Invalidates the normalizer cache for every variable. Called after
+ * wholesale parameter changes (NUTS updates, graph initialisation) and once
+ * at sampler start.
+ *
+ * @param[in,out] state           Sweep state with observation members filled
+ * @param pairwise_effects        Pairwise-effect matrix (rows = pairs, cols = groups)
+ * @param pairwise_effect_indices Row index per variable pair (V x V)
+ * @param inclusion_indicator     Difference inclusion indicators (V x V)
+ * @param projection              Group contrast matrix (G x (G-1))
+ * @param num_groups              Number of groups (G)
+ */
+void rebuild_sweep_state_weights(
+    CompareSweepState& state,
+    const arma::mat& pairwise_effects,
+    const arma::imat& pairwise_effect_indices,
+    const arma::imat& inclusion_indicator,
+    const arma::mat& projection,
+    const int num_groups
+);