From 3413a742285664a97524f3e5fcd9d30242dd380a Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 10 Jun 2026 03:36:07 +0000 Subject: [PATCH] perf(ruview-gamma): factorize GP once per recommendation (-81%) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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 --- CHANGELOG.md | 1 + v2/crates/ruview-gamma/README.md | 6 +- v2/crates/ruview-gamma/src/optimizer.rs | 133 +++++++++++++++++------- 3 files changed, 99 insertions(+), 41 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 347c0a84..d678006f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -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. diff --git a/v2/crates/ruview-gamma/README.md b/v2/crates/ruview-gamma/README.md index 0d7679aa..f6dbc750 100644 --- a/v2/crates/ruview-gamma/README.md +++ b/v2/crates/ruview-gamma/README.md @@ -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) diff --git a/v2/crates/ruview-gamma/src/optimizer.rs b/v2/crates/ruview-gamma/src/optimizer.rs index 3be50618..43e80d89 100644 --- a/v2/crates/ruview-gamma/src/optimizer.rs +++ b/v2/crates/ruview-gamma/src/optimizer.rs @@ -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> { 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 = (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, + /// `alpha = K⁻¹ y`. + alpha: Vec, + 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 = (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)]