Closed form approximations

crates/geo_filters/Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "geo_filters"
-version = "0.1.1"
+version = "0.2.0"
 edition = "2021"
 description = "Geometric filters for set cardinality estimation."
 repository = "https://github.com/github/rust-gems"

crates/geo_filters/README.md

@@ -161,6 +161,47 @@ estimate), whereas a `near` candidate is below the bound and therefore scanned i
 and exact costs match. Since a nearest-neighbor scan rejects far more candidates than it keeps, the
 early abandon dominates the search cost. Reproduce with `cargo bench --bench nearest_neighbor`.
 
+### Precision and configuration choice
+
+Ranking candidates by their one-bit distance is a *noisy* estimate of the true symmetric-difference
+size (calibrating it to an item count with `Metric::to_f32` is order-preserving, so it does not
+change the ranking). The **relative error** below is the standard deviation of `estimate / true − 1`
+— i.e. relative to the true number of differing items. It is smallest for small differences and grows
+**slowly, like √(ln n)**, as the difference size `n` increases: the dense low bits saturate to a
+random ½ (encoding only the *parity* of many items, so they add variance without carrying size
+information). Modelling the buckets as independent Bernoulli variables, the count's variance is
+exactly `½·expected_diff_buckets(2n)`, which is logarithmic in `n`; the resulting relative error is
+about `√((γ + ln(4·(1−φ)·n)) / (2·S))` with `S = 2^b / (2·ln 2)` and `γ` the Euler–Mascheroni
+constant. For `GeoDiffConfig7` that formula gives ~8.6% at `n = 100`, ~18% at `10⁴`, ~24% at `10⁶`,
+close to the measured ~8% / ~18% / ~27% (real-filter bucket correlations push it a little higher at
+the top end).
+
+| configuration     | relative error (1σ), small … 1M-item difference | may reorder candidates within |
+| ----------------- | ----------------------------------------------- | ----------------------------- |
+| `GeoDiffConfig7`  | ~8% … ~27%                                      | ~2.4×                         |
+| `GeoDiffConfig10` | ~2.5% … ~6%                                     | ~1.3×                         |
+| `GeoDiffConfig13` | ~0.8% … ~2%                                     | ~1.1×                         |
+
+The one-bit count is *not* exact even for a few items — two items can hash to the same bucket — but
+such collisions are rare while the difference is small, so it is a low-variance estimate there.
+Compared to the standard windowed estimator `GeoDiffCount::size()` (and `size_with_sketch`), the
+one-bit count matches it for small differences but is noisier for large ones: for `GeoDiffConfig7` at
+a one-million-item difference it is ~27% here versus ~14% for `size()`, which reads only a fixed
+window near the fringe and ignores the saturated bits. Use `size_with_sketch` when you need an
+accurate absolute difference size; the one-bit metric exists for fast *ranking* of near neighbors,
+where its early-abandon speed matters more than the last few percent of precision.
+
+Because of the coarse resolution, **the metric is best used with `GeoDiffConfig10` or
+`GeoDiffConfig13` (`b ≥ 10`)**. `GeoDiffConfig7` reorders candidates whose true distances lie within
+roughly a factor of two of each other, so over a large candidate set it returns an *approximate*
+nearest neighbor rather than the exact one — though it never promotes a genuinely far candidate over
+a close one, since the noise cannot bridge more than the factor in the last column, even across
+hundreds of millions of candidates.
+
+For an exact result, use the metric as a cheap first pass to shortlist the top *k* candidates, then
+re-rank the shortlist with the more robust `Count::size_with_sketch`, which does not accumulate this
+variance.
+
 ## Evaluation
 
 Accuracy and performance evaluations for the predefined filter configurations are included in the repository:

crates/geo_filters/src/diff_count/config.rs

@@ -103,6 +103,79 @@ fn expected_diff_buckets_fast(phi: f64, items: f64) -> (f64, f64) {
     (sum, derivative)
 }
 
+/// Euler–Mascheroni constant.
+const GAMMA: f64 = 0.577_215_664_901_532_9;
+/// `e^gamma`, precomputed.
+const E_GAMMA: f64 = 1.781_072_417_990_198;
+
+// Closed-form approximation of [`expected_diff_buckets`] and its inverse that avoids the
+// per-bucket summation (and the Newton iteration on it). It is exact in both limits and stays
+// within ~0.35% of the exact model for the predefined configurations (b = 7, 10, 13).
+//
+// Writing `p_k = (1 - phi) phi^k`, the expected number of one-bits is
+//     E(N) = 1/2 sum_k (1 - (1 - 2 p_k)^N) ≈ S · Ein(x) + 1/4 (1 - e^-x),
+// with `x = 2(1 - phi) N`, `S = 1 / (2 ln(1/phi))`, and `Ein(x) = gamma + ln x + E1(x)` the
+// entire exponential-integral bridge (`Ein(x) ≈ x` near 0, `≈ gamma + ln x` for large x). The
+// special function `Ein` is replaced by a rational "scaling factor" `rho(x)` inside the logarithm,
+//     Ein(x) ≈ ln(1 + e^gamma x rho(x)),  rho(x) = (e^-gamma + A x + B x^2) / (1 + C x + B x^2),
+// where `rho(0) = e^-gamma` reproduces the linear small-N regime and `rho(inf) = 1` reproduces the
+// logarithmic asymptote. The coefficients are universal (independent of phi), fitted so the
+// rational matches `Ein` to < 0.07%.
+const RHO_A: f64 = 0.384_417;
+const RHO_B: f64 = 0.130_468;
+const RHO_C: f64 = 0.442_202;
+
+/// Closed-form approximation of [`expected_diff_buckets`]: the expected number of one-bits for
+/// `items` items, together with its derivative with respect to `items`. See the comment above for
+/// the derivation. Accurate to < 1% for the predefined configurations.
+pub(super) fn expected_diff_buckets_approx(phi: f64, items: f64) -> (f64, f64) {
+    if items <= 0.0 {
+        return (0.0, 1.0);
+    }
+    let s = 0.5 / -phi.ln(); // = 1 / (2 ln(1/phi))
+    let a = 2.0 * (1.0 - phi);
+    let x = a * items;
+    let e_neg_x = (-x).exp();
+
+    // Y = e^gamma x rho(x) = (x + e^gamma A x^2 + e^gamma B x^3) / (1 + C x + B x^2)
+    let num = x * (1.0 + E_GAMMA * RHO_A * x + E_GAMMA * RHO_B * x * x);
+    let den = 1.0 + RHO_C * x + RHO_B * x * x;
+    let y = num / den;
+    let value = s * (1.0 + y).ln() + 0.25 * (1.0 - e_neg_x);
+
+    // derivative w.r.t. items = a * d/dx [ s ln(1 + Y) + 1/4 (1 - e^-x) ]
+    let num_prime = 1.0 + 2.0 * E_GAMMA * RHO_A * x + 3.0 * E_GAMMA * RHO_B * x * x;
+    let den_prime = RHO_C + 2.0 * RHO_B * x;
+    let y_prime = (num_prime * den - num * den_prime) / (den * den);
+    let derivative = a * (s * y_prime / (1.0 + y) + 0.25 * e_neg_x);
+
+    (value, derivative)
+}
+
+/// Closed-form approximation of the inverse of [`expected_diff_buckets`]: the number of items that
+/// produce `ones` one-bits. Uses the analytic large-N seed of the model above followed by three
+/// Newton steps on the (cheap) closed-form forward, which converge to full `f64` precision.
+pub(super) fn estimate_count_approx(phi: f64, ones: f64) -> f64 {
+    if ones <= 0.0 {
+        return 0.0;
+    }
+    let s = 0.5 / -phi.ln();
+    let a = 2.0 * (1.0 - phi);
+    // Invert the large-N asymptote `m ≈ S (gamma + ln(a N)) + 1/4` for the initial guess.
+    let mut items = ((ones - 0.25) / s - GAMMA).exp() / a;
+    if !items.is_finite() || items <= 0.0 {
+        items = ones;
+    }
+    for _ in 0..3 {
+        let (m, dm) = expected_diff_buckets_approx(phi, items);
+        items -= (m - ones) / dm;
+        if items < f64::MIN_POSITIVE {
+            items = f64::MIN_POSITIVE;
+        }
+    }
+    items
+}
+
 static LOOKUPS: [Lazy<Lookup>; 16] = [
     Lazy::new(|| build_lookup(0)),
     Lazy::new(|| build_lookup(1)),
@@ -180,6 +253,52 @@ mod tests {
         println!("{:?}", estimate_count(phi, 78.450, expected_diff_buckets));
     }
 
+    #[test]
+    fn test_expected_diff_buckets_approx() {
+        // The closed-form forward stays within 1% of the exact model for the predefined configs.
+        for b in [7u32, 10, 13] {
+            let phi = 0.5f64.powf(1.0 / (1u64 << b) as f64);
+            let mut items = 1.0f64;
+            loop {
+                let exact = expected_diff_buckets(phi, items).0;
+                if exact > 2000.0 {
+                    break;
+                }
+                if exact > 0.3 {
+                    let approx = expected_diff_buckets_approx(phi, items).0;
+                    let rel = (approx - exact).abs() / exact;
+                    assert!(
+                        rel < 0.01,
+                        "b={b} items={items}: exact={exact} approx={approx} rel={rel}"
+                    );
+                }
+                items *= 1.5;
+            }
+        }
+    }
+
+    #[test]
+    fn test_estimate_count_approx() {
+        // The closed-form inverse recovers the item count within 1% of the exact model.
+        for b in [7u32, 10, 13] {
+            let phi = 0.5f64.powf(1.0 / (1u64 << b) as f64);
+            let mut items = 1.0f64;
+            loop {
+                let ones = expected_diff_buckets(phi, items).0;
+                if ones > 2000.0 {
+                    break;
+                }
+                let est = estimate_count_approx(phi, ones);
+                let rel = (est - items).abs() / items;
+                assert!(
+                    rel < 0.01,
+                    "b={b} items={items} ones={ones}: est={est} rel={rel}"
+                );
+                items *= 1.5;
+            }
+        }
+    }
+
     #[test]
     fn test_estimation_lut_7() {
         let c = GeoDiffConfig7::<UnstableDefaultBuildHasher>::default();

crates/geo_filters/src/diff_count/metric.rs

@@ -5,10 +5,8 @@ use std::cmp::Ordering;
 use std::fmt;
 use std::ops::Add;
 
-use super::config::expected_diff_buckets;
-use crate::config::{
-    estimate_count, xor_bit_chunks, BitChunk, GeoConfig, IsBucketType, BITS_PER_BLOCK,
-};
+use super::config::{estimate_count_approx, expected_diff_buckets_approx};
+use crate::config::{xor_bit_chunks, BitChunk, GeoConfig, IsBucketType, BITS_PER_BLOCK};
 use crate::diff_count::GeoDiffCount;
 use crate::{Diff, Metric, MetricSpace};
 
@@ -48,24 +46,23 @@ impl<C: GeoConfig<Diff> + Default> Metric for OnesMetric<C> {
     }
 
     /// The calibrated size (item count) that this number of one-bits represents, obtained by
-    /// inverting the expected-buckets function with the same Newton iteration used to build the
-    /// estimation lookup table.
+    /// inverting the closed-form approximation of the expected-buckets function.
     fn to_f32(&self) -> f32 {
         if *self == Self::infinite() {
             return f32::INFINITY;
         }
         let phi = C::default().phi_f64();
-        estimate_count(phi, self.get() as f64, expected_diff_buckets).0 as f32
+        estimate_count_approx(phi, self.get() as f64) as f32
     }
 
     /// The number of one-bits expected for the given calibrated size, obtained by evaluating the
-    /// forward expected-buckets function.
+    /// closed-form approximation of the forward expected-buckets function.
     fn from_f32(value: f32) -> Self {
         if !value.is_finite() {
             return Self::infinite();
         }
         let phi = C::default().phi_f64();
-        let buckets = expected_diff_buckets(phi, value.max(0.0) as f64).0;
+        let buckets = expected_diff_buckets_approx(phi, value.max(0.0) as f64).0;
         Self::new(buckets.round() as usize)
     }
 }
@@ -121,6 +118,12 @@ impl<C: GeoConfig<Diff>> fmt::Debug for OnesMetric<C> {
 /// search - and, given a bound, abandoned early. The cached count (see [`Self::size`]) also yields
 /// an O(1) reverse-triangle lower bound on the distance. Once the nearest neighbor is found,
 /// [`Self::filter`] gives access to the underlying filter for a calibrated size estimate.
+///
+/// The one-bit distance is a noisy estimate of the true difference size, so this metric is best
+/// used with `GeoDiffConfig10` or `GeoDiffConfig13` (`b >= 10`). The coarser `GeoDiffConfig7` can
+/// reorder candidates whose true distances are within roughly a factor of two of each other,
+/// returning an approximate rather than exact nearest neighbor; for an exact result, shortlist with
+/// this metric and re-rank the top candidates with [`Count::size_with_sketch`](crate::Count).
 pub struct GeoDiffMetric<'a, C: GeoConfig<Diff>> {
     filter: GeoDiffCount<'a, C>,
     ones: OnesMetric<C>,