[MOD-17526][MOD-17527] Make SQ8 metadata exact and fix the symmetric L2 formulation

[dor-forer]

First of two stacked PRs closing MOD-17526. Found while reviewing #1007, which is what first makes any of this reachable: on main nothing constructs an index that uses QuantPreprocessor, so all of it is latent today.

Describe the changes in the pull request

SQ8 stored metadata now describes the vector the index actually holds, the reconstruction
min + delta * a[i], instead of the original input, and the symmetric L2 formulation no longer
cancels the answer away for vectors sharing a large offset. Two integer accumulators that could
overflow are widened, the plain uint8 L2 and scalar paths get the same signedness fix the inner
product side already received, and the quantizer no longer produces undefined behaviour for input
whose range is not representable in FP32. Blob size and slot count are unchanged.

What was wrong

1. Stored sums described the wrong vector. QuantPreprocessor accumulated x_sum and x_sum_squares over the original input, while every kernel that consumed them computed its other terms from the reconstruction min + delta * a[i]. Mixing two vectors into one formula gives wrong distances, and they can be negative: storing [0, 0.25, 1] and querying [0, 0.2501, 1] returned -0.000490427, so a radius-0 range query returned a vector that is not identical. This affected both the symmetric IP kernels and all the L2 kernels.

2. The symmetric L2 kernels cancelled away the answer. They used ||x||^2 + ||y||^2 - 2*IP in FP32. When the vectors share a large offset the two large terms nearly cancel and the result lives in bits rounding already discarded: storing [100000, 100008] against [100000, 100000] returned 0 where the truth is 64. Quantization is exact for that input, so this is the summation form, not the 8 bits.

3. The integer dot product was not exact. NEON and SVE accumulated exactly in integer lanes and then discarded it by returning float, exact only to dimension 258. AVX512 reduced 16 int32 lanes with _mm512_reduce_add_epi32, which overflows from dimension 33,027 (255*255*dim > INT_MAX). That is MOD-17527, folded in here because the L2 fix depends on the dot being exact, and doing it separately would mean writing the L2 expansion twice.

4. The plain uint8 kernels kept the same signed accumulator, on the paths MOD-17527 did not reach. L2_NEON_UINT8.h assigned vaddvq_u32 to an int32_t, L2_AVX512F_BW_VL_VNNI_UINT8.h returned _mm512_reduce_add_epi32 as a signed int, and ret_t for a 1-byte element type was int in both IP.cpp and L2.cpp, which is signed-overflow UB from dimension 33,026 while its own comment claimed support to 2^16. L2_NEON_DOTPROD_UINT8.h and L2_SVE_UINT8.h were already unsigned, which is why the defect was easy to miss. Fixing IP but not L2 for the same data type was not defensible, so both are in scope here.

5. Self-distance was still not exactly zero after fixing 1 to 3. The six-term L2 cancels exactly in real arithmetic, but the compiler contracts some of the products into fused multiply-adds and not others, so they stop rounding identically and stop cancelling. The AVX512 kernel returned about -6e-14 at dimension 64 and -1.7e-12 at 512. Harmless to ranking, but a caller taking a square root gets NaN, and a function named L2Sqr should not return a negative number.

Which issues this PR fixes

1. MOD-17526: SQ8 L2 kernels lose the answer through the ||x||^2 + ||y||^2 - 2*IP identity in FP32, and stored sums that describe the input rather than the reconstruction.
2. MOD-17527: the SQ8 to SQ8 integer dot product overflows a signed 32-bit accumulator past dimension 33,026. The plain uint8 L2 and scalar paths carry the identical defect and are fixed here too; they have no ticket of their own, and splitting them out would leave one data type half fixed.
3. MOD-17528: QuantPreprocessor converts NaN to uint8 when the derived range is not representable, which is undefined behaviour reachable from AddVector for input the API accepts.

Main objects this PR modified

1. vecsim_types::sq8 (src/VecSim/types/sq8.h): metadata slots redefined as Q_SUM and Q_SUM_SQUARES, plus unaligned readers and the four reconstructed_* derivation helpers.
2. QuantPreprocessor (src/VecSim/spaces/computer/preprocessors.h): sums the quantized bytes, computes the range in double, and drops a libm call per element.
3. The SQ8 to SQ8, SQ8 to FP32 and SQ8 to FP16 distance kernels across AVX512, AVX2, SSE4, NEON, NEON_DOTPROD, SVE and SVE2, plus their scalar fallbacks in IP.cpp and L2.cpp.
4. The plain uint8 kernels: UINT8_InnerProductImp on all four ISAs, the L2 reduces in L2_AVX512F_BW_VL_VNNI_UINT8.h and L2_NEON_UINT8.h, and ret_t in IP.cpp and L2.cpp.
5. tests/unit/test_spaces.cpp, tests/unit/test_components.cpp and the blob-building helpers in tests/utils/tests_utils.h and tests/unit/unit_test_utils.h, plus the uint8 spaces benchmark fixture.

Mark if applicable

• This PR introduces API changes
• This PR introduces serialization changes

Neither box is checked, deliberately, and the second one deserves a sentence. The stored SQ8 blob keeps its exact size and slot offsets, but two of those slots change meaning, so a blob written by older code would be misinterpreted by this one. That is not a serialization change in practice because no code path on main constructs an index that uses QuantPreprocessor, so no persisted index can contain an SQ8 blob. SQ8 index creation arrives with #1007, which is stacked on this.

What this changes

sq8.h redefines two metadata slots as Q_SUM and Q_SUM_SQUARES, uint32 sums over the quantized bytes. Blob size and slot count are unchanged: for an L2 index the SUM slot was dead weight, since only the symmetric IP kernels read it. The header gains unaligned readers and four exact-derivation helpers so the layout has exactly one interpretation, in the spirit of storage_bytes_count.

The symmetric L2 kernels now use the expanded difference, which keeps any common offset in a (min1 - min2) term formed once and exactly, and combine in double. Measured relative error is ~2e-8, including at dim 16384 where the FP32 combination was 16% off. Two alternatives were measured and rejected: combining the exact sums in FP32 (16% error), and widening only the metadata and final subtraction while leaving the IP in FP32, which returns -1472 for the case above, worse than the current 0.

The quadratic part of that expansion is grouped so the integer combination is formed before any conversion to double:

d1^2*S1 + d2^2*S2 - 2*d1*d2*Q  ==  d1*d2*(S1 + S2 - 2Q) + (d1 - d2)*(d1*S1 - d2*S2)

S1 + S2 - 2Q is sum((a[i] - b[i])^2), an exact non-negative integer. For two blobs sharing a delta the answer is therefore a non-negative integer scaled by delta^2 and cannot come out negative, and for a blob against itself both factors are exactly zero, so the distance is exactly zero regardless of how the compiler contracts the surrounding arithmetic. This is deliberately not fixed by clamping at zero: a clamp would also have turned the -201.5 the old metadata produced for a vector against itself into a clean 0, hiding the very defect this PR fixes.

The asymmetric L2 kernels derive ||x||^2 exactly, which removes the negative distances. Their remaining cancellation needs the difference taken inside the loop and is PR 2 in this stack.

UINT8_InnerProductImp returns the exact integer dot on all four ISAs. Callers that wanted a float now convert explicitly, including one that computed 1 - imp() in int arithmetic. On the plain uint8 side, ret_t is now 64-bit for every element type and stays signed, so 1 - imp() cannot underflow, and the two L2 reduces are read back unsigned, which costs no instructions on either ISA.

MAX_EXACT_DIM is removed rather than moved. It was not load-bearing, appearing only in its own definition and one debug assert, so release builds behaved identically with or without it, and its comment was wrong: 33,026 is where a signed 32-bit accumulator wraps, whereas Q_SUM_SQUARES stays exact in uint32 through dimension 66,051. For the record the real per-field ceilings are 66,051 for Q_SUM_SQUARES (L2 only), 16,843,009 for Q_SUM, and roughly 528,000 for the AVX512 per-lane accumulation. The first binds, and it is far above any dimension this library is used at.

Performance

This costs something, and the number should be visible rather than discovered later. Measured on an Ice Lake-SP host (Xeon Platinum 8375C, AVX512-VNNI), pinned core, 10 repetitions, every claim replicated across two independent runs per side and aggregated over dim >= 128, because that box produces reproducible 50% swings on small-dim scalar cases from binary layout alone.

Path	Effect
SQ8 and uint8 SIMD distance kernels	+1 to +4.8 ns per call, flat in ns rather than proportional, so about 7% at dim >= 128 after the mitigation below
SQ8_SQ8_NAIVE_* scalar kernels	23% to 49% faster, from integer accumulation
quantization path (QuantPreprocessor::quantize)	26% faster than before this PR

Two mitigations are included. sq8::reconstructed_* are marked always_inline, because GCC outlined reconstructed_l2_sqr and then built a realigned 64-byte stack frame plus a vzeroupper inside an otherwise leaf SIMD kernel just to call it; the AVX512 L2 kernel went from 15 to 54 instructions for that reason. Removing all 96 out-of-line calls and 33 of the stack realignments recovers 24% to 30% of the regression, with no case regressing against the unmitigated branch. Separately, std::round in the quantizer was an out-of-line call to libm round() per element, because that translation unit compiles at the x86-64 baseline and has no roundsd; since the value is already clamped to [0, 255] and therefore non-negative, adding 0.5 and truncating is std::round and compiles to a single cvttsd2si. That turns an 18% regression on the quantize path into a 26% improvement over main.

The residual 7% is the FP64 conversions and the dependency chain in the tail, which is the price of the exactness the rest of this PR establishes. Whether it is observable in real query latency is unknown: everything above is a hot-cache tight loop, whereas in an HNSW search each distance call follows a pointer chase that may stall on memory. There is no SQ8 index-level benchmark in the suite to answer that (MOD-14960), which is the second PR in this series to need one.

Also declined: narrowing the reduce back to 32 bits, which is safe at every reachable dimension and worth about 2% of a distance call, but costs an ARM-fork change and new ISA-specific test seams. Recorded for later rather than done here.

Verification

Original commits, built and run on an AVX512 + VNNI host (dorer-intel), debug:

Check	Result
library + unit tests compile	clean, 0 warnings
test_spaces	1461/1461, including 195 SQ8-to-SQ8 L2 and IP parameterizations
test_components	43/43
full ctest unit suite	2611/2611
check-format.sh	clean

The six later commits (uint8 accumulators, the integer regrouping, both performance changes, MAX_EXACT_DIM removal, the new tests) are clang-format clean and every touched translation unit passes -fsyntax-only, but were not built locally; CI is the first full build of them, and it is also the first execution of the NEON change on ARM hardware.

Test fixtures needed updating because they hand-build blobs. Two expectations moved, both because the quantity changed meaning rather than to make a test pass: the degenerate min == max case now expects both sums to be 0, since every value collapses to byte 0, and additionally asserts reconstructed_sum still recovers 3.5 * dim; and the FP16 preprocessor test's query-metadata assertions got their own baseline, having previously been compared against the storage sums, which only matched because both were taken over the input.

Regressions pinning each defect

Written against an FP64 reference computed over the reconstructed vectors, which is the quantity the kernels are trying to produce, rather than against numbers recorded from a run.

Test	Was
SQ8_SQ8_L2_survives_large_common_offset	0 instead of 64
SQ8_FP32_L2_is_non_negative_for_off_grid_query	-0.000490427, matching a radius-0 query
SQ8_SQ8_L2_matches_fp64_reference_at_high_dim	16% relative error at dim 16384
QuantizationHandlesNonRepresentableRange	undefined behaviour, NaN converted to uint8

Six further tests were added and run on both revisions to confirm each fails before this series and passes after it. They use only the quantizer test helper, sq8::MIN_VAL, sq8::DELTA and the public kernels, so the same source compiles on either side: a test written against Q_SUM or sq8::reconstructed_* would not compile before the change and so could not demonstrate anything.

Test	Was
SQ8_SQ8_L2_self_distance_is_exactly_zero	nonzero at every dimension, scalar and SIMD, sign varying with how individual elements rounded
SQ8_SQ8_L2_self_distance_whiteboard_case	exactly -201.5 for [0, 100.5, 255] against itself, with numbers small enough to check by hand
SQ8_SQ8_distance_ignores_sub_quantum_input_changes	two inputs quantizing to identical bytes answered the same query differently, by 2.1e-3 for L2 and 1.5e-3 for IP
SQ8_SQ8_IP_matches_reconstruction_dot	2.3e-3 relative error against the dot product of the reconstructions
UINT8_L2Sqr_and_InnerProduct_are_exact_past_int32	UB in the scalar path and negative L2 distances from the AVX512 and NEON reduces, at dimensions 33,026 and 40,000 with worst-case all-255 bytes

Dimensions are chosen so both kernel families run: 3, 5 and 8 dispatch to the naive implementations and 64 and above to AVX512F_BW_VL_VNNI on an Ice Lake host. The identical-bytes test keeps min nonzero on purpose, because with min == 0 the old IP formula collapses to delta_x * delta_y * q_dot, the input-derived sums drop out, and that half of the test would be inert.

The uint8 spaces benchmark fixture is also fixed: it paired new[] with delete and stored the trailing norms through unaligned float casts, so any measurement taken from it was suspect.

Also folded in: MOD-17528

The quantization UB lives in the same function and is the same question of whether the preprocessor's output can be trusted. [-FLT_MAX, +FLT_MAX] made max - min overflow to inf, then delta inf, inv_delta 0, and inf * 0 NaN, whose conversion to an integer is UB, reachable from AddVector for input the API accepts. The range and the per-element scaling now happen in double, which cannot overflow for any pair of finite floats, so the class disappears rather than being special-cased, and the result is clamped before conversion so an unexpected value would be a wrong byte rather than UB.

The alternative reformulation x * inv_delta - min * inv_delta was measured and rejected: it avoids the overflow but reintroduces cancellation for ordinary vectors with a large offset, quantizing them coarsely.

One behavioural note on the std::round replacement: the two forms differ only where x + 0.5 itself double-rounds, at scaled == 0.49999999999999994, where the new form yields byte 1 and std::round yields 0. There is no divergence at any of the 255 .5 boundaries, across 2M uniform random values, or across 3M realistic (x - min) * inv_delta draws. lrint and nearbyint are not substitutes: they round ties to even, which would move 178.5 to 178 and change bytes in common cases.

Not in this PR

• The asymmetric L2 cancellation, which needs the difference taken inside the loop. That is PR 2 in this stack, 13 new SIMD kernels.
• An SQ8 index-level benchmark to say whether the 7% per-call cost is observable end to end (MOD-14960).
• Narrowing the integer reduce back to 32 bits, worth about 2% of a distance call.
• ARM was not executed locally; only x86, and the later commits not even that. CI covers NEON and SVE.

Related: #1007 is blocked on this landing first, since it is what makes these paths reachable.

🤖 Generated with Claude Code

---

Note

Medium Risk
Touches core VecSim distance and quantization paths on every ISA with subtle numeric behavior; mitigated by broad unit regressions but asymmetric L2 cancellation and later commits may still need CI on ARM.

Overview
Makes SQ8 stored metadata describe the quantized reconstruction (Q_SUM / Q_SUM_SQUARES over bytes, same blob layout) and routes distance math through shared sq8::reconstructed_* helpers in double instead of duplicated inline formulas.

QuantPreprocessor now scales in double, accumulates integer sums over quantized bytes, and hardens edge cases (huge range, subnormal delta, centered ±inf) so quantization no longer hits NaN→uint8 UB.

Symmetric SQ8–SQ8 inner product and L2 use exact uint64_t byte dot products and reconstructed_l2_sqr (expanded difference) rather than ||x||² + ||y||² − 2·IP, fixing cancellation (e.g. large common offset → 0 distance) and negative self-distances. Asymmetric SQ8–FP32/FP16 L2 derives storage ||x||² via reconstructed_sum_squares and combines in double (full asymmetric cancellation deferred).

Plain uint8 paths widen accumulators to 64-bit (ret_t = long long, UINT8_InnerProductImp → uint64_t, unsigned L2 horizontal sums) so dimensions above ~33k do not wrap signed 32-bit totals.

Tests and quantizer fixtures updated; regressions added against FP64 references and prior failure modes.

^{Reviewed by Cursor Bugbot for commit 729f621. Bugbot is set up for automated code reviews on this repo. Configure here.}

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 97.17%. Comparing base (efd63da) to head (729f621).

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #1011      +/-   ##
==========================================
- Coverage   97.17%   97.17%   -0.01%     
==========================================
  Files         141      141              
  Lines        8328     8383      +55     
==========================================
+ Hits         8093     8146      +53     
- Misses        235      237       +2     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[cursor]

Cursor Bugbot has reviewed your changes using high effort and found 1 potential issue.

^{❌ Bugbot Autofix is OFF. To automatically fix reported issues with cloud agents, have a team admin enable autofix in the Cursor dashboard.}

^{Reviewed by Cursor Bugbot for commit 2b66b5f. Configure here.}

[lerman25]

Nice - 1 blocking comment

[lerman25]

Blocking: WithNorm can still make min_val/max_val non-finite before this double range calculation. Both find_min_max() and transformed_value() compute input[i] - mean[i] in FP32. For example, finite FP32 input [FLT_MAX, 0] with finite mean [-FLT_MAX, 0] centers to [+Inf, 0]; this then gives diff = Inf, delta = Inf, and inv_delta = 0, so to_byte(+Inf) evaluates Inf * 0 as NaN. std::clamp preserves NaN, and the following conversion to uint32_t is undefined behavior. This is reachable by the mean-centred SQ8 configuration introduced by #1007, so the finite-input safety claim is incomplete. Please perform/check centering in a representation that cannot overflow here (or reject non-finite derived values before quantization) and add a UBSan regression for this case.

[dor-forer]

Agreed, and fixed in 6c8d6e0. You are right that the double range does not help when the value is already lost upstream in FP32.

find_min_max now clamps its endpoints to the float range. That is also a storage constraint rather than only a guard: min_val is stored as FP32, so a non-finite endpoint could not be stored under any arithmetic. Only the endpoints are clamped, so the per-element loop stays FP32 and pays nothing.

With min_val finite and inv_delta finite and positive, NaN is unreachable in to_byte: a centered inf element gives inf * finite = inf, not 0 * inf, and lands on 255. I also moved the bound from std::clamp to std::fmin/std::fmax, since clamp propagates NaN, so the conversion is defined without relying on a proof that spans two functions.

I went with clamping rather than centering in double. The tradeoff: it puts overflowing elements at the ends instead of placing them proportionally. Happy to switch to double centering if you want proportionality, but it costs a conversion per element and the FP32 min slot still cannot hold the range.

UBSan regression added: QuantizationHandlesNonRepresentableCenteredRange, with your exact input.