Performance follow-up: GGM chol deferral, NUTS gradient workspaces, mixed AM caching, bgmCompare trims

src/ggm_gradient_interface.cpp

@@ -58,7 +58,7 @@ Rcpp::List ggm_test_forward_map(
     GGMGradientEngine engine;
     engine.rebuild(cs, 100, dummy_S, ip, dp);
 
-    ForwardMapResult fm = engine.forward_map(theta);
+    const ForwardMapResult& fm = engine.forward_map(theta);
 
     return Rcpp::List::create(
         Rcpp::Named("Phi") = Rcpp::wrap(fm.Phi),

src/models/bgmCompare/bgmCompare_logp_and_grad.cpp

@@ -835,16 +835,15 @@ double log_pseudoposterior_main_component(
   }
 
   const int num_variables = observations.n_cols;
-  const int max_num_categories = num_categories.max();
   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];
-    const arma::imat blume_capel_stats = blume_capel_stats_group[group];
-
-    arma::mat main_group(num_variables, max_num_categories, arma::fill::zeros);
-    arma::mat pairwise_group(num_variables, num_variables, arma::fill::zeros);
+    const arma::imat& counts_per_category = counts_per_category_group[group];
+    const arma::imat& blume_capel_stats = blume_capel_stats_group[group];
 
     const arma::vec proj_g = projection.row(group).t(); // length = num_groups-1
 
@@ -853,27 +852,21 @@ double log_pseudoposterior_main_component(
       variable, num_groups, main_effects, main_effect_indices, proj_g
     );
 
-    // store into row v
-    main_group(variable, arma::span(0, me.n_elem - 1)) = me.t();
-
-    // upper triangle incl. base value; mirror to keep symmetry
+    // pairwise weights with the other variables; entry `variable` stays zero
     for (int u = 0; u < num_variables; u++) {
       if(u == variable) continue;
-      double w = compute_group_pairwise_effects(
+      pairwise_col(u) = compute_group_pairwise_effects(
         variable, u, num_groups, pairwise_effects, pairwise_effect_indices,
         inclusion_indicator, proj_g
       );
-      pairwise_group(variable, u) = w;
-      pairwise_group(u, variable) = w;
     }
 
     // ---- data contribution pseudolikelihood (linear terms) ----
     if (is_ordinal_variable(variable)) {
-      const double val = main_group(variable, category);
-      log_pp += static_cast<double>(counts_per_category(category, variable)) * val;
+      log_pp += static_cast<double>(counts_per_category(category, variable)) *
+        me(category);
     } else {
-      log_pp += static_cast<double>(blume_capel_stats(par, variable)) *
-        main_group(variable, par);
+      log_pp += static_cast<double>(blume_capel_stats(par, variable)) * me(par);
     }
 
     // ---- data contribution pseudolikelihood (quadratic terms) ----
@@ -882,7 +875,7 @@ double log_pseudoposterior_main_component(
     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_group.col(variable);
+    const arma::vec rest_score = obs * pairwise_col;
     const int num_cats = num_categories(variable);
 
     // bound to stabilize exp; clamp at 0 (the reference-category exponent) so
@@ -892,14 +885,12 @@ double log_pseudoposterior_main_component(
     arma::vec denom(rest_score.n_elem, arma::fill::zeros);
 
     if (is_ordinal_variable(variable)) {
-      arma::vec main_eff = main_group.row(variable).cols(0, num_cats - 1).t();
       denom = compute_denom_ordinal(
-        rest_score, main_eff, bound
+        rest_score, me, bound
       );
     } else {
-      // linear/quadratic main effects from main_group
-      const double lin_effect  = main_group(variable, 0);
-      const double quad_effect = main_group(variable, 1);
+      const double lin_effect  = me(0);
+      const double quad_effect = me(1);
       const int ref = baseline_category(variable);
 
       denom = compute_denom_blume_capel(
@@ -959,8 +950,7 @@ double log_pseudoposterior_main_component(
 //  - main_effect_indices: Index ranges [row_start, row_end] for each variable.
 //  - pairwise_effect_indices: Lookup table mapping (var1, var2) to 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.
 //  - num_categories: Number of categories per variable.
 //  - pairwise_stats_group: Per-group pairwise sufficient statistics.
 //  - residual_matrices: Per-group residual matrices (persons × variables).
@@ -987,8 +977,7 @@ double log_pseudoposterior_pair_component(
     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 arma::ivec& num_categories,
     const std::vector<arma::mat>& pairwise_stats_group,
     const std::vector<arma::mat>& residual_matrices,
@@ -1007,7 +996,7 @@ double log_pseudoposterior_pair_component(
     return 0.0;
   }
 
-  const int num_variables = observations.n_cols;
+  const int num_variables = obs_double_groups[0].n_cols;
   const int max_num_categories = num_categories.max();
   double log_pp = 0.0;
   int idx = pairwise_effect_indices(variable1, variable2);
@@ -1042,19 +1031,14 @@ double log_pseudoposterior_pair_component(
     }
 
     // ---- pseudolikelihood normalizing constants (using residual matrix + delta) ----
-    const int r0 = group_indices(group, 0);
-    const int r1 = group_indices(group, 1);
-
-    // Pre-convert observation columns for the two variables (needed for delta adjustment)
-    const arma::vec obs_var1 = arma::conv_to<arma::vec>::from(observations.col(variable1).rows(r0, r1));
-    const arma::vec obs_var2 = arma::conv_to<arma::vec>::from(observations.col(variable2).rows(r0, r1));
+    const arma::mat& obs_g = obs_double_groups[group];
 
     for (int v : {variable1, variable2}) {
       const int num_cats = num_categories(v);
-      const arma::vec& obs_other = (v == variable1) ? obs_var2 : obs_var1;
+      const int other = (v == variable1) ? variable2 : variable1;
 
       // Use residual_matrix with delta adjustment: O(n) instead of O(n*p)
-      arma::vec rest_score = residual_matrices[group].col(v) + obs_other * delta_g;
+      arma::vec rest_score = residual_matrices[group].col(v) + obs_g.col(other) * delta_g;
 
       // bound to stabilize exp; clamp at 0 so exp cannot overflow (the
       // Blume-Capel branch overwrites bound with its own max).

src/models/bgmCompare/bgmCompare_logp_and_grad.h

@@ -211,6 +211,7 @@ double log_pseudoposterior_main_component(
  * Uses pre-computed residual matrices adjusted by delta to avoid full
  * recomputation. Used by element-wise Metropolis updates.
  *
+ * @param obs_double_groups  Per-group observation matrices converted to double
  * @param residual_matrices  Pre-computed residual matrices per group
  * @param variable1          First variable index
  * @param variable2          Second variable index
@@ -224,8 +225,7 @@ double log_pseudoposterior_pair_component(
     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 arma::ivec& num_categories,
     const std::vector<arma::mat>& pairwise_stats_group,
     const std::vector<arma::mat>& residual_matrices,

src/models/bgmCompare/bgmCompare_sampler.cpp

@@ -58,9 +58,10 @@
 //
 // Notes:
 //  - The function updates both raw data and sufficient statistics in-place.
-//  - Group-specific pairwise effects are recomputed per missing entry.
-//  - For efficiency, you may consider incremental updates to `pairwise_stats`
-//    instead of full recomputation (`obs.t() * obs`) after each change.
+//  - Group-specific pairwise effects are built once per group; parameters
+//    are constant across the imputation pass.
+//  - `pairwise_stats` is updated incrementally: a changed cell only alters
+//    row and column `variable` of the group crossproduct.
 void impute_missing_bgmcompare(
     const arma::mat& main_effects,
     const arma::mat& pairwise_effects,
@@ -89,6 +90,24 @@ void impute_missing_bgmcompare(
   double exponent, cumsum, u;
   int score, person, variable, new_value, old_value, group;
 
+  // Group-specific pairwise effect matrices; the parameters are constant
+  // across the imputation pass, so one build per group suffices.
+  std::vector<arma::mat> group_pairwise_effects(num_groups);
+  for(int g = 0; g < num_groups; g++) {
+    const arma::vec proj_g = projection.row(g).t();
+    group_pairwise_effects[g].zeros(num_variables, num_variables);
+    for(int v1 = 0; v1 < num_variables-1; v1++) {
+      for(int v2 = v1 + 1; v2 < num_variables; v2++) {
+        double w = compute_group_pairwise_effects(
+            v1, v2, num_groups, pairwise_effects, pairwise_effect_indices,
+            inclusion_indicator, proj_g
+        );
+        group_pairwise_effects[g](v1, v2) = w;
+        group_pairwise_effects[g](v2, v1) = w;
+      }
+    }
+  }
+
   //Impute missing data
   for(int missing = 0; missing < num_missings; missing++) {
     // Identify the observation to impute
@@ -101,21 +120,8 @@ void impute_missing_bgmcompare(
     arma::vec group_main_effects = compute_group_main_effects(
       variable, num_groups, main_effects,  main_effect_indices, proj_g);
 
-    // Generate a new observation based on the model
-    arma::mat group_pairwise_effects(num_variables, num_variables, arma::fill::zeros);
-    for(int v1 = 0; v1 < num_variables-1; v1++) {
-      for(int v2 = v1 + 1; v2 < num_variables; v2++) {
-        double w = compute_group_pairwise_effects(
-            v1, v2, num_groups, pairwise_effects, pairwise_effect_indices,
-            inclusion_indicator, proj_g
-        );
-        group_pairwise_effects(v1, v2) = w;
-        group_pairwise_effects(v2, v1) = w;
-      }
-    }
-
     double rest_score =
-      arma::as_scalar(observations.row(person) * group_pairwise_effects.col(variable));
+      arma::as_scalar(observations.row(person) * group_pairwise_effects[group].col(variable));
 
     if(is_ordinal_variable[variable] == true) {
       // For regular binary or ordinal variables
@@ -159,27 +165,32 @@ void impute_missing_bgmcompare(
 
       // Update sufficient statistics for main effects
       if(is_ordinal_variable[variable] == true) {
-        arma::imat counts_per_category_group = counts_per_category[group];
+        arma::imat& counts_per_category_group = counts_per_category[group];
         if(old_value > 0)
           counts_per_category_group(old_value-1, variable)--;
         if(new_value > 0)
           counts_per_category_group(new_value-1, variable)++;
-        counts_per_category[group] = counts_per_category_group;
       } else {
-        arma::imat blume_capel_stats_group = blume_capel_stats[group];
+        arma::imat& blume_capel_stats_group = blume_capel_stats[group];
         blume_capel_stats_group(0, variable) -= old_value;
         blume_capel_stats_group(0, variable) += new_value;
         blume_capel_stats_group(1, variable) -= old_value * old_value;
         blume_capel_stats_group(1, variable) += new_value * new_value;
-        blume_capel_stats[group] = blume_capel_stats_group;
       }
 
-      // Update sufficient statistics for pairwise effects
-      const int r0 = group_indices(group, 0);
-      const int r1 = group_indices(group, 1);
-      arma::mat obs = arma::conv_to<arma::mat>::from(observations.rows(r0, r1));
-      arma::mat pairwise_stats_group = obs.t() * obs; // crossprod
-      pairwise_stats[group] = pairwise_stats_group;
+      // Update sufficient statistics for pairwise effects. In the group
+      // crossproduct X.t() * X only row and column `variable` depend on the
+      // changed cell: with delta = new_value - old_value and x the person's
+      // row (which already holds new_value), entry (u, variable) gains
+      // delta * x(u), symmetrically; the diagonal gain 2 * delta * new_value
+      // - delta^2 equals new_value^2 - old_value^2.
+      const double delta = static_cast<double>(new_value - old_value);
+      const arma::rowvec obs_row =
+        arma::conv_to<arma::rowvec>::from(observations.row(person));
+      arma::mat& pairwise_stats_group = pairwise_stats[group];
+      pairwise_stats_group.col(variable) += delta * obs_row.t();
+      pairwise_stats_group.row(variable) += delta * obs_row;
+      pairwise_stats_group(variable, variable) -= delta * delta;
     }
   }
   return;
@@ -435,7 +446,7 @@ void update_pairwise_effects_metropolis_bgmcompare (
       double delta = theta - current;
       return log_pseudoposterior_pair_component(
         main_effects, pairwise_effects, main_effect_indices,
-        pairwise_effect_indices, projection, observations, group_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,
@@ -991,7 +1002,7 @@ void tune_proposal_sd_bgmcompare(
           return log_pseudoposterior_pair_component(
             main_effects, pairwise_effects,
             main_effect_indices, pairwise_effect_indices,
-            projection, observations, group_indices,
+            projection, obs_double_groups,
             num_categories, pairwise_stats, residual_matrices, num_groups,
             inclusion_indicator, is_ordinal_variable, baseline_category,
             v1, v2, h, delta,