perf(ruview-gamma): factorize GP once per recommendation (-81%)

BayesianOptimizer::recommend evaluated Expected Improvement at every 0.1 Hz grid candidate, and each predict() rebuilt the kernel matrix and re-ran Cholesky — ~82 factorizations per call — though K and alpha=K^-1 y depend only on the observations, not the query point. Fit the GP once (GpFit: cached L + alpha, lower-triangle-only K build) and reuse it across the grid. Bit-identical arithmetic: the pinned deterministic witness and all 64 tests are unchanged; pure work elimination. Measured (criterion, paired): gamma_bayesian_recommend 105us -> 19us (-81%); gamma_calibration_sweep 466us -> 135us (-71%). https://claude.ai/code/session_01MjBucx95K4BuUxZi8NWwRH
2026-07-20 17:03:24 +00:00 · 2026-06-10 03:36:07 +00:00
parent 3408883f79
commit 3413a74228
3 changed files with 99 additions and 41 deletions
@@ -12,6 +12,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - **RuView beyond-SOTA research series** (`docs/research/ruview-beyond-sota/`, 6 docs) — research-swarm output defining the beyond-SOTA bar and the path to it: system capability audit (role→crate maturity matrix, gap analysis, risk register), web-verified 2026 SOTA landscape per capability axis (incl. ratified IEEE 802.11bf-2025), 8-pillar target architecture on the ADR-136 contract spine (no rewrite), 6-layer benchmark/validation methodology (all 15 criterion bench targets inventoried; ADR-149 statistical protocol), and a determinism-safe optimization roadmap. Includes session validation evidence: 2,797 workspace tests / 0 failed, Python proof PASS (bit-exact), paired pre/post criterion runs.

 ### Performance
+- **`ruview-gamma` Bayesian optimizer factorized once per recommendation.** `BayesianOptimizer::recommend` evaluated Expected Improvement at every 0.1 Hz grid candidate, and each `predict` rebuilt the full kernel matrix and re-ran Cholesky — ~82 factorizations per recommendation, though `K` and `alpha = K⁻¹y` depend only on the observations, not the query point. Now the GP is fit once (`GpFit`: cached `L` + `alpha`, lower-triangle-only `K` build) and reused across the grid. Bit-identical arithmetic (the pinned deterministic witness `13cb164c…` and all 64 tests are unchanged), pure work elimination. Measured (criterion, paired): `gamma_bayesian_recommend` 105 µs → 19 µs (**−81%**), `gamma_calibration_sweep` (9 recommends) 466 µs → 135 µs (**−71%**).
 - **CIR estimator warm-start precompute** — the diagonal Tikhonov preconditioner `diag(Φ^H Φ)+λI` and its CSR matrix were rebuilt every frame although they depend only on Φ and λ (fixed at `CirEstimator::new`); now precomputed at construction (`ruvsense/cir.rs`). Bit-identical floats (summation order unchanged, witness chain unaffected). Measured: `cir_estimate/he40` −3.9% (p<0.01), multiband groups −1.2/−1.4%; smaller configs within container noise.
 - **RF tomography solver hoisting** — ISTA gradient buffer no longer allocated inside the 100-iteration loop, and the Frobenius Lipschitz bound moved from per-`reconstruct` to construction (`ruvsense/tomography.rs`). Bit-identical results.

@@ -88,9 +88,9 @@ the optimizer, simulator, response update, or session hashing.
 | Bench | Median | Note |
 |-------|--------|------|
 | `gamma_safety_tick` | ~9.3 ns | vs ADR-250 §17 < 500 ms hard-stop latency bound |
-| `gamma_bandit_select` | ~73 ns | LinUCB decision |
-| `gamma_bayesian_recommend` | ~105 µs | GP + EI over the 0.1 Hz envelope grid |
-| `gamma_calibration_sweep` | ~486 µs | full 9-session enroll → simulate → score → update → witness |
+| `gamma_bandit_select` | ~74 ns | LinUCB decision |
+| `gamma_bayesian_recommend` | ~19 µs | GP + EI over the 0.1 Hz envelope grid (was ~105 µs: the GP is now factorized once per recommend, not once per grid candidate — −81%, bit-identical) |
+| `gamma_calibration_sweep` | ~135 µs | full 9-session enroll → simulate → score → update → witness (was ~486 µs, −71%) |

 ## Roadmap (ADR-250 §21)

@@ -102,17 +102,23 @@ impl BayesianOptimizer {
        best
    }

-    /// GP posterior `(mean, variance)` at `x`. Falls back to `(0, signal_var)`
-    /// (the prior) when there are no observations or K is not SPD.
-    pub fn predict(&self, x: f64) -> (f64, f64) {
+    /// Factorize the GP once: Cholesky `L` of `K = RBF(X,X)+noise·I` and the
+    /// weight vector `alpha = K⁻¹ y`. Both depend only on the observations, not
+    /// on any query point, so a single fit serves the whole acquisition grid.
+    /// Returns `None` when there are no observations or `K` is not SPD.
+    ///
+    /// The per-query arithmetic in [`GpFit::predict`] is identical to the old
+    /// inline path, so predictions (and therefore the session witness) are
+    /// bit-for-bit unchanged — this is a pure work-elimination optimization.
+    fn fit(&self) -> Option<GpFit<'_>> {
        let n = self.obs_x.len();
        if n == 0 {
-            return (0.0, self.signal_var);
+            return None;
        }
-        // K = RBF(X,X) + noise·I
+        // K (lower triangle is all Cholesky reads) = RBF(X,X) + noise·I.
        let mut k = vec![0.0f64; n * n];
        for i in 0..n {
-            for j in 0..n {
+            for j in 0..=i {
                let mut v = rbf_kernel(
                    &[self.obs_x[i]],
                    &[self.obs_x[j]],
@@ -125,22 +131,27 @@ impl BayesianOptimizer {
                k[i * n + j] = v;
            }
        }
-        let l = match cholesky(&k, n) {
-            Some(l) => l,
-            None => return (0.0, self.signal_var),
-        };
-        // k* = RBF(X, x)
-        let kstar: Vec<f64> = (0..n)
-            .map(|i| rbf_kernel(&[self.obs_x[i]], &[x], self.length_scale, self.signal_var))
-            .collect();
-        // alpha = K^-1 y  (solve L Lᵀ alpha = y)
+        let l = cholesky(&k, n)?;
+        // alpha = K⁻¹ y  (solve L Lᵀ alpha = y) — computed once.
        let y1 = forward_subst(&l, &self.obs_y, n);
        let alpha = back_subst_transpose(&l, &y1, n);
-        let mean = dot(&kstar, &alpha);
-        // var = k(x,x) - v·v, where L v = k*
-        let v = forward_subst(&l, &kstar, n);
-        let var = (self.signal_var - dot(&v, &v)).max(0.0);
-        (mean, var)
+        Some(GpFit {
+            obs_x: &self.obs_x,
+            l,
+            alpha,
+            n,
+            length_scale: self.length_scale,
+            signal_var: self.signal_var,
+        })
+    }
+
+    /// GP posterior `(mean, variance)` at `x`. Falls back to `(0, signal_var)`
+    /// (the prior) when there are no observations or K is not SPD.
+    pub fn predict(&self, x: f64) -> (f64, f64) {
+        match self.fit() {
+            Some(fit) => fit.predict(x),
+            None => (0.0, self.signal_var),
+        }
    }

    /// Expected Improvement (for maximization) at `x`.
@@ -149,44 +160,48 @@ impl BayesianOptimizer {
            Some((_, by)) => by,
            None => return self.signal_var.sqrt(), // pure exploration
        };
-        let (mu, var) = self.predict(x);
-        let sigma = var.sqrt();
-        if sigma <= 1e-12 {
-            return 0.0;
+        match self.fit() {
+            Some(fit) => fit.expected_improvement(x, best, self.xi),
+            None => self.signal_var.sqrt(),
        }
-        let imp = mu - best - self.xi;
-        let z = imp / sigma;
-        imp * normal_cdf(z) + sigma * normal_pdf(z)
    }

    /// Recommend the next stimulus by maximizing EI over the envelope's 0.1 Hz
    /// grid, holding `base`'s intensity (clamped). The result is guaranteed
    /// inside the envelope. With no observations it returns the 40 Hz prior.
+    ///
+    /// Fits the GP **once** and reuses the factorization across every grid
+    /// candidate (was: a full Cholesky per candidate, ~82× the work).
    pub fn recommend(
        &self,
        envelope: &SafetyEnvelope,
        base: &StimulusParameters,
    ) -> Recommendation {
-        if self.obs_x.is_empty() {
-            let s = envelope.clamp(*base);
-            return Recommendation {
-                stimulus: s,
-                expected_improvement: 0.0,
-                predicted_score: 0.0,
-                confidence: 0.0,
-            };
-        }
+        let fit = match self.fit() {
+            Some(f) => f,
+            None => {
+                let s = envelope.clamp(*base);
+                return Recommendation {
+                    stimulus: s,
+                    expected_improvement: 0.0,
+                    predicted_score: 0.0,
+                    confidence: 0.0,
+                };
+            }
+        };
+        // best() is Some here (fit exists ⇒ ≥1 observation).
+        let best = self.best().map(|(_, by)| by).unwrap_or(0.0);
        let grid = fine_grid(envelope);
        let mut best_f = base.frequency_hz;
        let mut best_ei = f64::NEG_INFINITY;
        for &f in &grid {
-            let ei = self.expected_improvement(f);
+            let ei = fit.expected_improvement(f, best, self.xi);
            if ei > best_ei {
                best_ei = ei;
                best_f = f;
            }
        }
-        let (mu, var) = self.predict(best_f);
+        let (mu, var) = fit.predict(best_f);
        let mut s = *base;
        s.frequency_hz = best_f;
        let s = envelope.clamp(s);
@@ -201,6 +216,48 @@ impl BayesianOptimizer {
    }
 }

+/// A cached GP factorization (Cholesky `L` + weights `alpha`) over a fixed set
+/// of observations, reused across many query points in one acquisition pass.
+struct GpFit<'a> {
+    obs_x: &'a [f64],
+    /// Cholesky factor of `K` (lower-triangular, `n×n` row-major).
+    l: Vec<f64>,
+    /// `alpha = K⁻¹ y`.
+    alpha: Vec<f64>,
+    n: usize,
+    length_scale: f64,
+    signal_var: f64,
+}
+
+impl GpFit<'_> {
+    /// Posterior `(mean, variance)` at `x`. Bit-identical to the former inline
+    /// computation; only the shared `L`/`alpha` are now precomputed.
+    fn predict(&self, x: f64) -> (f64, f64) {
+        let n = self.n;
+        // k* = RBF(X, x)
+        let kstar: Vec<f64> = (0..n)
+            .map(|i| rbf_kernel(&[self.obs_x[i]], &[x], self.length_scale, self.signal_var))
+            .collect();
+        let mean = dot(&kstar, &self.alpha);
+        // var = k(x,x) - v·v, where L v = k*
+        let v = forward_subst(&self.l, &kstar, n);
+        let var = (self.signal_var - dot(&v, &v)).max(0.0);
+        (mean, var)
+    }
+
+    /// Expected Improvement at `x` given the incumbent `best` and margin `xi`.
+    fn expected_improvement(&self, x: f64, best: f64, xi: f64) -> f64 {
+        let (mu, var) = self.predict(x);
+        let sigma = var.sqrt();
+        if sigma <= 1e-12 {
+            return 0.0;
+        }
+        let imp = mu - best - xi;
+        let z = imp / sigma;
+        imp * normal_cdf(z) + sigma * normal_pdf(z)
+    }
+}
+
 /// A recommendation plus its explainability fields (ADR-250 §14 response,
 /// §16 "explainable recommendation").
 #[derive(Debug, Clone, Copy, PartialEq)]