[MOD-17526][MOD-17527] Make SQ8 metadata exact and fix the symmetric L2 formulation - #1011
[MOD-17526][MOD-17527] Make SQ8 metadata exact and fix the symmetric L2 formulation#1011dor-forer wants to merge 17 commits into
Conversation
WIP: source complete, test fixtures not yet updated, not yet compiled. Stored SQ8 metadata now describes the reconstruction rather than the input, and the sums are integers rather than FP32. * sq8.h: SUM and SUM_SQUARES become Q_SUM and Q_SUM_SQUARES, uint32 sums over the quantized bytes. Same slot count and blob size. Adds unaligned readers and four exact-derivation helpers (reconstructed_sum, reconstructed_sum_squares, reconstructed_ip, reconstructed_l2_sqr) so the layout has one interpretation. MAX_EXACT_DIM records the dimension past which sum(a^2) leaves uint32. * preprocessors.h: accumulates the quantized bytes in integers and packs mixed FP32/uint32 metadata. * The symmetric IP kernels previously combined sums taken over the original input with a dot product of the reconstructions, mixing two different vectors into one formula. All five now derive every term from the integer sums. * The symmetric L2 kernels used ||x||^2 + ||y||^2 - 2*IP in FP32, which cancels away the answer when the vectors share a large offset: stored [100000, 100008] against [100000, 100000] returned 0 where the truth is 64. They now use the expanded difference, which keeps the offset in a (min1 - min2) term formed once and exactly, and combine in double. Measured relative error drops to ~2e-8, including at dim 16384 where the FP32 combination was 16% off. * The asymmetric L2 kernels derive ||x||^2 exactly from the integer sums, which removes the negative distances (stored [0, 0.25, 1] against [0, 0.2501, 1] returned -0.000490427). Their remaining cancellation needs the difference taken inside the loop and is the next change in this stack. * UINT8_InnerProductImp now returns the exact integer dot on all four ISAs. 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 into an int, which overflows from dimension 33,027. Callers that wanted a float now convert explicitly, and one that computed 1 - imp() as int is now float arithmetic. This is the widening asked for separately, needed here because the L2 expansion depends on the dot being exact. Verified so far: clang-format clean; -Wall -Werror -fsyntax-only clean on both scalar TUs and on 17 of 19 x86 SQ8 kernels. The other two are not standalone compilable in an ad-hoc TU, on main as much as here, so they are covered only by a real build. Nothing has been executed yet. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The fixtures hand-build SQ8 blobs, so they encoded the old semantics. * tests_utils.h: the blob builder sums the quantized bytes as integers, and the FP16 reference L2 derives ||x||^2 through sq8::reconstructed_sum_squares_of and combines in double, matching the kernel it is a reference for. * test_spaces.cpp: the unaligned-metadata case stores uint32 sums. Its values are unchanged, since min 0 and delta 1 make the reconstruction equal to the bytes, so the expected distance does not move. It still asserts the metadata addresses are not float-aligned, which the uint32 loads must also tolerate. * test_components.cpp: adds a uint32 metadata reader beside the float one, and compares the FP16 path against the FP32 baseline on the integer sums. * The degenerate min == max case now expects both sums to be 0, because every value collapses to byte 0, and additionally asserts reconstructed_sum still recovers 3.5 * dim. That expectation moved because the quantity changed meaning, not because the number was adjusted to fit. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The FP16 preprocessor test compared query metadata against the storage sums. That only worked because both were taken over the input values; the storage sums now describe the reconstruction, so the query side needs its own expectation computed from the widened query values. Queries are never quantized, so their metadata stays FP32. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
ComputeSQ8Quantization builds the expected blob that the preprocessor tests byte-compare against, and it still wrote FP32 sums over the input. It did not reference the layout enum, so renaming the slots did not catch it; the byte comparison did. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Folds in the quantization UB, since it lives in the same function and is the same question of whether the preprocessor's output can be trusted. Finite FP32 input whose range is not representable, [-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 undefined behaviour. The range and the per-element scaling are now computed in double, which cannot overflow for any pair of finite floats, so the class disappears instead of being special-cased. The result is also clamped to [0, 255] before conversion, so a value the arithmetic should never produce would be a wrong byte rather than UB. Note the reformulation x * inv_delta - min * inv_delta was rejected: it avoids the overflow but reintroduces cancellation for ordinary vectors with a large offset, quantizing them coarsely. Regressions, all written against an FP64 reference over the reconstructed vectors rather than against numbers recorded from a run: * SQ8_SQ8_L2_survives_large_common_offset: stored [100000, 100008] against [100000, 100000], which returned 0 where the truth is 64. * SQ8_FP32_L2_is_non_negative_for_off_grid_query: stored [0, 0.25, 1] against [0, 0.2501, 1], which returned -0.000490427 and so matched a radius-0 query. * SQ8_SQ8_L2_matches_fp64_reference_at_high_dim: dim 16384 near-identical vectors, the regime where combining even the exact sums in FP32 was 16% off. * QuantizationHandlesNonRepresentableRange: the UB case above, under UBSan in CI. Also brought the QuantPreprocessor class documentation back in line with the code: it still described sums over the original values and presented the L2 identity as the symmetric formulation. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The large-offset regression asserted that the FP64 reference equals 64 within 1e-6 before comparing the kernel against it. The reference is 64.0000076: reconstructing through FP32 min and delta does not land exactly on 100008, so that residue is quantization error and belongs in the bound. The assertion that matters, kernel against reference, passed unchanged. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Codecov Report✅ All modified and coverable lines are covered by tests. 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. 🚀 New features to boost your workflow:
|
The uint8 kernels accumulate products of bytes, so the total reaches
255*255*dim = 65025*dim. Three paths could not hold that:
* IP.cpp / L2.cpp: ret_t for a 1-byte element type was int, so the
scalar UINT8_InnerProduct and UINT8_L2Sqr hit signed-overflow UB
from dimension 33,026, while the comment claimed support to 2^16.
ret_t is now 64-bit for every element type. Keeping it signed means
the "1 - ip" in the wrappers stays signed arithmetic and cannot
underflow, and the int8 paths are unaffected in behaviour.
* L2_AVX512F_BW_VL_VNNI_UINT8: the horizontal reduce was read back as
a signed int, so the distance went negative from dimension 33,026.
It is now read as uint32_t, which costs no instructions.
* L2_NEON_UINT8: same defect, int32_t receiving vaddvq_u32. Now
uint32_t. L2_NEON_DOTPROD_UINT8 and L2_SVE_UINT8 were already
unsigned and are unchanged.
This is the same defect class as MOD-17527, which this branch already
fixes on the inner product side; these are the L2 and scalar paths it
did not cover.
Also fixes the uint8 spaces benchmark fixture, which paired new[] with
delete and stored the trailing norms through unaligned float casts, so
that measurements taken from it are trustworthy.
The constant was not load-bearing: it appeared only in its own definition and in one debug assert, so release builds behaved identically with or without it, and it told the reader nothing that was true. Its comment claimed 33,026 was the uint32 exactness limit for sum(a[i]^2), but 65025 * dim stays exact in uint32 through dimension 66,051; 33,026 is where a *signed* 32-bit accumulator would have wrapped, which is a bound this code no longer has anywhere. For the record, the real per-field ceilings are: q_sum_squares (L2 only) exact through dim 66,051, q_sum through 16,843,009, and the AVX512 per-lane accumulation through roughly 528,000. The first is the binding one, and it is far above any dimension the library is used at.
std::round compiles to an out-of-line call to round() here, because this translation unit is built at the x86-64 baseline and so has no roundsd available; nm confirms the symbol and the object contains no rounding instruction. That call ran once per element. scaled is already clamped to [0, 255] and therefore non-negative, so adding 0.5 and truncating is std::round, and compiles to a single cvttsd2si. Measured on Ice Lake-SP, this makes the quantization path about 26% faster than before this branch, which also absorbs the ~18% this branch had added to it. The two forms differ only where x + 0.5 itself double-rounds, at scaled == 0.49999999999999994, where this yields byte 1 and std::round yields 0. No divergence at any of the 255 .5 boundaries, across 2M uniform random values, or across 3M realistic (x - min) * inv_delta draws. Note that lrint and nearbyint are NOT substitutes: they round ties to even, which would move 178.5 to 178 and change bytes in common cases.
reconstructed_l2_sqr returned a slightly negative squared distance for a
blob against itself: about -6e-14 at dimension 64 and -1.7e-12 at 512 on
the AVX512 kernel. The six terms cancel exactly in real arithmetic, but
the compiler contracts some of the products into fused multiply-adds and
not others, so they no longer round identically and no longer cancel.
Harmless to ranking, but a caller taking a square root of it gets NaN,
and a function named L2Sqr should not return a negative number.
Regrouping the quadratic part fixes it structurally:
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, so
for equal deltas the answer is a non-negative integer scaled by delta^2
and cannot be negative, and for identical blobs both factors are exactly
zero. Multiplying by an exact zero is exact, so contraction cannot break
it. Measured against a long double reference with -mfma: self-distance
exactly zero at dimensions 3 through 4096, mean relative error on random
pairs 5.6e-17 versus 2.9e-16 before, and in the regime this branch exists
to fix, near-duplicates sharing a large offset, exactly zero error
versus 1e-9 before.
Deliberately not fixed by clamping at zero. A clamp would have turned the
-201.5 the old metadata produced into a clean 0 and hidden the very defect
this branch fixes.
It also shortens the dependency chain, since the integer combination is
independent of the delta terms.
GCC outlines reconstructed_l2_sqr, plausibly because the 64 residual instantiations of each SQ8 kernel all call it. The cost is not just the call: in the AVX512 kernel it also sets up a realigned 64-byte stack frame and a vzeroupper inside what is otherwise a leaf SIMD function. SQ8_SQ8_L2SqrSIMD64_AVX512F_BW_VL_VNNI<0> went from 15 instructions to 54 for this reason. Marking the six helpers always_inline removes all 96 out-of-line calls in that translation unit and 33 of the 48 stack realignments, and recovers 24% to 30% of the per-call regression this branch introduced (median +2.15 ns to +1.44 ns on spaces_sq8_sq8, replicated across two runs per side on Ice Lake-SP). No case regressed against the unpatched branch. The cost is about 11 extra instructions per instantiation, since a 24-instruction body is now duplicated rather than shared.
Six tests, each of which fails on the commit before this series and
passes after it. They were run on both revisions to confirm that.
They use only the quantizer test helper, sq8::MIN_VAL, sq8::DELTA and
the public kernels, so the same source compiles on either side. That is
deliberate: a test written against Q_SUM or sq8::reconstructed_* would
not compile before the change and so could not demonstrate anything.
* self-distance is exactly zero, scalar and dispatched, five
dimensions. Was nonzero at every dimension, sign varying with the
rounding of individual elements.
* the same property with hand-checkable numbers, where the old result
is exactly -201.5 for a vector against itself.
* two inputs that quantize to identical bytes produce identical
distances. The distance must be a function of the stored blob and
nothing else; it differed by 2.1e-3 for L2 and 1.5e-3 for IP.
* the inner product equals the dot product of the reconstructions,
which was off by 2.3e-3 relative.
* plain uint8 L2 and inner product stay exact past INT_MAX, at
dimensions 33,026 and 40,000 with worst-case all-255 bytes. The
scalar path was UB there and the AVX512 and NEON L2 reduces returned
negative distances.
The dimensions are chosen so that both the scalar kernels and the SIMD
kernels are exercised: 3, 5 and 8 dispatch to the naive implementation
and 64 and above to AVX512F_BW_VL_VNNI on an Ice Lake host.
The third test keeps min nonzero on purpose. With min == 0 the old inner
product formula collapses to delta_x * delta_y * q_dot, the input-derived
sums drop out, and the inner product half of the test cannot detect
anything.
make check-format was failing on it.
A list of the source paths touched by this change set, committed by accident in 2b84cc2. Nothing in the repo, build or tests references it. Reported by Cursor Bugbot.
There was a problem hiding this comment.
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.
The test freed both blobs but never deleted the QuantPreprocessor it allocated, leaking 88 bytes in 3 allocations. LeakSanitizer failed it in the asan and coverage jobs; a plain debug ctest run does not catch this, which is why it passed locally. Every sibling test in the file already deletes its preprocessor.
lerman25
left a comment
There was a problem hiding this comment.
Nice - 1 blocking comment
| // FP32 arithmetic could leave the representable range for input that is entirely valid. | ||
| // [-FLT_MAX, +FLT_MAX] made max - min overflow to inf, then delta inf, inv_delta 0, and | ||
| // finally inf * 0 = NaN, whose conversion to an integer is undefined behaviour. Doubles | ||
| // cannot overflow for any pair of finite floats, so the whole class disappears rather than |
There was a problem hiding this comment.
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.
There was a problem hiding this comment.
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.
Widening the range to double closed the path MOD-17528 described and left two others open, both reported in review and both reaching the same conversion of NaN to an integer. * find_min_max centers in FP32, so input FLT_MAX against mean -FLT_MAX centers to 6.8e38, i.e. inf, before the double range is ever computed. min_val is stored as FP32, so a non-finite endpoint could not be stored under any arithmetic; the endpoints are now clamped to the float range, which is both the storage limit and the guard. Only the endpoints are clamped, so the per-element loop stays FP32 and pays nothing. * delta is stored as FP32, and (float)(diff / 255.0) underflows to zero for any diff below about 1.8e-43. diff itself is not zero there, so testing diff did not catch it. 1/delta was then inf and the minimum element, whose numerator is exactly zero, scaled to 0 * inf = NaN. The guard now tests delta, which is the value that gets stored and inverted, and subsumes the old min == max check. With min_val finite and inv_delta finite and positive, NaN is unreachable in to_byte. The bound there is still moved from std::clamp to std::fmin/std::fmax, which return the non-NaN operand where clamp propagates NaN, so the conversion is defined on its own terms rather than on a proof spanning two functions. This function has now produced three separate conversion-UB defects, which is what makes the backstop worth one instruction. Both cases are pinned by UBSan regressions next to the existing one. Also drops the claim that widening the subtraction removed the whole class, which is what the review was disputing, and names the two places that finish the job instead. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
The reasoning had accumulated as a block above nearly every statement, which buried the code it was explaining. It is now one description above quantize(), covering why the scaling is in double and which three guards keep NaN out of the byte conversion. Two inline comments are kept: the +0.5 trick, which reads like an oversight and invites being changed back to std::round, and one naming what transformed_value does. Also drops the pre-existing "We know (input - min) => 0", which stopped being true once the endpoints were clamped and would have someone read the fmax bound in to_byte as dead code. Comments only. With comments and blank lines stripped, the file is identical to the previous commit. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
The uint8 kernels accumulate products or squared differences of bytes, so
the total reaches 255 * 255 * dim = 65025 * dim. Four paths could not hold
that, and all of them are on the plain int8/uint8 index paths that ship
today.
* IP.cpp / L2.cpp: ret_t for a 1-byte element type was int, so the
scalar UINT8_InnerProduct and UINT8_L2Sqr executed signed-overflow UB
from dimension 33,026, while the comment claimed support to 2^16.
ret_t is now 64-bit for every element type. Keeping it signed means the
"1 - ip" in the wrappers stays signed arithmetic and cannot underflow,
and the int8 paths are unaffected. The L2 comment still carried the old
"at least 2 bytes wider" rationale and is corrected.
* UINT8_InnerProductImp returned float on NEON and SVE, which
accumulated exactly in integer lanes and then discarded it, exact only
to dimension 258 since 2^24 / 65025 = 258.
* AVX512 reduced 16 int32 lanes with _mm512_reduce_add_epi32 into a
signed int, wrapping from dimension 33,026.
* L2_AVX512F_BW_VL_VNNI_UINT8 and L2_NEON_UINT8 read an unsigned
horizontal reduce back into a signed int, so the distance went negative
from the same dimension. L2_NEON_DOTPROD_UINT8 and L2_SVE_UINT8 were
already unsigned and are unchanged.
Note the accumulation itself was never the problem. The SIMD adds wrap
modulo 2^32 and are bit-exact, so the bit pattern was already correct; the
top bit was being read as a sign. An unsigned 32-bit reduce therefore costs
nothing over the original and is exact through dimension 66,051, twice the
old signed limit of 33,025.
Above 66,051 a 32-bit total genuinely does run out, so each kernel gains a
`bool Wide` template parameter selecting the epilogue: the narrow unsigned
32-bit reduce, or a widening one that zero-extends the lanes to 64 bits
first and cannot wrap at any dimension. The lanes are accumulated
identically either way. The choosers pick once per index, so no branch
enters the kernel.
Widening unconditionally would have been simpler and was measured rather
than assumed. On an Ice Lake-SP Xeon it costs 4 extra uops in the epilogue:
dim 32 +20%
dim 55-200 +8 to +11%
dim 256 +7%
dim 900-1024 +4 to +5%
15 repetitions, pinned core, two passes with the A/B order reversed; sign
and magnitude hold across both. The loop bodies are instruction-for-
instruction identical with the loop tops at the same 32-byte offset, so
this is the epilogue alone. Instruction count understates it, because the
widening reduce lengthens a dependency chain rather than adding throughput
work; that is also why the absolute delta grows at high dim, where fewer
calls overlap to hide the latency.
Selecting per dimension keeps that cost off every ordinary index. The price
is instantiating both variants: on the AVX512F_BW_VL_VNNI translation unit
at -O2, object size goes from 514,792 to 657,592 bytes, +27.7%, with 197
extra exported symbols. The narrow instantiation is unchanged at 40
instructions for residual 32, before and after, so the common case keeps
the full benefit.
To avoid a second case ladder, CHOOSE_IMPLEMENTATION now forwards trailing
arguments as further template arguments using __VA_OPT__, so the same
ladder serves kernels templated on <residual> and on <residual, Wide>
alike and every existing call site is untouched.
CHOOSE_UINT8_IMPLEMENTATION wraps the dimension test so each of the 15
uint8 call sites is a one-word change, and
CHOOSE_SVE_UINT8_IMPLEMENTATION does the same for the SVE ladder.
The SQ8-to-SQ8 inner product kernels reuse this helper, on main as much as
here, so they now pass Wide explicitly. They pass false: SQ8 is capped at
the same dimension independently, because its q_sum_squares metadata slot
is a uint32 holding 65025 * dim, so a widening reduce there would exceed
what the metadata itself can represent.
Split out of #1011 because none of this depends on the SQ8 metadata
contract that PR is changing, while all of it affects code reachable today.
#1011 does depend on this, through the helper above.
The regressions use all-255 bytes, the worst case, which makes the expected
value an exact integer: dimensions 33,026 and 40,000 for the narrow path,
which the old signed reduce got wrong, and 66,052 and 80,000 for the wide
path. The existing UINT8 suites stop at dimension 128, which is why all of
this went unseen; being SIMD-versus-scalar comparisons they would also have
agreed with each other wherever both wrapped.
Also fixes the uint8 spaces benchmark fixture, which paired new[] with
delete and stored the trailing norms through unaligned float casts, so
measurements taken from it can be trusted.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
QuantPreprocessor::quantize could execute undefined behaviour for input
the API accepts, in three ways that all end at the same cast.
The scale was derived entirely in FP32. For [-FLT_MAX, +FLT_MAX], every
component finite and accepted without validation, max - min overflowed to
inf, delta became inf, inv_delta 0, and the per-element product
inf * 0 = NaN, whose conversion to an integer is undefined. UBSan:
"-nan is outside the range of representable values of type 'unsigned
char'". The existing diff == 0 guard covers equal values, not overflow of
the subtraction. (MOD-17528)
Two more paths reach the same cast and were found in review of the
follow-up work:
* With WithNorm, centering is an FP32 subtraction, so finite input
against a finite mean can produce inf before any range is computed:
FLT_MAX against mean -FLT_MAX centers to 6.8e38. Widening the range
to double does not help, because the value is already lost upstream.
* delta is stored as FP32, and (float)(diff / 255) underflows to zero
for any diff below about 1.8e-43 while diff itself is nonzero, so
testing diff does not catch it. 1/delta was then inf, and the minimum
element, whose numerator is exactly zero, scaled to 0 * inf.
There is also no rejection path: nothing in VecSim validates finiteness,
and AddVector has no way to report "unquantizable", so the contract has to
be saturation rather than an error.
All three are closed by normalizing the endpoints once, immediately after
find_min_max, covering both the plain and WithNorm branches. Both
endpoints get a two-sided clamp behind an order check that catches NaN. A
one-sided clamp is not enough: for an all-+inf vector inf <= inf passes
the order check, so std::max(+inf, -FLT_MAX) would leave min at +inf and
store it. min is stored as FP32, so a non-finite endpoint could not be
represented under any arithmetic; this is the storage limit as much as a
guard.
That bounds diff at 6.8e38, so delta can be neither inf nor NaN and only
the underflow guard remains, as one comparison. inv_delta stays double,
not for precision, which needs only +/-0.5 in 255, but because an FP32
reciprocal overflows for a subnormal delta: [0, 7e-37] gives delta
2.7e-39, whose FP64 reciprocal is finite and correctly maps the top
element to 255.
The per-element bound is written by hand rather than with std::clamp or
std::fmin/std::fmax. Each of those breaks something: std::clamp is
comparisons and propagates NaN into the cast, while fmin/fmax are
NaN-correct but compile to two out-of-line libm calls per element at this
translation unit's baseline. Measured at -O3: 8 instructions for this
form against 9 for std::clamp and 10 plus two calls for fmin/fmax. It
also subsumes the std::round that was there, since bounding first makes
+0.5 and truncation equivalent, and round() is likewise out-of-line here.
Tests: the MOD-17528 reproduction, the WithNorm centering case, and a
table-driven domain matrix covering constant vectors positive negative and
zero, a subnormal but representable delta, a range that underflows and
collapses, the full FP32 range, all +inf, all -inf, mixed infinities, and
NaN first, middle and last. Expected bytes and metadata are asserted
rather than a range check, which is vacuous for uint8_t. The three NaN
cases pin that position matters: std::minmax_element compares with < and
every comparison against NaN is false, so a NaN at either end reaches an
endpoint and trips the order check while one in the middle is skipped and
the finite values set a real range. Expectations were derived by
simulating the pipeline, which corrected three of them.
Metadata meaning is unchanged: the sums are still FP32 over the input
values. Making them exact integer sums over the quantized bytes is a
storage-contract change that has to move together with every kernel that
reads them, and stays in #1011.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
The uint8 kernels accumulate products or squared differences of bytes, so
the total reaches 255 * 255 * dim = 65025 * dim. Four paths could not hold
that, and all of them are on the plain int8/uint8 index paths that ship
today.
* IP.cpp / L2.cpp: ret_t for a 1-byte element type was int, so the
scalar UINT8_InnerProduct and UINT8_L2Sqr executed signed-overflow UB
from dimension 33,026, while the comment claimed support to 2^16.
ret_t is now 64-bit for every element type. Keeping it signed means the
"1 - ip" in the wrappers stays signed arithmetic and cannot underflow,
and the int8 paths are unaffected. The L2 comment still carried the old
"at least 2 bytes wider" rationale and is corrected.
* UINT8_InnerProductImp returned float on NEON and SVE, which
accumulated exactly in integer lanes and then discarded it, exact only
to dimension 258 since 2^24 / 65025 = 258.
* AVX512 reduced 16 int32 lanes with _mm512_reduce_add_epi32 into a
signed int, wrapping from dimension 33,026.
* L2_AVX512F_BW_VL_VNNI_UINT8 and L2_NEON_UINT8 read an unsigned
horizontal reduce back into a signed int, so the distance went negative
from the same dimension. L2_NEON_DOTPROD_UINT8 and L2_SVE_UINT8 were
already unsigned and are unchanged.
The accumulation itself was never the problem. The SIMD adds wrap modulo
2^32 and are bit-exact, so the bit pattern was already correct; the top bit
was being read as a sign. An unsigned 32-bit reduce therefore costs nothing
over the original and is exact through dimension 66,051, twice the old
signed limit of 33,025. Verified: the AVX512 object is byte-for-byte the
same size as before at 514,792, with the same 33 and 40 instructions for
residual 0 and 32.
Above 66,051 the choosers hand back the scalar kernel, which after the
ret_t change is exact to roughly dimension 2.8e14. That is one comparison
at index creation, reusing the "if (dim < 32) return ret_dist_func" idiom
the choosers already had, and it leaves every kernel untouched.
Two alternatives were explored and rejected, both recorded on the constant:
* Widening the horizontal reduce. Measured on an Ice Lake-SP Xeon it
costs 4 extra uops in the epilogue: +20% at dimension 32, +8-11% across
55-200, +4-5% at 900-1024, on byte-identical loop code. It also only
moves the limit, and to a different place per ISA, since NEON combines
four accumulators with vaddq_u32 in 32 bits before any widening reduce
sees them, capping it at 264,204 rather than the 1,056,816 AVX512 gets.
* Chunking the accumulation and flushing into a 64-bit total. Exact at
any dimension, and cheap when the chunk loop lives in the wrapper
rather than the kernel: +2 instructions on the fast path against +12 to
+21 when placed inside. Deferred rather than dismissed, since it is
only worth the restructuring if such dimensions become real.
Nothing comparable supports that range today, which is what settles it.
Lucene caps its scalar-quantized format at 1,024 dimensions and
Elasticsearch caps dense vectors at 4,096, both keeping a 32-bit
accumulator safe by contract. Faiss's QT_8bit_direct accumulates
full-range bytes into 32-bit lanes with no widening and carries the same
theoretical limit. Qdrant quantizes to 0..127, lowering the per-element cap
to 16,129, and its raw uint8 metric still sums into i32. The scalar
fallback here is already stricter than any of them.
Split out of #1011 because none of this depends on the SQ8 metadata
contract that PR is changing, while all of it affects code reachable today.
#1011 does depend on this, since its SQ8_SQ8 kernels call
UINT8_InnerProductImp. Note SQ8 is independently capped at the same
dimension, because q_sum_squares is a uint32 slot holding 65025 * dim, so
that PR needs its own fence regardless.
The regressions use all-255 bytes, the worst case, which makes the expected
value an exact integer, at dimensions 33,026 and 40,000 for the SIMD path
and 66,052 for the fallback. The fallback test asserts the returned
function pointer, not just the distance: on a host with no uint8 SIMD tier
the value comparison would pass either way, but the pointer identity would
not. The existing UINT8 suites stop at dimension 128, which is why all of
this went unseen; being SIMD-versus-scalar comparisons they would also have
agreed with each other wherever both wrapped.
Also fixes the uint8 spaces benchmark fixture, which paired new[] with
delete and stored the trailing norms through unaligned float casts, so
measurements taken from it can be trusted.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
QuantPreprocessor::quantize could execute undefined behaviour for input
the API accepts, in three ways that all end at the same cast.
The scale was derived entirely in FP32. For [-FLT_MAX, +FLT_MAX], every
component finite and accepted without validation, max - min overflowed to
inf, delta became inf, inv_delta 0, and the per-element product
inf * 0 = NaN, whose conversion to an integer is undefined. UBSan:
"-nan is outside the range of representable values of type 'unsigned
char'". The existing diff == 0 guard covers equal values, not overflow of
the subtraction. (MOD-17528)
Two more paths reach the same cast and were found in review of the
follow-up work:
* With WithNorm, centering is an FP32 subtraction, so finite input
against a finite mean can produce inf before any range is computed:
FLT_MAX against mean -FLT_MAX centers to 6.8e38. Widening the range
to double does not help, because the value is already lost upstream.
* delta is stored as FP32, and (float)(diff / 255) underflows to zero
for any diff below about 1.8e-43 while diff itself is nonzero, so
testing diff does not catch it. 1/delta was then inf, and the minimum
element, whose numerator is exactly zero, scaled to 0 * inf.
There is also no rejection path: nothing in VecSim validates finiteness,
and AddVector has no way to report "unquantizable", so the contract has to
be saturation rather than an error.
All three are closed by normalizing the endpoints once, immediately after
find_min_max, covering both the plain and WithNorm branches. Both
endpoints get a two-sided clamp behind an order check that catches NaN. A
one-sided clamp is not enough: for an all-+inf vector inf <= inf passes
the order check, so std::max(+inf, -FLT_MAX) would leave min at +inf and
store it. min is stored as FP32, so a non-finite endpoint could not be
represented under any arithmetic; this is the storage limit as much as a
guard.
That bounds diff at 6.8e38, so delta can be neither inf nor NaN and only
the underflow guard remains, as one comparison. inv_delta stays double,
not for precision, which needs only +/-0.5 in 255, but because an FP32
reciprocal overflows for a subnormal delta: [0, 7e-37] gives delta
2.7e-39, whose FP64 reciprocal is finite and correctly maps the top
element to 255.
The per-element bound is written by hand rather than with std::clamp or
std::fmin/std::fmax. Each of those breaks something: std::clamp is
comparisons and propagates NaN into the cast, while fmin/fmax are
NaN-correct but compile to two out-of-line libm calls per element at this
translation unit's baseline. Measured at -O3: 8 instructions for this
form against 9 for std::clamp and 10 plus two calls for fmin/fmax. It
also subsumes the std::round that was there, since bounding first makes
+0.5 and truncation equivalent, and round() is likewise out-of-line here.
Tests: the MOD-17528 reproduction, the WithNorm centering case, and a
table-driven domain matrix covering constant vectors positive negative and
zero, a subnormal but representable delta, a range that underflows and
collapses, the full FP32 range, all +inf, all -inf, mixed infinities, and
NaN first, middle and last. Expected bytes and metadata are asserted
rather than a range check, which is vacuous for uint8_t. The three NaN
cases pin that position matters: std::minmax_element compares with < and
every comparison against NaN is false, so a NaN at either end reaches an
endpoint and trips the order check while one in the middle is skipped and
the finite values set a real range. Expectations were derived by
simulating the pipeline, which corrected three of them.
Metadata meaning is unchanged: the sums are still FP32 over the input
values. Making them exact integer sums over the quantized bytes is a
storage-contract change that has to move together with every kernel that
reads them, and stays in #1011.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
…nput QuantPreprocessor::quantize could execute undefined behaviour for input the API accepts. Three paths, all ending at the same conversion to a byte, and all reachable with entirely finite components. The scale was derived in FP32, so for [-FLT_MAX, +FLT_MAX] the range overflowed: max - min became inf, delta inf, inv_delta 0, and the per-element product inf * 0 = NaN, whose conversion to an integer is undefined. UBSan: "-nan is outside the range of representable values of type 'unsigned char'". The existing diff == 0 guard covers equal values, not overflow of the subtraction. (MOD-17528) With WithNorm, centering is an FP32 subtraction, so finite input against a finite mean reaches inf before any range is computed: FLT_MAX against mean -FLT_MAX centers to 6.8e38. Widening the range does not help, because the value is already lost upstream. And delta is stored as FP32, so (float)(diff / 255) underflows to zero while diff itself is nonzero, which leaves 1/delta inf and scales the minimum element, numerator exactly zero, to 0 * inf. All three are closed by normalizing the endpoints once, after find_min_max, covering the plain and WithNorm branches together, after which a single delta comparison suffices. Both endpoints get a two-sided clamp: for an all-+inf vector inf <= inf passes the order check, so a one-sided std::max would leave min at +inf and store it. min is stored as FP32, so a non-finite endpoint could not be represented under any arithmetic. inv_delta stays double, not for precision, which needs only +/-0.5 in 255, but because an FP32 reciprocal overflows for a subnormal delta: [0, 7e-37] gives delta 2.7e-39, whose FP64 reciprocal is finite and correctly maps the top element to 255. The per-element bound and rounding are written by hand. That replaces std::round, an out-of-line libm call at this translation unit's baseline that ran once per element: bounding first makes adding 0.5 and truncating equivalent for non-negative values. Measured at -O3 it is 8 instructions against 9 for std::clamp plus std::round's call. Deliberately NOT claimed: that the function is defined for non-finite components. It is not, and cannot be made so here. find_min_max uses std::minmax_element, whose precondition is that the comparison induce a strict weak ordering, and floating-point < is not one once a NaN is present: incomparability must be transitive, yet 1.0 is incomparable with NaN and NaN with 2.0 while 1.0 and 2.0 are comparable. The undefined behaviour is therefore inside that algorithm, before quantize() ever sees a range, so no normalization afterwards can define a portable result. A partial contract would also be misleading, since x_mean_ip, the quantized sums and the whole query metadata path are untouched and can still produce non-finite values. Non-finite components are treated as unsupported. The order check remains as a defensive fallback so that a NaN endpoint cannot be stored, which degrades a caller error into meaningless-but-finite metadata rather than poisoning every distance computed against that vector. Validating at the public ingestion boundary belongs in its own change; nothing in VecSim does it today. This matches how comparable systems handle it. Lucene validates vector components and throws on NaN or infinity, and Elasticsearch rejects NaN, infinity and magnitudes that overflow before delegating to Lucene. Faiss guards only an exactly zero range and assumes finite input otherwise, with the same float-to-integer concern at its final cast. Qdrant gets defined bytes from Rust's saturating cast, which C++ does not have, and can still store an unusable scale. Tests: the MOD-17528 reproduction, the WithNorm centering case, and a table-driven domain matrix over constant vectors positive negative and zero, a single element, a subnormal but representable delta, a range that underflows and collapses, the full FP32 range, and all-+inf, all--inf and mixed infinities. Infinities are pinned exactly, since < remains a strict weak ordering over finites and +/-inf; only NaN breaks it. Expected bytes and metadata are asserted rather than a range check, which is vacuous for uint8_t, and every expectation was derived by simulating the pipeline rather than predicted. NaN input gets one test that asserts only that the stored min and delta stay finite and delta positive, which is position-independent and portable, and which under UBSan also covers the conversion. Metadata meaning is unchanged: the sums are still FP32 over the input values. Making them exact integer sums over the quantized bytes is a storage-contract change that has to move together with every kernel that reads them, and stays in #1011. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
…nput QuantPreprocessor::quantize could execute undefined behaviour for input the API accepts. Three paths, all ending at the same conversion to a byte, and all reachable with entirely finite components. The scale was derived in FP32, so for [-FLT_MAX, +FLT_MAX] the range overflowed: max - min became inf, delta inf, inv_delta 0, and the per-element product inf * 0 = NaN, whose conversion to an integer is undefined. UBSan: "-nan is outside the range of representable values of type 'unsigned char'". The existing diff == 0 guard covers equal values, not overflow of the subtraction. (MOD-17528) With WithNorm, centering is an FP32 subtraction, so finite input against a finite mean reaches inf before any range is computed: FLT_MAX against mean -FLT_MAX centers to 6.8e38. Widening the range does not help, because the value is already lost upstream. And delta is stored as FP32, so (float)(diff / 255) underflows to zero while diff itself is nonzero, which leaves 1/delta inf and scales the minimum element, numerator exactly zero, to 0 * inf. find_min_max now guarantees that both endpoints are finite and ordered, and everything else follows from that. The guarantee lives there rather than at the call site because that is where it can be broken: the WithNorm branch creates the inf itself, from two valid operands, and the plain branch passes through whatever the input holds. Both endpoints get a two-sided clamp, since for an all-+inf vector inf <= inf passes the order check and a one-sided std::max would leave min at +inf. Clamping to the float range is not merely defensive either: min is stored as an FP32 field, so a non-finite endpoint could not be represented under any arithmetic. With that established, quantize needs one delta comparison, and inv_delta stays double. Not for precision, which needs only +/-0.5 in 255, but because an FP32 reciprocal overflows for a subnormal delta: [0, 7e-37] gives delta 2.7e-39, whose FP64 reciprocal is finite and correctly maps the top element to 255. The per-element bound and rounding are written by hand. That replaces std::round, an out-of-line libm call at this translation unit's baseline that ran once per element: bounding first makes adding 0.5 and truncating equivalent for non-negative values. Measured at -O3 it is 8 instructions against 9 for std::clamp plus std::round's call. Deliberately NOT claimed: that the function is defined for non-finite components. It is not, and cannot be made so here. find_min_max uses std::minmax_element, whose precondition is that the comparison induce a strict weak ordering, and floating-point < is not one once a NaN is present: incomparability must be transitive, yet 1.0 is incomparable with NaN and NaN with 2.0 while 1.0 and 2.0 are comparable. The undefined behaviour is therefore inside that algorithm, before any range exists, so no normalization afterwards can define a portable result. A partial contract would also be misleading, since x_mean_ip, the quantized sums and the whole query metadata path are untouched and can still produce non-finite values. Non-finite components are treated as unsupported. The order check remains as a defensive fallback so that a NaN endpoint cannot be stored, which degrades a caller error into meaningless-but-finite metadata rather than poisoning every distance computed against that vector. Validating at the public ingestion boundary belongs in its own change; nothing in VecSim does it today. This matches how comparable systems handle it. Lucene validates vector components and throws on NaN or infinity, and Elasticsearch rejects NaN, infinity and magnitudes that overflow before delegating to Lucene. Faiss guards only an exactly zero range and assumes finite input otherwise, with the same float-to-integer concern at its final cast. Qdrant gets defined bytes from Rust's saturating cast, which C++ does not have, and can still store an unusable scale. Tests: the MOD-17528 reproduction, the WithNorm centering case, and a table-driven domain matrix over constant vectors positive negative and zero, a single element, a subnormal but representable delta, a range that underflows and collapses, the full FP32 range, and all-+inf, all--inf and mixed infinities. Infinities are pinned exactly, since < remains a strict weak ordering over finites and +/-inf; only NaN breaks it. Expected bytes and metadata are asserted rather than a range check, which is vacuous for uint8_t, and every expectation was derived by simulating the pipeline rather than predicted. NaN input gets one test that asserts only that the stored min and delta stay finite and delta positive, which is position-independent and portable, and which under UBSan also covers the conversion. Metadata meaning is unchanged: the sums are still FP32 over the input values. Making them exact integer sums over the quantized bytes is a storage-contract change that has to move together with every kernel that reads them, and stays in #1011. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
…nput
QuantPreprocessor::quantize could execute undefined behaviour for input the
API accepts. Three paths, all ending at the same conversion to a byte, and
all reachable with entirely finite components.
The scale was derived in FP32, so for [-FLT_MAX, +FLT_MAX] the range
overflowed: max - min became inf, delta inf, inv_delta 0, and the
per-element product inf * 0 = NaN, whose conversion to an integer is
undefined. UBSan: "-nan is outside the range of representable values of
type 'unsigned char'". The existing diff == 0 guard covers equal values,
not overflow of the subtraction. (MOD-17528)
With WithNorm, centering is an FP32 subtraction, so finite input against a
finite mean reaches inf before any range is computed: FLT_MAX against mean
-FLT_MAX centers to 6.8e38. Widening the range does not help, because the
value is already lost upstream.
And delta is stored as FP32, so (float)(diff / 255) underflows to zero
while diff itself is nonzero, which leaves 1/delta inf and scales the
minimum element, numerator exactly zero, to 0 * inf.
find_min_max now guarantees that both endpoints are finite and ordered, and
everything else follows from that. The guarantee lives there rather than at
the call site because that is where it can be broken: the WithNorm branch
creates the inf itself, from two valid operands, and the plain branch passes
through whatever the input holds. Both endpoints get a two-sided clamp,
since for an all-+inf vector inf <= inf passes the order check and a
one-sided std::max would leave min at +inf. Clamping to the float range is
not merely defensive either: min is stored as an FP32 field, so a non-finite
endpoint could not be represented under any arithmetic.
With that established, quantize needs one delta comparison, and inv_delta
stays double. Not for precision, which needs only +/-0.5 in 255, but because
an FP32 reciprocal overflows for a subnormal delta: [0, 7e-37] gives delta
2.7e-39, whose FP64 reciprocal is finite and correctly maps the top element
to 255.
The per-element bound and rounding are written by hand. That replaces
std::round, an out-of-line libm call at this translation unit's baseline
that ran once per element: bounding first makes adding 0.5 and truncating
equivalent for non-negative values. Measured at -O3 it is 8 instructions
against 9 for std::clamp plus std::round's call.
Scope, stated precisely, because it is narrower than it might look:
* The byte conversion is defined for all finite components. That is the
goal and it is met.
* The stored min and delta are finite with delta positive. The *sums* are
not covered: they are accumulated in FP32 over the input values, so
[-FLT_MAX, +FLT_MAX] stores sum_squares as inf even though this
function's own arithmetic is now well defined. Separate problem,
separate change.
* Non-finite components are unsupported, and nothing here is tested
against them. std::minmax_element requires its comparison to induce a
strict weak ordering, and floating-point < is not one once a NaN is
present: incomparability must be transitive, yet 1.0 is incomparable
with NaN and NaN with 2.0 while 1.0 and 2.0 are comparable. Violating
that precondition is undefined behaviour inside the algorithm, before
any range exists, so no assertion about the outcome would be portable,
including a weak one about metadata finiteness. UBSan reporting nothing
does not establish otherwise. Rejecting non-finite components at the
public ingestion boundary is the actual fix and belongs in its own
change; nothing in VecSim does it today.
The order check on the endpoints stays, as an invariant guard rather than a
NaN policy: everything downstream is written assuming min <= max, and
asserting that once per vector costs less than reasoning about whether it
can be violated.
This matches how comparable systems handle it. Lucene validates vector
components and throws on NaN or infinity, and Elasticsearch rejects NaN,
infinity and magnitudes that overflow before delegating to Lucene. Faiss
guards only an exactly zero range and assumes finite input otherwise, with
the same float-to-integer concern at its final cast. Qdrant gets defined
bytes from Rust's saturating cast, which C++ does not have, and can still
store an unusable scale.
Tests: the MOD-17528 reproduction, the WithNorm centering case, and a
table-driven matrix over finite input only: constant vectors positive
negative and zero, a single element, a subnormal but representable delta, a
range that underflows and collapses, and the full FP32 range. Expected bytes
and scale metadata are asserted rather than a range check, which is vacuous
for uint8_t, and every expectation was derived by simulating the pipeline
rather than predicted. The parameterized test is named for what it checks,
ScaleMetadataAndBytesAreAsExpected, since the sums are out of scope.
Metadata meaning is unchanged: the sums are still FP32 over the input
values. Making them exact integer sums over the quantized bytes is a
storage-contract change that has to move together with every kernel that
reads them, and stays in #1011.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
…nput
QuantPreprocessor::quantize could execute undefined behaviour for input the
API accepts. Three paths, all ending at the same conversion to a byte, and
all reachable with entirely finite components.
The scale was derived in FP32, so for [-FLT_MAX, +FLT_MAX] the range
overflowed: max - min became inf, delta inf, inv_delta 0, and the
per-element product inf * 0 = NaN, whose conversion to an integer is
undefined. UBSan: "-nan is outside the range of representable values of
type 'unsigned char'". The existing diff == 0 guard covers equal values,
not overflow of the subtraction. (MOD-17528)
With WithNorm, centering is an FP32 subtraction, so finite input against a
finite mean reaches inf before any range is computed: FLT_MAX against mean
-FLT_MAX centers to 6.8e38. Widening the range does not help, because the
value is already lost upstream.
And delta is stored as FP32, so (float)(diff / 255) underflows to zero
while diff itself is nonzero, which leaves 1/delta inf and scales the
minimum element, numerator exactly zero, to 0 * inf.
find_min_max now guarantees that both endpoints are finite and ordered, and
everything else follows from that. The guarantee lives there rather than at
the call site because that is where it can be broken: the WithNorm branch
creates the inf itself, from two valid operands, and the plain branch passes
through whatever the input holds. Both endpoints get a two-sided clamp,
since for an all-+inf vector inf <= inf passes the order check and a
one-sided std::max would leave min at +inf. Clamping to the float range is
not merely defensive either: min is stored as an FP32 field, so a non-finite
endpoint could not be represented under any arithmetic.
With that established, quantize needs one delta comparison, and inv_delta
stays double. Not for precision, which needs only +/-0.5 in 255, but because
an FP32 reciprocal overflows for a subnormal delta: [0, 7e-37] gives delta
2.7e-39, whose FP64 reciprocal is finite and correctly maps the top element
to 255.
The per-element bound and rounding are written by hand. That replaces
std::round, an out-of-line libm call at this translation unit's baseline
that ran once per element: bounding first makes adding 0.5 and truncating
equivalent for non-negative values. Measured at -O3 it is 8 instructions
against 9 for std::clamp plus std::round's call.
Scope, stated precisely, because it is narrower than it might look:
* The byte conversion is defined for all finite components. That is the
goal and it is met.
* The stored min and delta are finite with delta positive. The *sums* are
not covered: they are accumulated in FP32 over the input values, so
[-FLT_MAX, +FLT_MAX] stores sum_squares as inf even though this
function's own arithmetic is now well defined. Separate problem,
separate change.
* Non-finite components are unsupported, and nothing here is tested
against them. std::minmax_element requires its comparison to induce a
strict weak ordering, and floating-point < is not one once a NaN is
present: incomparability must be transitive, yet 1.0 is incomparable
with NaN and NaN with 2.0 while 1.0 and 2.0 are comparable. Violating
that precondition is undefined behaviour inside the algorithm, before
any range exists, so no assertion about the outcome would be portable,
including a weak one about metadata finiteness. UBSan reporting nothing
does not establish otherwise. Rejecting non-finite components at the
public ingestion boundary is the actual fix and belongs in its own
change; nothing in VecSim does it today.
The order check on the endpoints stays, as an invariant guard rather than a
NaN policy: everything downstream is written assuming min <= max, and
asserting that once per vector costs less than reasoning about whether it
can be violated.
This matches how comparable systems handle it. Lucene validates vector
components and throws on NaN or infinity, and Elasticsearch rejects NaN,
infinity and magnitudes that overflow before delegating to Lucene. Faiss
guards only an exactly zero range and assumes finite input otherwise, with
the same float-to-integer concern at its final cast. Qdrant gets defined
bytes from Rust's saturating cast, which C++ does not have, and can still
store an unusable scale.
Tests: the MOD-17528 reproduction, the WithNorm centering case, and a
table-driven matrix over finite input only: constant vectors positive
negative and zero, a single element, a subnormal but representable delta, a
range that underflows and collapses, and the full FP32 range. Expected bytes
and scale metadata are asserted rather than a range check, which is vacuous
for uint8_t, and every expectation was derived by simulating the pipeline
rather than predicted. The parameterized test is named for what it checks,
ScaleMetadataAndBytesAreAsExpected, since the sums are out of scope.
Metadata meaning is unchanged: the sums are still FP32 over the input
values. Making them exact integer sums over the quantized bytes is a
storage-contract change that has to move together with every kernel that
reads them, and stays in #1011.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
…nput
QuantPreprocessor::quantize could execute undefined behaviour for input the
API accepts. Three paths, all ending at the same conversion to a byte, and
all reachable with entirely finite components.
The scale was derived in FP32, so for [-FLT_MAX, +FLT_MAX] the range
overflowed: max - min became inf, delta inf, inv_delta 0, and the
per-element product inf * 0 = NaN, whose conversion to an integer is
undefined. UBSan: "-nan is outside the range of representable values of
type 'unsigned char'". The existing diff == 0 guard covers equal values,
not overflow of the subtraction. (MOD-17528)
With WithNorm, centering is an FP32 subtraction, so finite input against a
finite mean reaches inf before any range is computed: FLT_MAX against mean
-FLT_MAX centers to 6.8e38. Widening the range does not help, because the
value is already lost upstream.
And delta is stored as FP32, so (float)(diff / 255) underflows to zero
while diff itself is nonzero, which leaves 1/delta inf and scales the
minimum element, numerator exactly zero, to 0 * inf.
find_min_max now guarantees that both endpoints are finite and ordered, and
everything else follows from that. The guarantee lives there rather than at
the call site because that is where it can be broken: the WithNorm branch
creates the inf itself, from two valid operands, and the plain branch passes
through whatever the input holds. Both endpoints get a two-sided clamp,
since for an all-+inf vector inf <= inf passes the order check and a
one-sided std::max would leave min at +inf. Clamping to the float range is
not merely defensive either: min is stored as an FP32 field, so a non-finite
endpoint could not be represented under any arithmetic.
With that established, quantize needs one delta comparison, and inv_delta
stays double. Not for precision, which needs only +/-0.5 in 255, but because
an FP32 reciprocal overflows for a subnormal delta: [0, 7e-37] gives delta
2.7e-39, whose FP64 reciprocal is finite and correctly maps the top element
to 255.
The per-element bound and rounding are written by hand. That replaces
std::round, an out-of-line libm call at this translation unit's baseline
that ran once per element: bounding first makes adding 0.5 and truncating
equivalent for non-negative values. Measured at -O3 it is 8 instructions
against 9 for std::clamp plus std::round's call.
Scope, stated precisely, because it is narrower than it might look:
* The byte conversion is defined for all finite components. That is the
goal and it is met.
* The stored min and delta are finite with delta positive. The *sums* are
not covered: they are accumulated in FP32 over the input values, so
[-FLT_MAX, +FLT_MAX] stores sum_squares as inf even though this
function's own arithmetic is now well defined. Separate problem,
separate change.
* Non-finite components are unsupported, and nothing here is tested
against them. std::minmax_element requires its comparison to induce a
strict weak ordering, and floating-point < is not one once a NaN is
present: incomparability must be transitive, yet 1.0 is incomparable
with NaN and NaN with 2.0 while 1.0 and 2.0 are comparable. Violating
that precondition is undefined behaviour inside the algorithm, before
any range exists, so no assertion about the outcome would be portable,
including a weak one about metadata finiteness. UBSan reporting nothing
does not establish otherwise. Rejecting non-finite components at the
public ingestion boundary is the actual fix and belongs in its own
change; nothing in VecSim does it today.
The order check on the endpoints stays, as an invariant guard rather than a
NaN policy: everything downstream is written assuming min <= max, and
asserting that once per vector costs less than reasoning about whether it
can be violated.
This matches how comparable systems handle it. Lucene validates vector
components and throws on NaN or infinity, and Elasticsearch rejects NaN,
infinity and magnitudes that overflow before delegating to Lucene. Faiss
guards only an exactly zero range and assumes finite input otherwise, with
the same float-to-integer concern at its final cast. Qdrant gets defined
bytes from Rust's saturating cast, which C++ does not have, and can still
store an unusable scale.
Tests: the MOD-17528 reproduction, the WithNorm centering case, and a
table-driven matrix over finite input only: constant vectors positive
negative and zero, a single element, a subnormal but representable delta, a
range that underflows and collapses, and the full FP32 range. Expected bytes
and scale metadata are asserted rather than a range check, which is vacuous
for uint8_t, and every expectation was derived by simulating the pipeline
rather than predicted. The parameterized test is named for what it checks,
ScaleMetadataAndBytesAreAsExpected, since the sums are out of scope.
Metadata meaning is unchanged: the sums are still FP32 over the input
values. Making them exact integer sums over the quantized bytes is a
storage-contract change that has to move together with every kernel that
reads them, and stays in #1011.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

First of two stacked PRs closing MOD-17526. Found while reviewing #1007, which is what first makes any of this reachable: on
mainnothing constructs an index that usesQuantPreprocessor, 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 longercancels 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
Stored sums described the wrong vector.
QuantPreprocessoraccumulatedx_sumandx_sum_squaresover the original input, while every kernel that consumed them computed its other terms from the reconstructionmin + 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.The symmetric L2 kernels cancelled away the answer. They used
||x||^2 + ||y||^2 - 2*IPin 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.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.The plain uint8 kernels kept the same signed accumulator, on the paths MOD-17527 did not reach.
L2_NEON_UINT8.hassignedvaddvq_u32to anint32_t,L2_AVX512F_BW_VL_VNNI_UINT8.hreturned_mm512_reduce_add_epi32as a signedint, andret_tfor a 1-byte element type wasintin bothIP.cppandL2.cpp, which is signed-overflow UB from dimension 33,026 while its own comment claimed support to 2^16.L2_NEON_DOTPROD_UINT8.handL2_SVE_UINT8.hwere 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.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
L2Sqrshould not return a negative number.Which issues this PR fixes
||x||^2 + ||y||^2 - 2*IPidentity in FP32, and stored sums that describe the input rather than the reconstruction.QuantPreprocessorconverts NaN to uint8 when the derived range is not representable, which is undefined behaviour reachable fromAddVectorfor input the API accepts.Main objects this PR modified
vecsim_types::sq8(src/VecSim/types/sq8.h): metadata slots redefined asQ_SUMandQ_SUM_SQUARES, plus unaligned readers and the fourreconstructed_*derivation helpers.QuantPreprocessor(src/VecSim/spaces/computer/preprocessors.h): sums the quantized bytes, computes the range in double, and drops a libm call per element.IP.cppandL2.cpp.UINT8_InnerProductImpon all four ISAs, the L2 reduces inL2_AVX512F_BW_VL_VNNI_UINT8.handL2_NEON_UINT8.h, andret_tinIP.cppandL2.cpp.tests/unit/test_spaces.cpp,tests/unit/test_components.cppand the blob-building helpers intests/utils/tests_utils.handtests/unit/unit_test_utils.h, plus the uint8 spaces benchmark fixture.Mark if applicable
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
mainconstructs an index that usesQuantPreprocessor, so no persisted index can contain an SQ8 blob. SQ8 index creation arrives with #1007, which is stacked on this.What this changes
sq8.hredefines two metadata slots asQ_SUMandQ_SUM_SQUARES, uint32 sums over the quantized bytes. Blob size and slot count are unchanged: for an L2 index theSUMslot 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 ofstorage_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:
S1 + S2 - 2Qissum((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 bydelta^2and 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 clean0, hiding the very defect this PR fixes.The asymmetric L2 kernels derive
||x||^2exactly, which removes the negative distances. Their remaining cancellation needs the difference taken inside the loop and is PR 2 in this stack.UINT8_InnerProductImpreturns the exact integer dot on all four ISAs. Callers that wanted a float now convert explicitly, including one that computed1 - imp()in int arithmetic. On the plain uint8 side,ret_tis now 64-bit for every element type and stays signed, so1 - imp()cannot underflow, and the two L2 reduces are read back unsigned, which costs no instructions on either ISA.MAX_EXACT_DIMis 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, whereasQ_SUM_SQUARESstays exact in uint32 through dimension 66,051. For the record the real per-field ceilings are 66,051 forQ_SUM_SQUARES(L2 only), 16,843,009 forQ_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.
SQ8_SQ8_NAIVE_*scalar kernelsQuantPreprocessor::quantize)Two mitigations are included.
sq8::reconstructed_*are markedalways_inline, because GCC outlinedreconstructed_l2_sqrand then built a realigned 64-byte stack frame plus avzeroupperinside 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::roundin the quantizer was an out-of-line call to libmround()per element, because that translation unit compiles at the x86-64 baseline and has noroundsd; since the value is already clamped to[0, 255]and therefore non-negative, adding 0.5 and truncating isstd::roundand compiles to a singlecvttsd2si. That turns an 18% regression on the quantize path into a 26% improvement overmain.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:test_spacestest_componentsctestunit suitecheck-format.shThe six later commits (uint8 accumulators, the integer regrouping, both performance changes,
MAX_EXACT_DIMremoval, the new tests) areclang-formatclean 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 == maxcase now expects both sums to be 0, since every value collapses to byte 0, and additionally assertsreconstructed_sumstill recovers3.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.
SQ8_SQ8_L2_survives_large_common_offsetSQ8_FP32_L2_is_non_negative_for_off_grid_querySQ8_SQ8_L2_matches_fp64_reference_at_high_dimQuantizationHandlesNonRepresentableRangeSix 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::DELTAand the public kernels, so the same source compiles on either side: a test written againstQ_SUMorsq8::reconstructed_*would not compile before the change and so could not demonstrate anything.SQ8_SQ8_L2_self_distance_is_exactly_zeroSQ8_SQ8_L2_self_distance_whiteboard_case[0, 100.5, 255]against itself, with numbers small enough to check by handSQ8_SQ8_distance_ignores_sub_quantum_input_changesSQ8_SQ8_IP_matches_reconstruction_dotUINT8_L2Sqr_and_InnerProduct_are_exact_past_int32Dimensions are chosen so both kernel families run: 3, 5 and 8 dispatch to the naive implementations and 64 and above to
AVX512F_BW_VL_VNNIon an Ice Lake host. The identical-bytes test keepsminnonzero on purpose, because withmin == 0the old IP formula collapses todelta_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[]withdeleteand stored the trailing norms through unalignedfloatcasts, 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]mademax - minoverflow to inf, thendeltainf,inv_delta0, andinf * 0NaN, whose conversion to an integer is UB, reachable fromAddVectorfor 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_deltawas 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::roundreplacement: the two forms differ only wherex + 0.5itself double-rounds, atscaled == 0.49999999999999994, where the new form yields byte 1 andstd::roundyields 0. There is no divergence at any of the 255.5boundaries, across 2M uniform random values, or across 3M realistic(x - min) * inv_deltadraws.lrintandnearbyintare not substitutes: they round ties to even, which would move 178.5 to 178 and change bytes in common cases.Not in this PR
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_SQUARESover bytes, same blob layout) and routes distance math through sharedsq8::reconstructed_*helpers in double instead of duplicated inline formulas.QuantPreprocessornow 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_tbyte dot products andreconstructed_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||²viareconstructed_sum_squaresand 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.