[MOD-17526] Make the SQ8 metadata describe the reconstruction

src/VecSim/spaces/computer/preprocessors.h

@@ -152,9 +152,14 @@ class CosinePreprocessor : public PreprocessorInterface {
  *
  * Storage layout:
  * | quantized_values[dim] | min_val | delta | x_sum | (x_sum_squares for L2 only) |
- * where:
- * x_sum = Σx_i: sum of the original values,
- * x_sum_squares = Σx_i²: sum of squares of the original values.
+ * where, writing x_r[i] = min_val + delta * a[i] for the value a blob actually represents:
+ * x_sum        = Σx_r[i]:  sum of the reconstructed values,
+ * x_sum_squares = Σx_r[i]²: sum of their squares.
+ *
+ * These describe x_r, not the input x. Every formula below is written in terms of x_r, because
+ * that is what a quantized blob holds, so sums over the input would not satisfy them: the gap is
+ * the quantization error, about 0.4% of ||x||², which is larger than the distance between two
+ * similar vectors and made L2 come out negative.
  *
  * Storage metadata is always FP32 (independent of DataType) to match the asymmetric distance
  * kernels. The quantized blob size is:
@@ -194,10 +199,11 @@ class CosinePreprocessor : public PreprocessorInterface {
  *   where y_sum = Σy_i is precomputed and stored in the query blob.
  *
  * For L2:
- *   ||x - y||² = Σx_i² - 2*Σ(x_i * y_i) + Σy_i²
- *              = x_sum_squares - 2 * IP(x, y) + y_sum_squares
+ *   ||x_r - y||² = Σx_r[i]² - 2*Σ(x_r[i] * y_i) + Σy_i²
+ *                = x_sum_squares - 2 * IP(x_r, y) + y_sum_squares
  *   where:
- *     - x_sum_squares = Σx_i² is precomputed and stored in the storage blob
+ *     - x_sum_squares = Σx_r[i]² is precomputed and stored in the storage blob. A query is not
+ *       quantized, so y_sum_squares is Σy_i² over the query itself.
  *     - IP(x, y) is computed using the formula above
  *     - y_sum_squares = Σy_i² is precomputed and stored in the query blob
  *   For normalized L2, x and y in this formula are the centered values x' and y'; their distance
@@ -217,7 +223,7 @@ class CosinePreprocessor : public PreprocessorInterface {
  *              + delta_x * delta_y * Σ(qx_i * qy_i)
  *   where:
  *     - sum_x, sum_y are precomputed sums of the values represented by each blob
- *     - Σqx_i = (sum_x - dim * min_x) / delta_x  (sum of quantized values, derived from stored sum)
+ *     - Σqx_i = (sum_x - dim * min_x) / delta_x  (exact, since sum_x is derived from Σqx_i)
  *     - Σqy_i = (sum_y - dim * min_y) / delta_y
  *
  * For L2:
@@ -293,13 +299,19 @@ class QuantPreprocessor : public PreprocessorInterface {
             return static_cast<OUTPUT_TYPE>(scaled + MetadataType{0.5});
         };
 
-        // Compute sum (and sum of squares for L2) while quantizing.
-        // Accumulators are FP32 to preserve metadata precision for FP16 inputs.
-        // 4 independent accumulators (sum)
-        float s0{}, s1{}, s2{}, s3{};
-
-        // 4 independent accumulators (sum of squares), only used for L2
-        float q0{}, q1{}, q2{}, q3{};
+        // Sum the quantized bytes, not the input values. Every kernel term is written in terms of
+        // the reconstruction x_r[i] = min + delta * a[i], so the stored sums have to describe x_r
+        // as well. Summing the input instead leaves a systematic mismatch of about 0.4% of ||x||^2,
+        // which is larger than a small true distance and made L2 come out negative: for two
+        // near-duplicate vectors at dim 128 the reconstruction distance is 2.93e-04 and the kernels
+        // returned -1.49e-01. The byte sums are exact integers, so the derivation below is the only
+        // place rounding enters.
+        //
+        // 4 independent accumulators each, so the unrolled loop keeps four dependency chains.
+        uint32_t s0{}, s1{}, s2{}, s3{};
+        // 64-bit: 65025 * 65536 leaves under 1% of UINT32_MAX, and overflow here would corrupt
+        // the metadata rather than round it. This is the write path, once per vector.
+        uint64_t q0{}, q1{}, q2{}, q3{}; // only used for L2
 
         size_t i = 0;
         // round dim down to the nearest multiple of 4
@@ -315,40 +327,59 @@ class QuantPreprocessor : public PreprocessorInterface {
             // We know (input - min) => 0
             // If min == max, all values are the same and should be quantized to 0.
             // reconstruction will yield the same original value for all vectors.
-            quantized[i] = to_byte((x0 - min_val) * inv_delta);
-            quantized[i + 1] = to_byte((x1 - min_val) * inv_delta);
-            quantized[i + 2] = to_byte((x2 - min_val) * inv_delta);
-            quantized[i + 3] = to_byte((x3 - min_val) * inv_delta);
-
-            // Accumulate sum for all metrics
-            s0 += x0;
-            s1 += x1;
-            s2 += x2;
-            s3 += x3;
-
-            // Accumulate sum of squares only for L2 metric
+            const OUTPUT_TYPE a0 = to_byte((x0 - min_val) * inv_delta);
+            const OUTPUT_TYPE a1 = to_byte((x1 - min_val) * inv_delta);
+            const OUTPUT_TYPE a2 = to_byte((x2 - min_val) * inv_delta);
+            const OUTPUT_TYPE a3 = to_byte((x3 - min_val) * inv_delta);
+            quantized[i] = a0;
+            quantized[i + 1] = a1;
+            quantized[i + 2] = a2;
+            quantized[i + 3] = a3;
+
+            s0 += a0;
+            s1 += a1;
+            s2 += a2;
+            s3 += a3;
+
             if constexpr (Metric == VecSimMetric_L2) {
-                q0 += x0 * x0;
-                q1 += x1 * x1;
-                q2 += x2 * x2;
-                q3 += x3 * x3;
+                q0 += uint64_t{a0} * a0;
+                q1 += uint64_t{a1} * a1;
+                q2 += uint64_t{a2} * a2;
+                q3 += uint64_t{a3} * a3;
             }
         }
 
         // Tail: 0..3 remaining elements (still the same pass, just finishing work).
-        // Sum/sum_squares become metadata, so they are MetadataType.
-        MetadataType sum = (s0 + s1) + (s2 + s3);
-        MetadataType sum_squares = (q0 + q1) + (q2 + q3);
+        uint32_t q_sum = (s0 + s1) + (s2 + s3);
+        uint64_t q_sum_squares{};
+        if constexpr (Metric == VecSimMetric_L2) {
+            q_sum_squares = (q0 + q1) + (q2 + q3);
+        }
 
         for (; i < this->dim; ++i) {
             const float x = transformed_value(input, i);
-            quantized[i] = to_byte((x - min_val) * inv_delta);
-            sum += x;
+            const OUTPUT_TYPE a = to_byte((x - min_val) * inv_delta);
+            quantized[i] = a;
+            q_sum += a;
             if constexpr (Metric == VecSimMetric_L2) {
-                sum_squares += x * x;
+                q_sum_squares += uint64_t{a} * a;
             }
         }
 
+        // Derive the reconstruction sums from the exact byte sums, in double so the expansion does
+        // not lose the cross term, then store FP32 as before. The slot layout and types are
+        // unchanged; only what the numbers describe changes.
+        //   sum         = sum(x_r[i])   = dim*min + delta*q_sum
+        //   sum_squares = sum(x_r[i]^2) = dim*min^2 + 2*min*delta*q_sum + delta^2*q_sum_squares
+        const double d_min = min_val, d_delta = delta, d_dim = static_cast<double>(this->dim);
+        const MetadataType sum = static_cast<MetadataType>(d_dim * d_min + d_delta * q_sum);
+        MetadataType sum_squares{};
+        if constexpr (Metric == VecSimMetric_L2) {
+            sum_squares =
+                static_cast<MetadataType>(d_dim * d_min * d_min + 2.0 * d_min * d_delta * q_sum +
+                                          d_delta * d_delta * q_sum_squares);
+        }
+
         // Metadata uses MetadataType. Use memcpy because the metadata offset
         // (dim * sizeof(uint8_t)) is not guaranteed to be sizeof(MetadataType)-aligned.
         void *meta_dst = quantized + this->dim;

tests/unit/test_components.cpp

@@ -1650,6 +1650,14 @@ class QuantPreprocessorFP16MetricTest : public testing::TestWithParam<VecSimMetr
         ComputeSQ8Quantization(widened_blob, dim, baseline_storage);
         const float baseline_min = load_meta(baseline_storage, dim + sq8::MIN_VAL * sizeof(float));
         const float baseline_delta = load_meta(baseline_storage, dim + sq8::DELTA * sizeof(float));
+        // The storage baseline's sums describe the reconstruction min + delta * a[i]. The query is
+        // not quantized, so its sums describe the input itself. Those are different quantities and
+        // must be checked separately: they were only ever equal because both used to be the input.
+        float input_sum = 0.0f, input_sum_sq = 0.0f;
+        for (size_t i = 0; i < dim; i++) {
+            input_sum += widened_blob[i];
+            input_sum_sq += widened_blob[i] * widened_blob[i];
+        }
         const float baseline_sum = load_meta(baseline_storage, dim + sq8::SUM * sizeof(float));
         const float baseline_sum_sq =
             load_meta(baseline_storage, dim + sq8::SUM_SQUARES * sizeof(float));
@@ -1702,11 +1710,11 @@ class QuantPreprocessorFP16MetricTest : public testing::TestWithParam<VecSimMetr
             // Query FP32 metadata: y_sum (and y_sum_squares for L2) match the FP32 baseline.
             ASSERT_FLOAT_EQ(
                 load_meta(query_blob, query_meta_offset + sq8::SUM_QUERY * sizeof(float)),
-                baseline_sum);
+                input_sum);
             if constexpr (Metric == VecSimMetric_L2) {
                 ASSERT_FLOAT_EQ(load_meta(query_blob, query_meta_offset +
                                                           sq8::SUM_SQUARES_QUERY * sizeof(float)),
-                                baseline_sum_sq);
+                                input_sum_sq);
             }
 
             allocator->free_allocation(storage_blob);
@@ -1724,11 +1732,11 @@ class QuantPreprocessorFP16MetricTest : public testing::TestWithParam<VecSimMetr
             EXPECT_NO_FATAL_FAILURE(
                 CompareVectors<float16>(static_cast<const float16 *>(blob), original_blob, dim));
             ASSERT_FLOAT_EQ(load_meta(blob, query_meta_offset + sq8::SUM_QUERY * sizeof(float)),
-                            baseline_sum);
+                            input_sum);
             if constexpr (Metric == VecSimMetric_L2) {
                 ASSERT_FLOAT_EQ(
                     load_meta(blob, query_meta_offset + sq8::SUM_SQUARES_QUERY * sizeof(float)),
-                    baseline_sum_sq);
+                    input_sum_sq);
             }
             allocator->free_allocation(blob);
         }

tests/unit/test_spaces.cpp

@@ -4908,3 +4908,154 @@ TEST_F(SpacesTest, SQ8_SQ8_DispatcherAlignmentHints) {
     check("Cosine", &spaces::Cosine_SQ8_SQ8_GetDistFunc);
 }
 #endif // CPU_FEATURES_ARCH_X86_64
+
+namespace {
+// The stored sum_squares is ||x_r||^2 by construction, so read it rather than recomputing.
+float sq8_norm_sq(const std::vector<uint8_t> &blob, size_t dim) {
+    return load_unaligned<float>(blob.data() + dim + sq8::SUM_SQUARES * sizeof(float));
+}
+
+std::vector<uint8_t> sq8_quantize_l2(const std::vector<float> &v) {
+    const size_t dim = v.size();
+    std::vector<uint8_t> blob(dim * sizeof(uint8_t) + 4 * sizeof(float));
+    test_utils::quantize_float_vec_to_sq8_with_metadata(v.data(), dim, blob.data());
+    return blob;
+}
+} // namespace
+
+namespace {
+// The reconstruction the kernels compute distances over: x_r[i] = min + delta * a[i]. Reference
+// values are accumulated in double so the comparison is against the algebra, not against another
+// FP32 implementation with the same rounding.
+double ReferenceL2SqrOverReconstruction(const uint8_t *a, const uint8_t *b, size_t dim) {
+    const float min_a = load_unaligned<float>(a + dim + sq8::MIN_VAL * sizeof(float));
+    const float delta_a = load_unaligned<float>(a + dim + sq8::DELTA * sizeof(float));
+    const float min_b = load_unaligned<float>(b + dim + sq8::MIN_VAL * sizeof(float));
+    const float delta_b = load_unaligned<float>(b + dim + sq8::DELTA * sizeof(float));
+    double acc = 0.0;
+    for (size_t i = 0; i < dim; i++) {
+        const double d = (static_cast<double>(min_a) + static_cast<double>(delta_a) * a[i]) -
+                         (static_cast<double>(min_b) + static_cast<double>(delta_b) * b[i]);
+        acc += d * d;
+    }
+    return acc;
+}
+} // namespace
+
+// The stored sums describe the reconstruction, not the input. Summing the input instead left a
+// systematic mismatch of about 0.4% of ||x||^2, which is larger than the distance between two
+// near-duplicate vectors, so L2 came back negative: at dim 128 the reconstruction distance is
+// 2.93e-04 and the kernels returned -1.49e-01. Near-duplicates are the case that exposes it,
+// because the true distance is small enough for the mismatch to dominate.
+TEST_F(SpacesTest, SQ8_SQ8_L2_is_non_negative_and_matches_reconstruction) {
+    std::mt19937 gen(20260820);
+    std::uniform_real_distribution<float> value_gen(-1.0f, 1.0f);
+    std::uniform_real_distribution<float> nudge_gen(-1e-3f, 1e-3f);
+
+    for (const size_t dim : {4UL, 15UL, 64UL, 128UL, 512UL}) {
+        std::vector<float> x(dim), y(dim);
+        for (size_t i = 0; i < dim; i++) {
+            x[i] = value_gen(gen);
+            y[i] = x[i] + nudge_gen(gen);
+        }
+        const std::vector<uint8_t> qx = sq8_quantize_l2(x);
+        const std::vector<uint8_t> qy = sq8_quantize_l2(y);
+
+        const double expected = ReferenceL2SqrOverReconstruction(qx.data(), qy.data(), dim);
+
+        // The error floor is FP32 cancellation in sum_sq_x + sum_sq_y - 2*IP, so it scales with the
+        // norms, not with the distance. A flat bound is thousands of times the distance it guards
+        // at small dimensions, which would let a kernel returning zero pass.
+        const float tol = 8.0f * std::numeric_limits<float>::epsilon() *
+                          (sq8_norm_sq(qx, dim) + sq8_norm_sq(qy, dim));
+        unsigned char alignment = 0;
+        auto dispatched = L2_SQ8_SQ8_GetDistFunc(dim, &alignment, nullptr);
+
+        for (const auto &probe :
+             {std::make_pair("scalar", SQ8_SQ8_L2Sqr), std::make_pair("dispatched", dispatched)}) {
+            const float got = probe.second(qx.data(), qy.data(), dim);
+            // Bounded by the noise floor, not by zero: this change does not make the result
+            // provably non-negative. The two sides of the subtraction are computed by different
+            // routes, and a blob against itself lands just below zero at dim 512.
+            EXPECT_GE(got, -tol) << probe.first << " negative beyond the noise floor, dim " << dim;
+            EXPECT_NEAR(got, expected, tol)
+                << probe.first << " does not match the reconstruction, dim " << dim;
+        }
+    }
+}
+
+// The asymmetric path is how the defect reached production unnoticed: its tests run dimensions 1,
+// 5, 7 and 15 against an absolute tolerance of 0.01, and 0.4% of ||x||^2 at dim 5 is about 0.007,
+// which fits underneath. Storage is quantized and the query is not, so this checks
+// sum(x_r^2) + sum(y^2) - 2*sum(x_r*y) against a double reference over the same two sides.
+TEST_F(SpacesTest, SQ8_FP32_L2_matches_reconstruction_against_float_query) {
+    std::mt19937 gen(20260822);
+    std::uniform_real_distribution<float> value_gen(-1.0f, 1.0f);
+    std::uniform_real_distribution<float> nudge_gen(-1e-3f, 1e-3f);
+
+    for (const size_t dim : {5UL, 15UL, 64UL, 128UL, 512UL}) {
+        std::vector<float> x(dim), y(dim);
+        for (size_t i = 0; i < dim; i++) {
+            x[i] = value_gen(gen);
+            y[i] = x[i] + nudge_gen(gen); // near-duplicate: small true distance exposes the offset
+        }
+        const std::vector<uint8_t> qx = sq8_quantize_l2(x);
+
+        // The query carries its own sum of squares, over the query itself, since it is not
+        // quantized. Reference is the distance between the reconstruction and that query.
+        std::vector<float> query(dim + 2, 0.0f);
+        double ref = 0.0, y_sum_sq = 0.0;
+        const float min_x = load_unaligned<float>(qx.data() + dim + sq8::MIN_VAL * sizeof(float));
+        const float delta_x = load_unaligned<float>(qx.data() + dim + sq8::DELTA * sizeof(float));
+        for (size_t i = 0; i < dim; i++) {
+            query[i] = y[i];
+            y_sum_sq += static_cast<double>(y[i]) * y[i];
+            const double d = (static_cast<double>(min_x) + static_cast<double>(delta_x) * qx[i]) -
+                             static_cast<double>(y[i]);
+            ref += d * d;
+        }
+        // Both query slots matter: the asymmetric inner product is
+        // min*sum(y) + delta*sum(a[i]*y[i]), so zeroing SUM_QUERY silently drops the min*sum(y)
+        // term. Fill them the way the production query preprocessor does.
+        test_utils::preprocess_sq8_fp32_query(query.data(), dim);
+
+        const float tol = 8.0f * std::numeric_limits<float>::epsilon() *
+                          (sq8_norm_sq(qx, dim) + static_cast<float>(y_sum_sq));
+        unsigned char alignment = 0;
+        auto dispatched = L2_SQ8_FP32_GetDistFunc(dim, &alignment, nullptr);
+        for (const auto &probe :
+             {std::make_pair("scalar", SQ8_FP32_L2Sqr), std::make_pair("dispatched", dispatched)}) {
+            // storage first, query second, as the other asymmetric tests call it.
+            const float got = probe.second(qx.data(), query.data(), dim);
+            EXPECT_GE(got, -tol) << probe.first << " negative beyond the noise floor, dim " << dim;
+            EXPECT_NEAR(got, ref, tol)
+                << probe.first << " does not match the reconstruction, dim " << dim;
+        }
+    }
+}
+
+// A blob against itself. Exact zero is not claimed here: the two sides of
+// sum_sq_x + sum_sq_y - 2*IP are computed by different routes and round differently, so this pins
+// the magnitude rather than the bit pattern.
+TEST_F(SpacesTest, SQ8_SQ8_L2_self_distance_is_near_zero) {
+    std::mt19937 gen(20260821);
+    std::uniform_real_distribution<float> value_gen(-3.0f, 5.0f);
+
+    for (const size_t dim : {3UL, 8UL, 64UL, 512UL}) {
+        std::vector<float> x(dim);
+        double norm_sq = 0.0;
+        for (size_t i = 0; i < dim; i++) {
+            x[i] = value_gen(gen);
+            norm_sq += static_cast<double>(x[i]) * x[i];
+        }
+        const std::vector<uint8_t> q = sq8_quantize_l2(x);
+        const float tol = static_cast<float>(1e-5 * norm_sq);
+
+        EXPECT_NEAR(SQ8_SQ8_L2Sqr(q.data(), q.data(), dim), 0.0f, tol)
+            << "scalar kernel, dim " << dim;
+        unsigned char alignment = 0;
+        auto dispatched = L2_SQ8_SQ8_GetDistFunc(dim, &alignment, nullptr);
+        EXPECT_NEAR(dispatched(q.data(), q.data(), dim), 0.0f, tol)
+            << "dispatched kernel, dim " << dim;
+    }
+}