Add HLL++ bias correction to the dense visitor estimate

packages/core/src/hll-bias.ts

@@ -0,0 +1,94 @@
+/**
+ * HyperLogLog++ bias correction for precision 11 (2,048 registers).
+ *
+ * Source: Heule, Nunkesser, Hall, "HyperLogLog in Practice" (EDBT 2013), the
+ * appendix data for p = 11: `rawEstimateData` (RAW_ESTIMATES), `biasData`
+ * (BIASES) and the linear counting threshold. Entry i is the mean raw estimate
+ * the authors measured at one cardinality and the mean bias of that estimate;
+ * raw − bias gives the cardinality, from 1 to 5·m = 10,240 in steps of 51 or 52.
+ */
+
+/** The paper's Threshold(11): linear counting is used while its estimate is at most this. */
+export const LINEAR_COUNTING_THRESHOLD = 1800;
+
+const NEIGHBOURS = 6;
+
+/** Mean bias of the 6 table points nearest to `raw` (nearest-neighbour interpolation, k = 6). */
+export function estimateBias(raw: number): number {
+  let lo = 0;
+  let hi = RAW_ESTIMATES.length;
+  while (lo < hi) {
+    const mid = (lo + hi) >>> 1;
+    if (RAW_ESTIMATES[mid]! < raw) lo = mid + 1;
+    else hi = mid;
+  }
+  let left = lo - 1;
+  let right = lo;
+  let sum = 0;
+  for (let k = 0; k < NEIGHBOURS; k++) {
+    const takeLeft =
+      right === RAW_ESTIMATES.length ||
+      (left >= 0 && raw - RAW_ESTIMATES[left]! <= RAW_ESTIMATES[right]! - raw);
+    sum += takeLeft ? BIASES[left--]! : BIASES[right++]!;
+  }
+  return sum / NEIGHBOURS;
+}
+
+export const RAW_ESTIMATES: readonly number[] = [
+  1477, 1501.6014, 1526.5802, 1551.7942, 1577.3042, 1603.2062, 1629.8402, 1656.2292, 1682.9462,
+  1709.9926, 1737.3026, 1765.4252, 1793.0578, 1821.6092, 1849.626, 1878.5568, 1908.527, 1937.5154,
+  1967.1874, 1997.3878, 2027.37, 2058.1972, 2089.5728, 2120.1012, 2151.9668, 2183.292, 2216.0772,
+  2247.8578, 2280.6562, 2313.041, 2345.714, 2380.3112, 2414.1806, 2447.9854, 2481.656, 2516.346,
+  2551.5154, 2586.8378, 2621.7448, 2656.6722, 2693.5722, 2729.1462, 2765.4124, 2802.8728, 2838.898,
+  2876.408, 2913.4926, 2951.4938, 2989.6776, 3026.282, 3065.7704, 3104.1012, 3143.7388, 3181.6876,
+  3221.1872, 3261.5048, 3300.0214, 3339.806, 3381.409, 3421.4144, 3461.4294, 3502.2286, 3544.651,
+  3586.6156, 3627.337, 3670.083, 3711.1538, 3753.5094, 3797.01, 3838.6686, 3882.1678, 3922.8116,
+  3967.9978, 4009.9204, 4054.3286, 4097.5706, 4140.6014, 4185.544, 4229.5976, 4274.583, 4316.9438,
+  4361.672, 4406.2786, 4451.8628, 4496.1834, 4543.505, 4589.1816, 4632.5188, 4678.2294, 4724.8908,
+  4769.0194, 4817.052, 4861.4588, 4910.1596, 4956.4344, 5002.5238, 5048.13, 5093.6374, 5142.8162,
+  5187.7894, 5237.3984, 5285.6078, 5331.0858, 5379.1036, 5428.6258, 5474.6018, 5522.7618, 5571.5822,
+  5618.59, 5667.9992, 5714.88, 5763.454, 5808.6982, 5860.3644, 5910.2914, 5953.571, 6005.9232,
+  6055.1914, 6104.5882, 6154.5702, 6199.7036, 6251.1764, 6298.7596, 6350.0302, 6398.061, 6448.4694,
+  6495.933, 6548.0474, 6597.7166, 6646.9416, 6695.9208, 6742.6328, 6793.5276, 6842.1934, 6894.2372,
+  6945.3864, 6996.9228, 7044.2372, 7094.1374, 7142.2272, 7192.2942, 7238.8338, 7288.9006, 7344.0908,
+  7394.8544, 7443.5176, 7490.4148, 7542.9314, 7595.6738, 7641.9878, 7694.3688, 7743.0448, 7797.522,
+  7845.53, 7899.594, 7950.3132, 7996.455, 8050.9442, 8092.9114, 8153.1374, 8197.4472, 8252.8278,
+  8301.8728, 8348.6776, 8401.4698, 8453.551, 8504.6598, 8553.8944, 8604.1276, 8657.6514, 8710.3062,
+  8758.908, 8807.8706, 8862.1702, 8910.4668, 8960.77, 9007.2766, 9063.164, 9121.0534, 9164.1354,
+  9218.1594, 9267.767, 9319.0594, 9372.155, 9419.7126, 9474.3722, 9520.1338, 9572.368, 9622.7702,
+  9675.8448, 9726.5396, 9778.7378, 9827.6554, 9878.1922, 9928.7782, 9978.3984, 10026.578,
+  10076.5626, 10137.1618, 10177.5244, 10229.9176,
+];
+
+export const BIASES: readonly number[] = [
+  1476, 1449.6014, 1423.5802, 1397.7942, 1372.3042, 1347.2062, 1321.8402, 1297.2292, 1272.9462,
+  1248.9926, 1225.3026, 1201.4252, 1178.0578, 1155.6092, 1132.626, 1110.5568, 1088.527, 1066.5154,
+  1045.1874, 1024.3878, 1003.37, 982.1972, 962.5728, 942.1012, 922.9668, 903.292, 884.0772,
+  864.8578, 846.6562, 828.041, 809.714, 792.3112, 775.1806, 757.9854, 740.656, 724.346, 707.5154,
+  691.8378, 675.7448, 659.6722, 645.5722, 630.1462, 614.4124, 600.8728, 585.898, 572.408, 558.4926,
+  544.4938, 531.6776, 517.282, 505.7704, 493.1012, 480.7388, 467.6876, 456.1872, 445.5048, 433.0214,
+  420.806, 411.409, 400.4144, 389.4294, 379.2286, 369.651, 360.6156, 350.337, 342.083, 332.1538,
+  322.5094, 315.01, 305.6686, 298.1678, 287.8116, 280.9978, 271.9204, 265.3286, 257.5706, 249.6014,
+  242.544, 235.5976, 229.583, 220.9438, 214.672, 208.2786, 201.8628, 195.1834, 191.505, 186.1816,
+  178.5188, 172.2294, 167.8908, 161.0194, 158.052, 151.4588, 148.1596, 143.4344, 138.5238, 133.13,
+  127.6374, 124.8162, 118.7894, 117.3984, 114.6078, 109.0858, 105.1036, 103.6258, 98.6018000000004,
+  95.7618000000002, 93.5821999999998, 88.5900000000001, 86.9992000000002, 82.8800000000001,
+  80.4539999999997, 74.6981999999998, 74.3644000000004, 73.2914000000001, 65.5709999999999,
+  66.9232000000002, 65.1913999999997, 62.5882000000001, 61.5702000000001, 55.7035999999998,
+  56.1764000000003, 52.7596000000003, 53.0302000000001, 49.0609999999997, 48.4694, 44.933,
+  46.0474000000004, 44.7165999999997, 41.9416000000001, 39.9207999999999, 35.6328000000003,
+  35.5276000000003, 33.1934000000001, 33.2371999999996, 33.3864000000003, 33.9228000000003,
+  30.2371999999996, 29.1373999999996, 25.2272000000003, 24.2942000000003, 19.8338000000003,
+  18.9005999999999, 23.0907999999999, 21.8544000000002, 19.5176000000001, 15.4147999999996,
+  16.9314000000004, 18.6737999999996, 12.9877999999999, 14.3688000000002, 12.0447999999997,
+  15.5219999999999, 12.5299999999997, 14.5940000000001, 14.3131999999996, 9.45499999999993,
+  12.9441999999999, 3.91139999999996, 13.1373999999996, 5.44720000000052, 9.82779999999912,
+  7.87279999999919, 3.67760000000089, 5.46980000000076, 5.55099999999948, 5.65979999999945,
+  3.89439999999922, 3.1275999999998, 5.65140000000065, 6.3062000000009, 3.90799999999945,
+  1.87060000000019, 5.17020000000048, 2.46680000000015, 0.770000000000437, -3.72340000000077,
+  1.16400000000067, 8.05340000000069, 0.135399999999208, 2.15940000000046, 0.766999999999825,
+  1.0594000000001, 3.15500000000065, -0.287399999999252, 2.37219999999979, -2.86620000000039,
+  -1.63199999999961, -2.22979999999916, -0.15519999999924, -1.46039999999994, -0.262199999999211,
+  -2.34460000000036, -2.8078000000005, -3.22179999999935, -5.60159999999996, -8.42200000000048,
+  -9.43740000000071, 0.161799999999857, -10.4755999999998, -10.0823999999993,
+];

packages/core/src/hll.ts

@@ -1,6 +1,10 @@
 import { hash32 } from './hash.js';
+import { LINEAR_COUNTING_THRESHOLD, estimateBias } from './hll-bias.js';
 
-/** Register index bits. 2^11 registers → ~2.3% standard error, 2 KiB dense. */
+/**
+ * Register index bits. 2^11 registers → ~2.3% standard error, 2 KiB dense.
+ * The bias data in hll-bias.ts is for this precision only.
+ */
 export const HLL_PRECISION = 11;
 const M = 1 << HLL_PRECISION;
 /** Sparse sets switch to dense registers once they would be larger than them. */
@@ -72,9 +76,13 @@ export class Hll {
       sum += POW2_NEG[r]!;
       if (r === 0) zeros++;
     }
-    let estimate = (ALPHA * M * M) / sum;
-    if (estimate <= 2.5 * M && zeros > 0) estimate = M * Math.log(M / zeros);
-    return Math.round(estimate);
+    if (zeros > 0) {
+      const linear = M * Math.log(M / zeros);
+      if (linear <= LINEAR_COUNTING_THRESHOLD) return Math.round(linear);
+    }
+    // HLL++: the raw estimate overshoots up to 5·m; subtract the paper's empirical bias there.
+    const raw = (ALPHA * M * M) / sum;
+    return Math.round(raw <= 5 * M ? raw - estimateBias(raw) : raw);
   }
 
   clone(): Hll {

packages/core/test/fixtures/stored-sketches.json

@@ -0,0 +1,11 @@
+{
+  "description": "Hll.toBytes() output written at commit 5239238, before the bias correction. The ids are stored-0, stored-1, ...",
+  "sparse": {
+    "ids": 20,
+    "hex": "0053d11ed4eaee97de2eded8d5ceeea3d329b23af50c7e3045169195d0a1d37ec11fffdc212f44b4bec7861012e206779a442cb0753a244deec44d397b011277c5834de891bc8361f56f3ae2e2f39cd024"
+  },
+  "dense": {
+    "ids": 5000,
+    "hex": "01030101010205020403040204020800020204010103040104000201000301010202010206020502010102060602050504030203000502040203000402030301040304020400050103010604010602010800030402030303020101040302020102010402010002020104050102020103040102010403010103020b0603020103020403000300080303030501070402020308000103020501050306020203080206030305030102010202020405030501010202010102020503040202040303000604030104010203020401030003000204020200030402020104040601020701040302070303040304050002040301020201030401050203030301020005060302030103050104010500010202040402030102020404000304060902030206020203040503040302040002020207020404010404020302010002010302000105000303030306040500010300040301040202010204070203030402020202040202010001070302020603000004020302000405030101050306010301020304020102010200010003000601010202000302040402040204040500060501010101020502040304020301030101060201030300000802050301020402010400040103030a0103010204020103050101030203020008020003020208030402020200070501010304060300010100070204010002060304020403050102020303020703020505010203000301030402030403080101050402020002080104040200010501020102010302050003060102000002070303030200020201020201040201010103020203040102050504090101020502010204030102030602010405060202030400020105010102040201020103010303020105040103010202020202000404000103050301030104020101060505040301000101030304040204010102040203010301000203020002020701020404080702060400030306050607020202020203030102010303020508040404050406010501020402040306070103020200020502080205010503060102040804000400010603020104020403020703020405020302020202030304080204020600030102010302020302040401030500030103010a01040103020604010304010005000101020405030403070200040301010202020303030202030002010201050401020500050101020104010203040309040202020203030501030202080401020101010202030202040302070103030005030400020101010302040302010303040102010402020705030302050601000203000202000301010503040302020b030102030103040003020703080300040505020704020102060002050403010004010001010002020304050001000104030500010001020406020103000303050403010801010103030002040202010400030200020205010003010401030804030401010802030401020202040502020104020502040001010102010100060101020102030002040202000308020605020302060001030100000205040205030306080403010503020302020201010203010200020000010205020203030102010101020404010401030200010303030002050201000603020002030002020005010203000502040101020102010309000b04030206040203020202010104050504020102030001010102050304010307000403010300040301010002020202010700010203030503010500060103010502010402010604050306050503010101010301030103030401020301040001000203020302020302040205030103020403000107010402040202020300000002020301000303010603030500030304030402020401030903020407020401040506010201040100050502040204020302040202020304030101010102030101030301030602010501040401010801030204030104050103030203040604020402030400030104030206020101020102040502030206010100000103000102010503000107010504040101010103010103040403040705030603040201010302020402060004020202010101010401040305030304040202070302010401010201010303030000020202030002070103000003020205000203010104020202050202030503010005030300020103010202040301020105030102050203030401020108000401020104010102010202010102050002020202030302010202030003030505010304030304030602050301050306020802040301010004020300010403050102030201040404050501010202020403040502010002050102040301030201030100020203030306030102050004050102030203010101040404020106010204050508000403070203020404020303020200010703010204010304040301030103060203050303020401010304010203030701030301020202030503030004000403010202040b0002010606030503010201010101000202020505020302000704030101010202030103000303000401020404010203060201010100010302020304020301020200010205030301010202030505000205020705030102070203010601010302030304030001040300010203030101040201080d030106040a0003000402010204030206000302000102000400010502030200010004010805010403040402030201030302030203000300040204010003010202010600000308020302030905050204030602030500020601020201020302030101020206020103020200010200050101010104010200010301060103000101010206000302010200000202040200000203010505000101020102020402020204"
+  }
+}

packages/core/test/hll-bias.test.ts

@@ -0,0 +1,28 @@
+import { describe, expect, it } from 'vitest';
+import { BIASES, LINEAR_COUNTING_THRESHOLD, RAW_ESTIMATES, estimateBias } from '../src/hll-bias.js';
+
+describe('HLL++ bias table (p = 11)', () => {
+  it('pairs each mean raw estimate with its bias at cardinalities from 1 to 5·m', () => {
+    expect(RAW_ESTIMATES).toHaveLength(201);
+    expect(BIASES).toHaveLength(201);
+    const cardinalities = RAW_ESTIMATES.map((raw, i) => raw - BIASES[i]!);
+    expect(cardinalities[0]).toBeCloseTo(1, 6);
+    expect(cardinalities[200]).toBeCloseTo(5 * 2048, 6);
+    for (let i = 1; i < RAW_ESTIMATES.length; i++) {
+      expect(RAW_ESTIMATES[i]!).toBeGreaterThan(RAW_ESTIMATES[i - 1]!);
+      expect(Math.round(cardinalities[i]! - cardinalities[i - 1]!)).toBeOneOf([51, 52]);
+    }
+    expect(LINEAR_COUNTING_THRESHOLD).toBe(1800);
+  });
+
+  it('averages the bias of the 6 nearest table points', () => {
+    const nearestSix = (raw: number) =>
+      RAW_ESTIMATES.map((point, i) => ({ distance: Math.abs(point - raw), bias: BIASES[i]! }))
+        .sort((a, b) => a.distance - b.distance)
+        .slice(0, 6)
+        .reduce((sum, p) => sum + p.bias / 6, 0);
+    const probes = [1_000, 20_000, ...RAW_ESTIMATES];
+    for (let raw = 1_400; raw < 10_400; raw += 13.7) probes.push(raw);
+    for (const raw of probes) expect(estimateBias(raw)).toBeCloseTo(nearestSix(raw), 9);
+  });
+});

packages/core/test/hll.test.ts

@@ -1,8 +1,37 @@
 import { describe, expect, it } from 'vitest';
 import { Hll, HLL_PRECISION, hash32 } from '../src/index.js';
+import stored from './fixtures/stored-sketches.json';
 
 const ids = (n: number, prefix = 'v') => Array.from({ length: n }, (_, i) => `${prefix}${i}`);
 
+/** 1.04 / √m, the textbook standard error at this precision (2.3%). */
+const STANDARD_ERROR = 1.04 / Math.sqrt(1 << HLL_PRECISION);
+const SEED = 2013;
+
+/** mulberry32: uniform 32-bit hashes, the same sequence on every run for a given seed. */
+function hashStream(seed: number): () => number {
+  let a = seed >>> 0;
+  return () => {
+    a = (a + 0x6d2b79f5) | 0;
+    let t = Math.imul(a ^ (a >>> 15), 1 | a);
+    t = (t + Math.imul(t ^ (t >>> 7), 61 | t)) ^ t;
+    return (t ^ (t >>> 14)) >>> 0;
+  };
+}
+
+/** Registers of a dense sketch or the sorted hashes of a sparse one. */
+function canonical(h: Hll): string {
+  const bytes = h.toBytes();
+  if (!h.isSparse) return bytes.join(',');
+  const view = new DataView(bytes.buffer);
+  const hashes: number[] = [];
+  for (let offset = 1; offset < bytes.length; offset += 4)
+    hashes.push(view.getUint32(offset, true));
+  return `sparse:${hashes.sort((a, b) => a - b).join(',')}`;
+}
+
+const fromHex = (hex: string) => Uint8Array.from(hex.match(/../g)!, (byte) => parseInt(byte, 16));
+
 describe('hash32', () => {
   it('is deterministic and seed-sensitive', () => {
     expect(hash32('visitor-1')).toBe(hash32('visitor-1'));
@@ -76,4 +105,110 @@ describe('Hll', () => {
     expect(() => Hll.fromBytes(new Uint8Array([1, 2, 3]))).toThrow(RangeError);
     expect(() => Hll.fromBytes(new Uint8Array([0, 1, 2]))).toThrow(RangeError);
   });
+
+  it('decodes sketches stored before the bias correction and writes the same bytes', () => {
+    for (const { ids: n, hex } of [stored.sparse, stored.dense]) {
+      const bytes = fromHex(hex);
+      const fresh = Hll.from(ids(n, 'stored-'));
+      const decoded = Hll.fromBytes(bytes);
+      expect(fresh.toBytes()).toEqual(bytes);
+      expect(decoded.toBytes()).toEqual(bytes);
+      expect(decoded.count()).toBe(fresh.count());
+    }
+    expect(Hll.fromBytes(fromHex(stored.sparse.hex)).count()).toBe(20);
+    const merged = Hll.fromBytes(fromHex(stored.dense.hex)).merge(Hll.from(ids(6_000, 'stored-')));
+    expect(canonical(merged)).toBe(canonical(Hll.from(ids(6_000, 'stored-'))));
+  });
+
+  it('merges associatively and commutatively across sparse and dense sketches', () => {
+    const inputs = [
+      ids(200, 'a'), // sparse
+      ids(250, 'b'), // sparse, still sparse when merged with the first
+      ids(300, 'a'), // sparse, overlaps the first; all three together are dense
+      ids(4_000, 'd'), // dense
+    ];
+    const sketch = (i: number) => Hll.from(inputs[i]!);
+    const union = (...picked: number[]) => Hll.from(picked.flatMap((i) => inputs[i]!));
+    for (let x = 0; x < inputs.length; x++) {
+      for (let y = 0; y < inputs.length; y++) {
+        expect(canonical(sketch(x).merge(sketch(y)))).toBe(canonical(union(x, y)));
+        for (let z = 0; z < inputs.length; z++) {
+          const left = sketch(x).merge(sketch(y)).merge(sketch(z));
+          const right = sketch(x).merge(sketch(y).merge(sketch(z)));
+          expect(canonical(left)).toBe(canonical(right));
+          expect(canonical(left)).toBe(canonical(union(x, y, z)));
+          expect(left.count()).toBe(sketch(z).merge(sketch(y)).merge(sketch(x)).count());
+        }
+      }
+    }
+    expect(union(0, 1).isSparse).toBe(true);
+    expect(union(0, 1, 2).isSparse).toBe(false);
+  });
+});
+
+describe('Hll accuracy (seeded hash streams)', () => {
+  it('is exact up to 512 and within 3% rms and 1% mean error from 513 to 200,000', () => {
+    const sizes = [
+      1, 2, 10, 100, 512, 513, 1_000, 1_500, 1_800, 2_000, 3_000, 4_000, 5_000, 6_000, 7_000, 8_000,
+      10_000, 12_000, 20_000, 50_000, 100_000, 200_000,
+    ];
+    const trials = 100;
+    const next = hashStream(SEED);
+    const errors = sizes.map(() => [] as number[]);
+    for (let t = 0; t < trials; t++) {
+      const h = new Hll();
+      let n = 0;
+      sizes.forEach((size, s) => {
+        for (; n < size; n++) h.addHash(next());
+        errors[s]!.push((h.count() - size) / size);
+      });
+    }
+    const outOfBound = sizes.flatMap((n, s) => {
+      const e = errors[s]!;
+      const rms = Math.sqrt(e.reduce((acc, x) => acc + x * x, 0) / trials);
+      const mean = e.reduce((acc, x) => acc + x, 0) / trials;
+      const ok = n <= 512 ? e.every((x) => x === 0) : rms < 0.03 && Math.abs(mean) < 0.01;
+      return ok ? [] : [{ n, rms: rms.toFixed(4), mean: mean.toFixed(4) }];
+    });
+    expect(outOfBound).toEqual([]);
+  });
+
+  it('does not jump where one estimator hands over to another', () => {
+    // Linear counting hands over near 1,800 and the bias table ends at a raw estimate of
+    // 5·m = 10,240; before the correction, the raw estimate took over near 5,120.
+    const every = 5;
+    const windows = [1_600, 4_600, 9_800].map((lo) =>
+      Array.from({ length: 201 }, (_, i) => lo + i * every),
+    );
+    const trials = 100;
+    const next = hashStream(SEED);
+    const sums = windows.map((points) => points.map(() => 0));
+    let largestChange = 0;
+    for (let t = 0; t < trials; t++) {
+      const h = new Hll();
+      let n = 0;
+      windows.forEach((points, w) => {
+        let previous: number | undefined;
+        points.forEach((point, i) => {
+          for (; n < point; n++) h.addHash(next());
+          const count = h.count();
+          if (previous !== undefined) {
+            largestChange = Math.max(largestChange, Math.abs(count - previous) / point);
+          }
+          previous = count;
+          sums[w]![i]! += count;
+        });
+      });
+    }
+    const worstMeanError = Math.max(
+      ...windows.flatMap((points, w) =>
+        points.map((point, i) => Math.abs(sums[w]![i]! / trials / point - 1)),
+      ),
+    );
+    // A systematic jump would push the mean of 100 sketches off n after the handover.
+    expect(worstMeanError).toBeLessThan(0.01);
+    // At a handover one sketch still switches between two estimates of about the same n, so
+    // it can step, but by less than two standard errors per 5 added hashes.
+    expect(largestChange).toBeLessThan(2 * STANDARD_ERROR);
+  });
 });