rqe-optimizer: fall back to each sketch's guarantee where the study measured nothing

docs/rqe_sketch_deployment_v1.md

@@ -240,6 +240,26 @@ lookback: series per group times scrapes per lookback.
   last N, the merged answer has no accuracy either (#158).
   `SaturationCurves::load` refuses a study with no merge curves for a
   candidate KLL or top-k sketch.
+- Where the study measured nothing for a bracketing point at or above its
+  first N (past an unsaturated curve, past the measured merge counts or N,
+  or with no merge curve), the accuracy falls back to the algorithm's
+  guarantee at 95% confidence (`rqe_optimizer::theory`; two-sided for
+  rank and relative errors, one-sided for CMS overestimates). It is never
+  better than the worst last measurement of the curves it extends. The fallbacks are:
+  - KLL: rank error `2.296 / k^0.9723` (DataSketches' 99% value), scaled
+    to 95% two-sided;
+  - DDSketch: relative value error `α` (unbounded store);
+  - HLL: `1.96 · 1.04 / sqrt(2^lg_k)`;
+  - CMS-heap top-k, unmerged: the share of the top `k` Zipf(θ, K) keys kept
+    ranked for every key at once. With Markov per row, independent rows,
+    and a union bound over the `K − k` other keys, key `i` fails with
+    probability at most `(cols · (p_i − p_{k+1}))^−rows`. Keys are counted
+    from the most frequent while the summed failure stays within 5%. It
+    assumes the heap ranks keys by their final CMS estimates.
+
+  A sketch with no published bound, a merged top-k answer, a point below
+  the first N, or a shape outside the grid stays unknown.
+  `accuracy_with_source` reports which points used a bound.
 
 The cost table's sketch accuracies are not read, but they are one point on
 each curve: a row's `measured_at` (items, Zipf θ and population, or Pareto

rqe-optimizer/src/lib.rs

@@ -16,6 +16,7 @@ pub mod enumerate;
 pub mod milp;
 pub mod pareto;
 pub mod saturation;
+pub mod theory;
 
 use std::collections::BTreeMap;
 use std::collections::BTreeSet;

rqe-optimizer/src/saturation.rs

@@ -19,6 +19,7 @@ use std::io;
 use std::path::Path;
 
 use crate::autosketch::window_adapter;
+use crate::theory;
 use crate::{
     accuracy_key, family_properties, table_accuracy, AccuracyDirection, AtomicCostEntry,
     Capability, Deployment, LabelSet, MetricFacts, Raqe, WorkloadFacts,
@@ -120,6 +121,47 @@ impl GridPoint {
         curve_error_at(&self.curve, self.n_sat, n, direction)
     }
 
+    /// The curves an answer merged from `merges` instances reads: the plain
+    /// curve alone, or the measured shard counts either side of `merges`
+    /// (one is the plain curve). Past the largest count, only the largest.
+    fn bracketing_curves(&self, merges: u64) -> Vec<&[(f64, f64, f64)]> {
+        if merges <= 1 {
+            return vec![&self.curve];
+        }
+        let below = self
+            .merged
+            .range(..=merges)
+            .next_back()
+            .map_or(&self.curve[..], |(_, c)| &c[..]);
+        match self.merged.range(merges..).next() {
+            Some((_, above)) => vec![below, above],
+            None => self
+                .merged
+                .values()
+                .next_back()
+                .map_or(vec![below], |c| vec![&c[..]]),
+        }
+    }
+
+    /// The smallest N the study measured the curves an answer merged from
+    /// `merges` instances reads.
+    fn first_measured_n(&self, merges: u64) -> Option<f64> {
+        self.bracketing_curves(merges)
+            .iter()
+            .map(|c| c.first().map(|p| p.0))
+            .try_fold(0.0, |max: f64, first| first.map(|f| max.max(f)))
+    }
+
+    /// The worst last error measured on the curves an answer merged from
+    /// `merges` instances reads: the floor a fallback extending them can't
+    /// beat.
+    fn last_measured_error(&self, merges: u64, direction: AccuracyDirection) -> Option<f64> {
+        self.bracketing_curves(merges)
+            .iter()
+            .filter_map(|c| c.last().map(|p| p.1))
+            .reduce(|a, b| worse(a, b, direction))
+    }
+
     /// The error at `n` items of the sketch merged from `merges` instances:
     /// the worse of the measured shard counts either side of `merges` (one
     /// is the plain sketch). `None` past the largest measured count, past the
@@ -196,6 +238,15 @@ fn worse(a: f64, b: f64, direction: AccuracyDirection) -> f64 {
     }
 }
 
+/// Where an accuracy came from.
+#[derive(Debug, Clone, Copy, PartialEq, Eq)]
+pub enum AccuracySource {
+    /// The cost table (exact accumulators) or the study's curves.
+    Measured,
+    /// The algorithm's published guarantee, where the study measured nothing.
+    Theory,
+}
+
 /// What [`SaturationCurves::check_cost_table`] found, one line per row.
 #[derive(Debug, Default, PartialEq)]
 pub struct CostTableCheck {
@@ -342,8 +393,24 @@ impl SaturationCurves {
         deployment: &Deployment,
         facts: &WorkloadFacts,
     ) -> Option<f64> {
+        self.accuracy_with_source(raqe, deployment, facts)
+            .map(|(value, _)| value)
+    }
+
+    /// [`SaturationCurves::accuracy`], and where the returned (worst) value
+    /// came from: measured, or the algorithm's guarantee
+    /// ([`crate::theory`]) for a bracketing point the study didn't measure
+    /// (past an unsaturated curve, past the measured merge counts or N, or
+    /// with no merge curve). Below the curves' first N, or outside the grid,
+    /// it stays unknown.
+    pub fn accuracy_with_source(
+        &self,
+        raqe: &Raqe,
+        deployment: &Deployment,
+        facts: &WorkloadFacts,
+    ) -> Option<(f64, AccuracySource)> {
         if family_properties(&deployment.config.sketch).exact {
-            return table_accuracy(raqe, deployment);
+            return table_accuracy(raqe, deployment).map(|v| (v, AccuracySource::Measured));
         }
         let metric_facts = &facts[&deployment.metric];
         let grouping = &deployment.grouping_labels;
@@ -356,10 +423,26 @@ impl SaturationCurves {
         } else {
             1
         };
-        let mut worst = None;
+        let params = config_params(&deployment.config)?;
+        let mut worst: Option<(f64, AccuracySource)> = None;
         for point in points {
-            let error = point.merged_error_at(merges, covered, direction)?;
-            worst = Some(worst.map_or(error, |w| worse(w, error, direction)));
+            let found = match point.merged_error_at(merges, covered, direction) {
+                Some(error) => (error, AccuracySource::Measured),
+                None if covered >= point.first_measured_n(merges)? => {
+                    let bound =
+                        theory::bound(&deployment.config.sketch, &params, point.shape, merges)?;
+                    // No better than the last measurement it extends.
+                    let error = match point.last_measured_error(merges, direction) {
+                        Some(last) => worse(bound, last, direction),
+                        None => bound,
+                    };
+                    (error, AccuracySource::Theory)
+                }
+                None => return None,
+            };
+            if worst.is_none_or(|(w, _)| worse(w, found.0, direction) != w) {
+                worst = Some(found);
+            }
         }
         worst
     }
@@ -805,6 +888,50 @@ mod tests {
         assert!(missing.unwrap_err().to_string().contains("no merge curve"));
     }
 
+    /// Past an unsaturated curve the study measured nothing: the answer
+    /// falls back to the algorithm's guarantee, and says so. A merged top-k
+    /// answer has no guarantee.
+    #[test]
+    fn an_unmeasured_point_falls_back_to_the_guarantee() {
+        // θ = 1.0, K = 1e5 never saturates; 10 series × 1e8 s = 1e9 items.
+        let facts = workload(10, shape(1.0, 1e5));
+        let r = topk_raqe(100_000_000);
+        let one_window = curves().accuracy_with_source(&r, &deployment(TOPK, 100_000_000), &facts);
+        let params = BTreeMap::from([("rows".to_string(), 3.0), ("cols".to_string(), 1024.0)]);
+        let zipf = MeasuredShape::Zipf {
+            skew: 1.0,
+            keys: 1e5,
+        };
+        let bound = theory::bound(TOPK, &params, zipf, 1).unwrap();
+        assert_eq!(one_window, Some((bound, AccuracySource::Theory)));
+        let merged = curves().accuracy_with_source(&r, &deployment(TOPK, 10_000_000), &facts);
+        assert_eq!(merged, None);
+        // Measured where the curve has it.
+        let measured = curves().accuracy_with_source(
+            &r,
+            &deployment(TOPK, 100_000_000),
+            &workload(10, shape(1.2, 1e3)),
+        );
+        assert_eq!(
+            measured.map(|(_, source)| source),
+            Some(AccuracySource::Measured)
+        );
+    }
+
+    /// A fallback's floor is the worse last measurement of both shard
+    /// counts around the merge count, not only the smaller.
+    #[test]
+    fn a_fallback_floor_reads_both_bracketing_merge_curves() {
+        let point = &curves().points_by_sketch[TOPK][1];
+        let higher = AccuracyDirection::HigherIsBetter;
+        // 4 shards end at 0.88, 16 at 0.80: 8 merges floor at 0.80.
+        assert!((point.last_measured_error(8, higher).unwrap() - 0.8).abs() < 1e-12);
+        assert!((point.last_measured_error(4, higher).unwrap() - 0.88).abs() < 1e-12);
+        assert!((point.last_measured_error(1, higher).unwrap() - 0.9).abs() < 1e-12);
+        // Past the largest count, the largest's.
+        assert!((point.last_measured_error(100, higher).unwrap() - 0.8).abs() < 1e-12);
+    }
+
     #[test]
     fn exact_accumulators_keep_the_cost_table_value() {
         let r = topk_raqe(1_000_000);

rqe-optimizer/src/theory.rs

@@ -0,0 +1,180 @@
+//! Each sketch algorithm's published accuracy guarantee, at 95% confidence,
+//! in the metric the optimizer checks it by. [`crate::saturation`] falls back
+//! to it where the study measured nothing: past an unsaturated curve, past
+//! the measured merge counts or N, or with no merge curve for a point.
+//!
+//! The bounds hold for any stream length, so they don't depend on N. They are
+//! looser than the measured curves.
+
+use std::collections::BTreeMap;
+
+use aqpbm_core::MeasuredShape;
+
+use crate::analytical_cost_model::TOPK_ENTRIES;
+
+/// Confidence the bounds are stated at. Rank and relative errors are
+/// absolute, so their bounds are two-sided; CMS only overestimates, so the
+/// top-k bound is one-sided.
+pub const CONFIDENCE: f64 = 0.95;
+/// The standard normal quantile at [`CONFIDENCE`], two-sided.
+const Z: f64 = 1.96;
+
+/// The bound on `sketch`'s error with `params` on data shaped `shape`,
+/// answering from `merges` merged instances, or `None` when the algorithm
+/// has none at [`CONFIDENCE`]:
+///
+/// - `kll-percall`: normalized rank error of one quantile. DataSketches
+///   publishes `2.296 / k^0.9723` at 99%; KLL's error scales with
+///   `sqrt(ln(2/δ))` (two-sided), which gives the 95% value. KLL merges without loss of
+///   this guarantee.
+/// - `dd`: relative value error `α` for an unbounded store (sketch-bench's
+///   DDSketch collapses no buckets); merges exactly. The measured metric can
+///   exceed `α` (zero-bucket values, rank conventions), which the caller's
+///   floor at the last measurement covers.
+/// - `hll`: relative error `z · 1.04 / sqrt(2^lg_k)`, `z` the two-sided
+///   normal quantile; merges exactly.
+/// - `cms-heap-topk-fastpath-vector2d`: precision@k on Zipf(θ) over `K`
+///   keys, for every key at once. A CMS row overshoots a key by more than
+///   `t` with probability at most `N / (w·t)` (Markov), and the `d` rows are
+///   independent, so a non-top key overtakes top key `i` with probability at
+///   most `(w·(p_i − p_{k+1}))^−d`. A union bound over the `K − k` non-top
+///   keys, then over the counted top keys, keeps the total failure within
+///   `1 − CONFIDENCE`: precision is the share of the top `k` counted,
+///   taken from the most frequent down. It assumes the heap ranks keys by
+///   their final CMS estimates. A merged answer has none: each shard's heap
+///   keeps only its own top `k`.
+pub fn bound(
+    sketch: &str,
+    params: &BTreeMap<String, f64>,
+    shape: MeasuredShape,
+    merges: u64,
+) -> Option<f64> {
+    match sketch {
+        "kll-percall" => {
+            let k = params.get("k")?;
+            let at_99 = 2.296 / k.powf(0.9723);
+            // Two-sided tails: ln(2/δ).
+            Some(at_99 * ((2.0 / (1.0 - CONFIDENCE)).ln() / 200f64.ln()).sqrt())
+        }
+        "dd" => params.get("alpha").copied(),
+        "hll" => {
+            let registers = 2f64.powf(*params.get("lg_k")?);
+            Some(Z * 1.04 / registers.sqrt())
+        }
+        "cms-heap-topk-fastpath-vector2d" if merges <= 1 => {
+            let (rows, cols) = (*params.get("rows")?, *params.get("cols")?);
+            let MeasuredShape::Zipf { skew, keys } = shape else {
+                return None;
+            };
+            Some(topk_precision(skew, keys, rows, cols))
+        }
+        _ => None,
+    }
+}
+
+/// The share of the top [`TOPK_ENTRIES`] keys of Zipf(`skew`) over `keys`
+/// keys that a CMS with `rows` × `cols` keeps ranked, all at once with
+/// probability [`CONFIDENCE`] (see [`bound`]).
+fn topk_precision(skew: f64, keys: f64, rows: f64, cols: f64) -> f64 {
+    let k = TOPK_ENTRIES as usize;
+    if keys <= k as f64 {
+        return 1.0;
+    }
+    let total = zipf_normalizer(skew, keys);
+    let share = |rank: usize| (rank as f64).powf(-skew) / total;
+    let next = share(k + 1);
+    let mut failure = 0.0;
+    let mut ranked = 0;
+    for i in 1..=k {
+        failure += (keys - k as f64) * (cols * (share(i) - next)).powf(-rows);
+        if failure > 1.0 - CONFIDENCE {
+            break;
+        }
+        ranked += 1;
+    }
+    ranked as f64 / k as f64
+}
+
+/// `Σ_{n=1}^{K} n^−θ`: the first terms exactly, the rest by its integral
+/// with the endpoint correction (Euler–Maclaurin), so large `K` costs O(1).
+fn zipf_normalizer(skew: f64, keys: f64) -> f64 {
+    const EXACT: f64 = 1000.0;
+    let head = keys.min(EXACT) as usize;
+    let mut total: f64 = (1..=head).map(|n| (n as f64).powf(-skew)).sum();
+    if keys > EXACT {
+        let integral = |x: f64| {
+            if (skew - 1.0).abs() < 1e-12 {
+                x.ln()
+            } else {
+                x.powf(1.0 - skew) / (1.0 - skew)
+            }
+        };
+        total += integral(keys) - integral(EXACT) + (keys.powf(-skew) - EXACT.powf(-skew)) / 2.0;
+    }
+    total
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    fn params(pairs: &[(&str, f64)]) -> BTreeMap<String, f64> {
+        pairs.iter().map(|&(k, v)| (k.to_string(), v)).collect()
+    }
+
+    const PARETO: MeasuredShape = MeasuredShape::Pareto { tail_index: 2.0 };
+
+    #[test]
+    fn kll_hll_and_dd_bounds_ignore_n_and_merges() {
+        let kll = bound("kll-percall", &params(&[("k", 200.0)]), PARETO, 64).unwrap();
+        // 1.32% at 99%, scaled by sqrt(ln 40 / ln 200).
+        assert!((kll - 0.01109).abs() < 1e-4, "{kll}");
+        assert_eq!(
+            bound("dd", &params(&[("alpha", 0.01)]), PARETO, 1000),
+            Some(0.01)
+        );
+        let hll = bound("hll", &params(&[("lg_k", 12.0)]), PARETO, 4).unwrap();
+        assert!((hll - 1.96 * 1.04 / 64.0).abs() < 1e-12);
+    }
+
+    #[test]
+    fn topk_precision_holds_for_every_key_at_once() {
+        let zipf = MeasuredShape::Zipf {
+            skew: 1.1,
+            keys: 1e4,
+        };
+        let topk = |rows: f64, cols: f64, merges| {
+            bound(
+                "cms-heap-topk-fastpath-vector2d",
+                &params(&[("rows", rows), ("cols", cols)]),
+                zipf,
+                merges,
+            )
+        };
+        let wide = topk(5.0, 16384.0, 1).unwrap();
+        let narrow = topk(3.0, 256.0, 1).unwrap();
+        // More rows and columns keep more of the top 32 ranked; the 32nd and
+        // 33rd keys of Zipf 1.1 are too close for any of them.
+        assert!(wide > narrow, "{wide} {narrow}");
+        assert!(wide < 1.0 && wide > 0.0, "{wide}");
+        // Merged heaps: no guarantee.
+        assert_eq!(topk(5.0, 16384.0, 4), None);
+    }
+
+    #[test]
+    fn the_normalizer_matches_the_exact_sum() {
+        for (skew, keys) in [(1.1, 1e4), (0.5, 5e4), (1.0, 2e3), (2.0, 1e5)] {
+            let exact: f64 = (1..=keys as usize).map(|n| (n as f64).powf(-skew)).sum();
+            let approx = zipf_normalizer(skew, keys);
+            assert!(
+                (approx - exact).abs() / exact < 1e-6,
+                "{skew} {keys}: {approx} {exact}"
+            );
+        }
+    }
+
+    #[test]
+    fn sketches_without_a_published_bound_have_none() {
+        assert_eq!(bound("univmon-topk", &BTreeMap::new(), PARETO, 1), None);
+    }
+}