diff --git a/src/VecSim/spaces/IP/IP.cpp b/src/VecSim/spaces/IP/IP.cpp index 2140c2345..b78ee4fe2 100644 --- a/src/VecSim/spaces/IP/IP.cpp +++ b/src/VecSim/spaces/IP/IP.cpp @@ -146,28 +146,18 @@ float SQ8_SQ8_InnerProduct_Impl(const void *pVect1v, const void *pVect2v, size_t const auto *pVect1 = static_cast(pVect1v); const auto *pVect2 = static_cast(pVect2v); - // Compute inner product of quantized values: Σ(q1[i]*q2[i]) - float product = 0; + // Inner product of the quantized bytes: Σ(a[i]*b[i]). Exact in integer arithmetic. + uint64_t q_dot = 0; for (size_t i = 0; i < dimension; i++) { - product += pVect1[i] * pVect2[i]; + q_dot += static_cast(pVect1[i]) * static_cast(pVect2[i]); } - // Metadata follows byte payloads and is not necessarily float-aligned. - const auto *params1 = pVect1 + dimension; - const float min_val1 = load_unaligned(params1 + sq8::MIN_VAL * sizeof(float)); - const float delta1 = load_unaligned(params1 + sq8::DELTA * sizeof(float)); - const float sum1 = load_unaligned(params1 + sq8::SUM * sizeof(float)); - - // Get quantization parameters from pVect2 - const auto *params2 = pVect2 + dimension; - const float min_val2 = load_unaligned(params2 + sq8::MIN_VAL * sizeof(float)); - const float delta2 = load_unaligned(params2 + sq8::DELTA * sizeof(float)); - const float sum2 = load_unaligned(params2 + sq8::SUM * sizeof(float)); - - // Apply the algebraic formula using precomputed sums: - // IP = min1*sum2 + min2*sum1 + delta1*delta2*Σ(q1[i]*q2[i]) - dim*min1*min2 - return min_val1 * sum2 + min_val2 * sum1 - static_cast(dimension) * min_val1 * min_val2 + - delta1 * delta2 * product; + // Both operands describe their reconstructions, so every term is derived from the stored + // integer sums rather than from sums over the original inputs. See sq8.h. + return static_cast(sq8::reconstructed_ip( + sq8::min_val_of(pVect1, dimension), sq8::delta_of(pVect1, dimension), + sq8::q_sum_of(pVect1, dimension), sq8::min_val_of(pVect2, dimension), + sq8::delta_of(pVect2, dimension), sq8::q_sum_of(pVect2, dimension), q_dot, dimension)); } // SQ8-to-SQ8: Both vectors are uint8 quantized with precomputed sum @@ -238,11 +228,12 @@ float FP16_InnerProduct(const void *pVect1, const void *pVect2, size_t dimension } // Return type for the inner product functions. -// The type should be able to hold `dimension * MAX_VAL(int_elem_t) * MAX_VAL(int_elem_t)`. -// To support dimension up to 2^16, we need the difference between the type and int_elem_t to be at -// least 2 bytes. We assert that in the implementation. +// The type must hold `dimension * MAX_VAL(int_elem_t) * MAX_VAL(int_elem_t)`. For uint8 that is +// 65025 * dimension, which overflows a 32-bit int from dimension 33,026, so this is 64-bit for +// every element type. Keeping it signed also means the `1 - ip` in the wrappers below stays +// signed arithmetic and cannot underflow. We assert the width in the implementation. template -using ret_t = std::conditional_t; +using ret_t = long long; template static inline ret_t @@ -273,7 +264,7 @@ float INT8_Cosine(const void *pVect1v, const void *pVect2v, size_t dimension) { float UINT8_InnerProduct(const void *pVect1v, const void *pVect2v, size_t dimension) { const auto *pVect1 = static_cast(pVect1v); const auto *pVect2 = static_cast(pVect2v); - return 1 - INTEGER_InnerProductImp(pVect1, pVect2, dimension); + return 1.0f - static_cast(INTEGER_InnerProductImp(pVect1, pVect2, dimension)); } float UINT8_Cosine(const void *pVect1v, const void *pVect2v, size_t dimension) { diff --git a/src/VecSim/spaces/IP/IP_AVX512F_BW_VL_VNNI_SQ8_SQ8.h b/src/VecSim/spaces/IP/IP_AVX512F_BW_VL_VNNI_SQ8_SQ8.h index ae6f96ea2..58480dcf5 100644 --- a/src/VecSim/spaces/IP/IP_AVX512F_BW_VL_VNNI_SQ8_SQ8.h +++ b/src/VecSim/spaces/IP/IP_AVX512F_BW_VL_VNNI_SQ8_SQ8.h @@ -38,29 +38,16 @@ using sq8 = vecsim_types::sq8; // Uses UINT8_InnerProductImp for efficient dot product computation with VNNI template // 0..63 float SQ8_SQ8_InnerProductImp(const void *pVec1v, const void *pVec2v, size_t dimension) { - // Compute raw dot product using efficient UINT8 AVX512 VNNI implementation - // UINT8_InnerProductImp uses _mm512_dpwssd_epi32 for native integer dot product - int dot_product = UINT8_InnerProductImp(pVec1v, pVec2v, dimension); + // Exact integer dot product of the quantized bytes. + const uint64_t q_dot = UINT8_InnerProductImp(pVec1v, pVec2v, dimension); - // Get dequantization parameters and precomputed values from the end of vectors - // Layout: [data (dim)] [min (float)] [delta (float)] [sum (float)] + // Every term describes the reconstructions, derived from the stored integer sums. See sq8.h. const uint8_t *pVec1 = static_cast(pVec1v); const uint8_t *pVec2 = static_cast(pVec2v); - - const auto *params1 = pVec1 + dimension; - const float min1 = load_unaligned(params1 + sq8::MIN_VAL * sizeof(float)); - const float delta1 = load_unaligned(params1 + sq8::DELTA * sizeof(float)); - const float sum1 = load_unaligned(params1 + sq8::SUM * sizeof(float)); - - const auto *params2 = pVec2 + dimension; - const float min2 = load_unaligned(params2 + sq8::MIN_VAL * sizeof(float)); - const float delta2 = load_unaligned(params2 + sq8::DELTA * sizeof(float)); - const float sum2 = load_unaligned(params2 + sq8::SUM * sizeof(float)); - - // Apply the algebraic formula using precomputed sums: - // IP = min1*sum2 + min2*sum1 + δ1*δ2 * Σ(q1[i]*q2[i]) - dim*min1*min2 - return min1 * sum2 + min2 * sum1 + delta1 * delta2 * static_cast(dot_product) - - static_cast(dimension) * min1 * min2; + return static_cast(sq8::reconstructed_ip( + sq8::min_val_of(pVec1, dimension), sq8::delta_of(pVec1, dimension), + sq8::q_sum_of(pVec1, dimension), sq8::min_val_of(pVec2, dimension), + sq8::delta_of(pVec2, dimension), sq8::q_sum_of(pVec2, dimension), q_dot, dimension)); } // SQ8-to-SQ8 Inner Product distance function diff --git a/src/VecSim/spaces/IP/IP_AVX512F_BW_VL_VNNI_UINT8.h b/src/VecSim/spaces/IP/IP_AVX512F_BW_VL_VNNI_UINT8.h index bd43bc901..3d95e53c9 100644 --- a/src/VecSim/spaces/IP/IP_AVX512F_BW_VL_VNNI_UINT8.h +++ b/src/VecSim/spaces/IP/IP_AVX512F_BW_VL_VNNI_UINT8.h @@ -31,8 +31,8 @@ static inline void InnerProductStep(uint8_t *&pVect1, uint8_t *&pVect2, __m512i } template // 0..63 -static inline int UINT8_InnerProductImp(const void *pVect1v, const void *pVect2v, - size_t dimension) { +static inline uint64_t UINT8_InnerProductImp(const void *pVect1v, const void *pVect2v, + size_t dimension) { uint8_t *pVect1 = (uint8_t *)pVect1v; uint8_t *pVect2 = (uint8_t *)pVect2v; @@ -87,19 +87,24 @@ static inline int UINT8_InnerProductImp(const void *pVect1v, const void *pVect2v } while (pVect1 < pEnd1); } - return _mm512_reduce_add_epi32(sum); + // Reduce in 64-bit. Each int32 lane is in range, but their total reaches 255*255*dim, which + // passes INT_MAX from dimension 33,027, so a 32-bit horizontal sum would wrap. Callers rely + // on this being the exact integer dot product. + const __m512i lo = _mm512_cvtepu32_epi64(_mm512_extracti64x4_epi64(sum, 0)); + const __m512i hi = _mm512_cvtepu32_epi64(_mm512_extracti64x4_epi64(sum, 1)); + return static_cast(_mm512_reduce_add_epi64(_mm512_add_epi64(lo, hi))); } template // 0..63 float UINT8_InnerProductSIMD64_AVX512F_BW_VL_VNNI(const void *pVect1v, const void *pVect2v, size_t dimension) { - return 1 - UINT8_InnerProductImp(pVect1v, pVect2v, dimension); + return 1.0f - static_cast(UINT8_InnerProductImp(pVect1v, pVect2v, dimension)); } template // 0..63 float UINT8_CosineSIMD64_AVX512F_BW_VL_VNNI(const void *pVect1v, const void *pVect2v, size_t dimension) { - float ip = UINT8_InnerProductImp(pVect1v, pVect2v, dimension); + float ip = static_cast(UINT8_InnerProductImp(pVect1v, pVect2v, dimension)); const float norm_v1 = load_unaligned(static_cast(pVect1v) + dimension); const float norm_v2 = load_unaligned(static_cast(pVect2v) + dimension); return 1.0f - ip / (norm_v1 * norm_v2); diff --git a/src/VecSim/spaces/IP/IP_NEON_DOTPROD_SQ8_SQ8.h b/src/VecSim/spaces/IP/IP_NEON_DOTPROD_SQ8_SQ8.h index d7f5b444e..b16b3fc99 100644 --- a/src/VecSim/spaces/IP/IP_NEON_DOTPROD_SQ8_SQ8.h +++ b/src/VecSim/spaces/IP/IP_NEON_DOTPROD_SQ8_SQ8.h @@ -39,29 +39,16 @@ using sq8 = vecsim_types::sq8; template // 0..63 float SQ8_SQ8_InnerProductSIMD64_NEON_DOTPROD_IMP(const void *pVec1v, const void *pVec2v, size_t dimension) { - // Compute raw dot product using efficient UINT8 DOTPROD implementation - // UINT8_InnerProductImp uses vdotq_u32 for native uint8 dot product - float dot_product = UINT8_InnerProductImp(pVec1v, pVec2v, dimension); + // Exact integer dot product of the quantized bytes. + const uint64_t q_dot = UINT8_InnerProductImp(pVec1v, pVec2v, dimension); - // Get dequantization parameters and precomputed values from the end of vectors - // Layout: [data (dim)] [min (float)] [delta (float)] [sum (float)] + // Every term describes the reconstructions, derived from the stored integer sums. See sq8.h. const uint8_t *pVec1 = static_cast(pVec1v); const uint8_t *pVec2 = static_cast(pVec2v); - - const auto *params1 = pVec1 + dimension; - const float min1 = load_unaligned(params1 + sq8::MIN_VAL * sizeof(float)); - const float delta1 = load_unaligned(params1 + sq8::DELTA * sizeof(float)); - const float sum1 = load_unaligned(params1 + sq8::SUM * sizeof(float)); - - const auto *params2 = pVec2 + dimension; - const float min2 = load_unaligned(params2 + sq8::MIN_VAL * sizeof(float)); - const float delta2 = load_unaligned(params2 + sq8::DELTA * sizeof(float)); - const float sum2 = load_unaligned(params2 + sq8::SUM * sizeof(float)); - - // Apply algebraic formula using precomputed sums: - // IP = min1*sum2 + min2*sum1 + δ1*δ2 * Σ(q1*q2) - dim*min1*min2 - return min1 * sum2 + min2 * sum1 + delta1 * delta2 * dot_product - - static_cast(dimension) * min1 * min2; + return static_cast(sq8::reconstructed_ip( + sq8::min_val_of(pVec1, dimension), sq8::delta_of(pVec1, dimension), + sq8::q_sum_of(pVec1, dimension), sq8::min_val_of(pVec2, dimension), + sq8::delta_of(pVec2, dimension), sq8::q_sum_of(pVec2, dimension), q_dot, dimension)); } // SQ8-to-SQ8 Inner Product distance function diff --git a/src/VecSim/spaces/IP/IP_NEON_DOTPROD_UINT8.h b/src/VecSim/spaces/IP/IP_NEON_DOTPROD_UINT8.h index 3abbd9bba..71620f816 100644 --- a/src/VecSim/spaces/IP/IP_NEON_DOTPROD_UINT8.h +++ b/src/VecSim/spaces/IP/IP_NEON_DOTPROD_UINT8.h @@ -27,7 +27,7 @@ InnerProductStep(uint8_t *&pVect1, uint8_t *&pVect2, uint32x4_t &sum) { } template // 0..63 -float UINT8_InnerProductImp(const void *pVect1v, const void *pVect2v, size_t dimension) { +uint64_t UINT8_InnerProductImp(const void *pVect1v, const void *pVect2v, size_t dimension) { uint8_t *pVect1 = (uint8_t *)pVect1v; uint8_t *pVect2 = (uint8_t *)pVect2v; @@ -97,20 +97,20 @@ float UINT8_InnerProductImp(const void *pVect1v, const void *pVect2v, size_t dim uint32x4_t total_sum = vaddq_u32(sum0, sum1); - int32_t result = vaddvq_u32(total_sum); - - return static_cast(result); + // Widening horizontal sum: the lane total reaches 255*255*dim, so reducing into 32 bits + // would wrap from dimension 33,027. Callers rely on this being the exact integer dot. + return vaddlvq_u32(total_sum); } template // 0..63 float UINT8_InnerProductSIMD16_NEON_DOTPROD(const void *pVect1v, const void *pVect2v, size_t dimension) { - return 1.0f - UINT8_InnerProductImp(pVect1v, pVect2v, dimension); + return 1.0f - static_cast(UINT8_InnerProductImp(pVect1v, pVect2v, dimension)); } template // 0..63 float UINT8_CosineSIMD_NEON_DOTPROD(const void *pVect1v, const void *pVect2v, size_t dimension) { - float ip = UINT8_InnerProductImp(pVect1v, pVect2v, dimension); + float ip = static_cast(UINT8_InnerProductImp(pVect1v, pVect2v, dimension)); const float norm_v1 = load_unaligned(static_cast(pVect1v) + dimension); const float norm_v2 = load_unaligned(static_cast(pVect2v) + dimension); return 1.0f - ip / (norm_v1 * norm_v2); diff --git a/src/VecSim/spaces/IP/IP_NEON_SQ8_SQ8.h b/src/VecSim/spaces/IP/IP_NEON_SQ8_SQ8.h index 3e931ee0e..cbdac623b 100644 --- a/src/VecSim/spaces/IP/IP_NEON_SQ8_SQ8.h +++ b/src/VecSim/spaces/IP/IP_NEON_SQ8_SQ8.h @@ -39,30 +39,16 @@ using sq8 = vecsim_types::sq8; template // 0..63 float SQ8_SQ8_InnerProductSIMD64_NEON_IMP(const void *pVec1v, const void *pVec2v, size_t dimension) { - // Compute raw dot product using efficient UINT8 implementation - // UINT8_InnerProductImp processes 16 elements at a time using native uint8 instructions - float dot_product = UINT8_InnerProductImp(pVec1v, pVec2v, dimension); + // Exact integer dot product of the quantized bytes. + const uint64_t q_dot = UINT8_InnerProductImp(pVec1v, pVec2v, dimension); - // Get dequantization parameters and precomputed values from the end of pVec1 - // Layout: [data (dim)] [min (float)] [delta (float)] [sum (float)] + // Every term describes the reconstructions, derived from the stored integer sums. See sq8.h. const uint8_t *pVec1 = static_cast(pVec1v); const uint8_t *pVec2 = static_cast(pVec2v); - - const auto *params1 = pVec1 + dimension; - const float min1 = load_unaligned(params1 + sq8::MIN_VAL * sizeof(float)); - const float delta1 = load_unaligned(params1 + sq8::DELTA * sizeof(float)); - const float sum1 = load_unaligned(params1 + sq8::SUM * sizeof(float)); - - // Get dequantization parameters and precomputed values from the end of pVec2 - const auto *params2 = pVec2 + dimension; - const float min2 = load_unaligned(params2 + sq8::MIN_VAL * sizeof(float)); - const float delta2 = load_unaligned(params2 + sq8::DELTA * sizeof(float)); - const float sum2 = load_unaligned(params2 + sq8::SUM * sizeof(float)); - - // Apply algebraic formula using precomputed sums: - // IP = min1*sum2 + min2*sum1 + δ1*δ2 * Σ(q1*q2) - dim*min1*min2 - return min1 * sum2 + min2 * sum1 + delta1 * delta2 * dot_product - - static_cast(dimension) * min1 * min2; + return static_cast(sq8::reconstructed_ip( + sq8::min_val_of(pVec1, dimension), sq8::delta_of(pVec1, dimension), + sq8::q_sum_of(pVec1, dimension), sq8::min_val_of(pVec2, dimension), + sq8::delta_of(pVec2, dimension), sq8::q_sum_of(pVec2, dimension), q_dot, dimension)); } // SQ8-to-SQ8 Inner Product distance function diff --git a/src/VecSim/spaces/IP/IP_NEON_UINT8.h b/src/VecSim/spaces/IP/IP_NEON_UINT8.h index 2d2b3f555..cb302a04e 100644 --- a/src/VecSim/spaces/IP/IP_NEON_UINT8.h +++ b/src/VecSim/spaces/IP/IP_NEON_UINT8.h @@ -35,7 +35,7 @@ InnerProductStep(uint8_t *&pVect1, uint8_t *&pVect2, uint32x4_t &sum) { } template // 0..63 -float UINT8_InnerProductImp(const void *pVect1v, const void *pVect2v, size_t dimension) { +uint64_t UINT8_InnerProductImp(const void *pVect1v, const void *pVect2v, size_t dimension) { uint8_t *pVect1 = (uint8_t *)pVect1v; uint8_t *pVect2 = (uint8_t *)pVect2v; @@ -105,20 +105,19 @@ float UINT8_InnerProductImp(const void *pVect1v, const void *pVect2v, size_t dim uint32x4_t total_sum = vaddq_u32(sum0, sum1); - // Horizontal sum of the 4 elements in the combined sum register - int32_t result = vaddvq_u32(total_sum); - - return static_cast(result); + // Widening horizontal sum: the lane total reaches 255*255*dim, so reducing into 32 bits + // would wrap from dimension 33,027. Callers rely on this being the exact integer dot. + return vaddlvq_u32(total_sum); } template // 0..15 float UINT8_InnerProductSIMD16_NEON(const void *pVect1v, const void *pVect2v, size_t dimension) { - return 1.0f - UINT8_InnerProductImp(pVect1v, pVect2v, dimension); + return 1.0f - static_cast(UINT8_InnerProductImp(pVect1v, pVect2v, dimension)); } template // 0..63 float UINT8_CosineSIMD_NEON(const void *pVect1v, const void *pVect2v, size_t dimension) { - float ip = UINT8_InnerProductImp(pVect1v, pVect2v, dimension); + float ip = static_cast(UINT8_InnerProductImp(pVect1v, pVect2v, dimension)); const float norm_v1 = load_unaligned(static_cast(pVect1v) + dimension); const float norm_v2 = load_unaligned(static_cast(pVect2v) + dimension); return 1.0f - ip / (norm_v1 * norm_v2); diff --git a/src/VecSim/spaces/IP/IP_SVE_SQ8_SQ8.h b/src/VecSim/spaces/IP/IP_SVE_SQ8_SQ8.h index 93ddb76cb..a61dc1c08 100644 --- a/src/VecSim/spaces/IP/IP_SVE_SQ8_SQ8.h +++ b/src/VecSim/spaces/IP/IP_SVE_SQ8_SQ8.h @@ -38,30 +38,17 @@ using sq8 = vecsim_types::sq8; // Uses UINT8_InnerProductImp for efficient dot product computation with SVE template float SQ8_SQ8_InnerProductSIMD_SVE_IMP(const void *pVec1v, const void *pVec2v, size_t dimension) { - // Compute raw dot product using efficient UINT8 SVE implementation - // UINT8_InnerProductImp uses svdot_u32 for native uint8 dot product - float dot_product = + // Exact integer dot product of the quantized bytes. + const uint64_t q_dot = UINT8_InnerProductImp(pVec1v, pVec2v, dimension); - // Get dequantization parameters and precomputed values from the end of vectors - // Layout: [data (dim)] [min (float)] [delta (float)] [sum (float)] + // Every term describes the reconstructions, derived from the stored integer sums. See sq8.h. const uint8_t *pVec1 = static_cast(pVec1v); const uint8_t *pVec2 = static_cast(pVec2v); - - const auto *params1 = pVec1 + dimension; - const float min1 = load_unaligned(params1 + sq8::MIN_VAL * sizeof(float)); - const float delta1 = load_unaligned(params1 + sq8::DELTA * sizeof(float)); - const float sum1 = load_unaligned(params1 + sq8::SUM * sizeof(float)); - - const auto *params2 = pVec2 + dimension; - const float min2 = load_unaligned(params2 + sq8::MIN_VAL * sizeof(float)); - const float delta2 = load_unaligned(params2 + sq8::DELTA * sizeof(float)); - const float sum2 = load_unaligned(params2 + sq8::SUM * sizeof(float)); - - // Apply algebraic formula with float conversion only at the end: - // IP = min1*sum2 + min2*sum1 + δ1*δ2 * Σ(q1*q2) - dim*min1*min2 - return min1 * sum2 + min2 * sum1 + delta1 * delta2 * dot_product - - static_cast(dimension) * min1 * min2; + return static_cast(sq8::reconstructed_ip( + sq8::min_val_of(pVec1, dimension), sq8::delta_of(pVec1, dimension), + sq8::q_sum_of(pVec1, dimension), sq8::min_val_of(pVec2, dimension), + sq8::delta_of(pVec2, dimension), sq8::q_sum_of(pVec2, dimension), q_dot, dimension)); } // SQ8-to-SQ8 Inner Product distance function diff --git a/src/VecSim/spaces/IP/IP_SVE_UINT8.h b/src/VecSim/spaces/IP/IP_SVE_UINT8.h index f6c6af3b8..9c5c0cb73 100644 --- a/src/VecSim/spaces/IP/IP_SVE_UINT8.h +++ b/src/VecSim/spaces/IP/IP_SVE_UINT8.h @@ -24,7 +24,7 @@ inline void InnerProductStep(const uint8_t *&pVect1, const uint8_t *&pVect2, siz } template -float UINT8_InnerProductImp(const void *pVect1v, const void *pVect2v, size_t dimension) { +uint64_t UINT8_InnerProductImp(const void *pVect1v, const void *pVect2v, size_t dimension) { const uint8_t *pVect1 = reinterpret_cast(pVect1v); const uint8_t *pVect2 = reinterpret_cast(pVect2v); @@ -83,20 +83,22 @@ float UINT8_InnerProductImp(const void *pVect1v, const void *pVect2v, size_t dim sum2 = svadd_u32_x(svptrue_b32(), sum2, sum3); // Perform vector addition in parallel and Horizontal sum - int32_t sum_all = svaddv_u32(svptrue_b32(), svadd_u32_x(svptrue_b32(), sum0, sum2)); - - return sum_all; + // svaddv_u32 already reduces into a 64-bit scalar; the previous int32_t truncated it, which + // wrapped from dimension 33,027 since the total reaches 255*255*dim. Callers rely on this + // being the exact integer dot product. + return svaddv_u32(svptrue_b32(), svadd_u32_x(svptrue_b32(), sum0, sum2)); } template float UINT8_InnerProductSIMD_SVE(const void *pVect1v, const void *pVect2v, size_t dimension) { - return 1.0f - - UINT8_InnerProductImp(pVect1v, pVect2v, dimension); + return 1.0f - static_cast(UINT8_InnerProductImp( + pVect1v, pVect2v, dimension)); } template float UINT8_CosineSIMD_SVE(const void *pVect1v, const void *pVect2v, size_t dimension) { - float ip = UINT8_InnerProductImp(pVect1v, pVect2v, dimension); + float ip = static_cast( + UINT8_InnerProductImp(pVect1v, pVect2v, dimension)); const float norm_v1 = load_unaligned(static_cast(pVect1v) + dimension); const float norm_v2 = load_unaligned(static_cast(pVect2v) + dimension); return 1.0f - ip / (norm_v1 * norm_v2); diff --git a/src/VecSim/spaces/L2/L2.cpp b/src/VecSim/spaces/L2/L2.cpp index 015f2200d..37f371174 100644 --- a/src/VecSim/spaces/L2/L2.cpp +++ b/src/VecSim/spaces/L2/L2.cpp @@ -31,17 +31,21 @@ float SQ8_FP32_L2Sqr(const void *pVect1v, const void *pVect2v, size_t dimension) // Get the raw inner product using the common implementation const float ip = SQ8_FP32_InnerProduct_Impl(pVect1v, pVect2v, dimension); - // Storage metadata follows a byte payload and is not necessarily float-aligned. + // ||x||² describes the reconstruction, derived exactly from the stored integer sums. const auto *pVect1 = static_cast(pVect1v); - const float x_sum_sq = - load_unaligned(pVect1 + dimension + sq8::SUM_SQUARES * sizeof(float)); + const double x_sum_sq = sq8::reconstructed_sum_squares( + sq8::min_val_of(pVect1, dimension), sq8::delta_of(pVect1, dimension), + sq8::q_sum_of(pVect1, dimension), sq8::q_sum_squares_of(pVect1, dimension), dimension); // Get precomputed sum of squares from query blob (pVect2 is FP32) const auto *pVect2 = static_cast(pVect2v); const float y_sum_sq = pVect2[dimension + sq8::SUM_SQUARES_QUERY]; - // L2² = ||x||² + ||y||² - 2*IP(x, y) - return x_sum_sq + y_sum_sq - 2.0f * ip; + // L2² = ||x||² + ||y||² - 2*IP(x, y). Note ip is accumulated in FP32 by the shared kernel, so + // this remains vulnerable to cancellation when ||x||² greatly exceeds the true distance; the + // asymmetric path needs the difference taken inside the loop to fix that. + return static_cast(x_sum_sq + static_cast(y_sum_sq) - + 2.0 * static_cast(ip)); } /* @@ -62,15 +66,18 @@ float SQ8_FP16_L2Sqr(const void *pVect1v, const void *pVect2v, size_t dimension) // byte offsets that are not guaranteed 4-byte aligned for odd `dimension`, so use // load_unaligned to avoid alignment UB. const auto *pVect1 = static_cast(pVect1v); - const float x_sum_sq = - load_unaligned(pVect1 + dimension + sq8::SUM_SQUARES * sizeof(float)); + const double x_sum_sq = sq8::reconstructed_sum_squares( + sq8::min_val_of(pVect1, dimension), sq8::delta_of(pVect1, dimension), + sq8::q_sum_of(pVect1, dimension), sq8::q_sum_squares_of(pVect1, dimension), dimension); const auto *pVect2 = static_cast(pVect2v); const auto *query_meta_bytes = reinterpret_cast(pVect2 + dimension); const float y_sum_sq = load_unaligned(query_meta_bytes + sq8::SUM_SQUARES_QUERY * sizeof(float)); - // L2² = ||x||² + ||y||² - 2*IP(x, y) - return x_sum_sq + y_sum_sq - 2.0f * ip; + // L2² = ||x||² + ||y||² - 2*IP(x, y). Still cancellation-prone for the same reason as the FP32 + // variant: ip is accumulated in FP32 by the shared kernel. + return static_cast(x_sum_sq + static_cast(y_sum_sq) - + 2.0 * static_cast(ip)); } float FP32_L2Sqr(const void *pVect1v, const void *pVect2v, size_t dimension) { @@ -137,7 +144,7 @@ float FP16_L2Sqr(const void *pVect1, const void *pVect2, size_t dimension) { // To support dimension up to 2^16, we need the difference between the type and int_elem_t to be at // least 2 bytes. We assert that in the implementation. template -using ret_t = std::conditional_t; +using ret_t = long long; // Difference type for the L2 functions. // The type should be able to hold `MIN_VAL(int_elem_t)-MAX_VAL(int_elem_t)`, and should be signed @@ -186,16 +193,18 @@ float SQ8_SQ8_L2Sqr(const void *pVect1v, const void *pVect2v, size_t dimension) const auto *pVect1 = static_cast(pVect1v); const auto *pVect2 = static_cast(pVect2v); - // Get precomputed sum of squares from both vectors - // Layout: [uint8_t values (dim)] [min_val] [delta] [sum] [sum_of_squares] - const float sum_sq_1 = - load_unaligned(pVect1 + dimension + sq8::SUM_SQUARES * sizeof(float)); - const float sum_sq_2 = - load_unaligned(pVect2 + dimension + sq8::SUM_SQUARES * sizeof(float)); - - // Use the common inner product implementation - const float ip = SQ8_SQ8_InnerProduct_Impl(pVect1v, pVect2v, dimension); + // Integer dot product of the quantized bytes, exact. + uint64_t q_dot = 0; + for (size_t i = 0; i < dimension; i++) { + q_dot += static_cast(pVect1[i]) * static_cast(pVect2[i]); + } - // L2² = ||x||² + ||y||² - 2*IP(x, y) - return sum_sq_1 + sum_sq_2 - 2.0f * ip; + // Expanded difference rather than ||x||² + ||y||² - 2*IP: see sq8::reconstructed_l2_sqr for + // why the two are not interchangeable in floating point. + return static_cast(sq8::reconstructed_l2_sqr( + sq8::min_val_of(pVect1, dimension), sq8::delta_of(pVect1, dimension), + sq8::q_sum_of(pVect1, dimension), sq8::q_sum_squares_of(pVect1, dimension), + sq8::min_val_of(pVect2, dimension), sq8::delta_of(pVect2, dimension), + sq8::q_sum_of(pVect2, dimension), sq8::q_sum_squares_of(pVect2, dimension), q_dot, + dimension)); } diff --git a/src/VecSim/spaces/L2/L2_AVX2_FMA_SQ8_FP16.h b/src/VecSim/spaces/L2/L2_AVX2_FMA_SQ8_FP16.h index c855b62ca..b03d8d8d4 100644 --- a/src/VecSim/spaces/L2/L2_AVX2_FMA_SQ8_FP16.h +++ b/src/VecSim/spaces/L2/L2_AVX2_FMA_SQ8_FP16.h @@ -21,12 +21,13 @@ float SQ8_FP16_L2SqrSIMD16_AVX2_FMA(const void *pVect1v, const void *pVect2v, si const uint8_t *pVect1 = static_cast(pVect1v); const uint8_t *params_bytes = pVect1 + dimension; - const float x_sum_sq = load_unaligned(params_bytes + sq8::SUM_SQUARES * sizeof(float)); + const double x_sum_sq = sq8::reconstructed_sum_squares_at(params_bytes, dimension); const float16 *pVect2 = static_cast(pVect2v); const auto *query_meta_bytes = reinterpret_cast(pVect2 + dimension); const float y_sum_sq = load_unaligned(query_meta_bytes + sq8::SUM_SQUARES_QUERY * sizeof(float)); - return x_sum_sq + y_sum_sq - 2.0f * ip; + return static_cast(x_sum_sq + static_cast(y_sum_sq) - + 2.0 * static_cast(ip)); } diff --git a/src/VecSim/spaces/L2/L2_AVX2_FMA_SQ8_FP32.h b/src/VecSim/spaces/L2/L2_AVX2_FMA_SQ8_FP32.h index 0f9d5bde9..d98b2c7a3 100644 --- a/src/VecSim/spaces/L2/L2_AVX2_FMA_SQ8_FP32.h +++ b/src/VecSim/spaces/L2/L2_AVX2_FMA_SQ8_FP32.h @@ -35,12 +35,12 @@ float SQ8_FP32_L2SqrSIMD16_AVX2_FMA(const void *pVect1v, const void *pVect2v, si // Get precomputed sum of squares from storage blob (pVect1v is SQ8 storage) const uint8_t *pVect1 = static_cast(pVect1v); - const float x_sum_sq = - load_unaligned(pVect1 + dimension + sq8::SUM_SQUARES * sizeof(float)); + const double x_sum_sq = sq8::reconstructed_sum_squares_of(pVect1, dimension); // Get precomputed sum of squares from query blob (pVect2v is FP32 query) const float y_sum_sq = static_cast(pVect2v)[dimension + sq8::SUM_SQUARES_QUERY]; // L2² = ||x||² + ||y||² - 2*IP(x, y) - return x_sum_sq + y_sum_sq - 2.0f * ip; + return static_cast(x_sum_sq + static_cast(y_sum_sq) - + 2.0 * static_cast(ip)); } diff --git a/src/VecSim/spaces/L2/L2_AVX2_SQ8_FP16.h b/src/VecSim/spaces/L2/L2_AVX2_SQ8_FP16.h index 7c2cbfcd8..7d5d35b8d 100644 --- a/src/VecSim/spaces/L2/L2_AVX2_SQ8_FP16.h +++ b/src/VecSim/spaces/L2/L2_AVX2_SQ8_FP16.h @@ -21,12 +21,13 @@ float SQ8_FP16_L2SqrSIMD16_AVX2(const void *pVect1v, const void *pVect2v, size_t const uint8_t *pVect1 = static_cast(pVect1v); const uint8_t *params_bytes = pVect1 + dimension; - const float x_sum_sq = load_unaligned(params_bytes + sq8::SUM_SQUARES * sizeof(float)); + const double x_sum_sq = sq8::reconstructed_sum_squares_at(params_bytes, dimension); const float16 *pVect2 = static_cast(pVect2v); const auto *query_meta_bytes = reinterpret_cast(pVect2 + dimension); const float y_sum_sq = load_unaligned(query_meta_bytes + sq8::SUM_SQUARES_QUERY * sizeof(float)); - return x_sum_sq + y_sum_sq - 2.0f * ip; + return static_cast(x_sum_sq + static_cast(y_sum_sq) - + 2.0 * static_cast(ip)); } diff --git a/src/VecSim/spaces/L2/L2_AVX2_SQ8_FP32.h b/src/VecSim/spaces/L2/L2_AVX2_SQ8_FP32.h index d08474f71..491ed2e98 100644 --- a/src/VecSim/spaces/L2/L2_AVX2_SQ8_FP32.h +++ b/src/VecSim/spaces/L2/L2_AVX2_SQ8_FP32.h @@ -35,12 +35,12 @@ float SQ8_FP32_L2SqrSIMD16_AVX2(const void *pVect1v, const void *pVect2v, size_t // Get precomputed sum of squares from storage blob (pVect1v is SQ8 storage) const uint8_t *pVect1 = static_cast(pVect1v); - const float x_sum_sq = - load_unaligned(pVect1 + dimension + sq8::SUM_SQUARES * sizeof(float)); + const double x_sum_sq = sq8::reconstructed_sum_squares_of(pVect1, dimension); // Get precomputed sum of squares from query blob (pVect2v is FP32 query) const float y_sum_sq = static_cast(pVect2v)[dimension + sq8::SUM_SQUARES_QUERY]; // L2² = ||x||² + ||y||² - 2*IP(x, y) - return x_sum_sq + y_sum_sq - 2.0f * ip; + return static_cast(x_sum_sq + static_cast(y_sum_sq) - + 2.0 * static_cast(ip)); } diff --git a/src/VecSim/spaces/L2/L2_AVX512F_BW_VL_VNNI_SQ8_FP32.h b/src/VecSim/spaces/L2/L2_AVX512F_BW_VL_VNNI_SQ8_FP32.h index f1c86d689..61cd7e45b 100644 --- a/src/VecSim/spaces/L2/L2_AVX512F_BW_VL_VNNI_SQ8_FP32.h +++ b/src/VecSim/spaces/L2/L2_AVX512F_BW_VL_VNNI_SQ8_FP32.h @@ -35,12 +35,12 @@ float SQ8_FP32_L2SqrSIMD16_AVX512F_BW_VL_VNNI(const void *pVect1v, const void *p // Get precomputed sum of squares from storage blob (pVect1v is SQ8 storage) const uint8_t *pVect1 = static_cast(pVect1v); - const float x_sum_sq = - load_unaligned(pVect1 + dimension + sq8::SUM_SQUARES * sizeof(float)); + const double x_sum_sq = sq8::reconstructed_sum_squares_of(pVect1, dimension); // Get precomputed sum of squares from query blob (pVect2v is FP32 query) const float y_sum_sq = static_cast(pVect2v)[dimension + sq8::SUM_SQUARES_QUERY]; // L2² = ||x||² + ||y||² - 2*IP(x, y) - return x_sum_sq + y_sum_sq - 2.0f * ip; + return static_cast(x_sum_sq + static_cast(y_sum_sq) - + 2.0 * static_cast(ip)); } diff --git a/src/VecSim/spaces/L2/L2_AVX512F_BW_VL_VNNI_SQ8_SQ8.h b/src/VecSim/spaces/L2/L2_AVX512F_BW_VL_VNNI_SQ8_SQ8.h index 6a862b2ed..7a94d59d5 100644 --- a/src/VecSim/spaces/L2/L2_AVX512F_BW_VL_VNNI_SQ8_SQ8.h +++ b/src/VecSim/spaces/L2/L2_AVX512F_BW_VL_VNNI_SQ8_SQ8.h @@ -8,6 +8,7 @@ */ #pragma once #include "VecSim/spaces/space_includes.h" +#include "VecSim/types/sq8.h" #include "VecSim/spaces/IP/IP_AVX512F_BW_VL_VNNI_SQ8_SQ8.h" /** @@ -26,19 +27,19 @@ template // 0..63 float SQ8_SQ8_L2SqrSIMD64_AVX512F_BW_VL_VNNI(const void *pVec1v, const void *pVec2v, size_t dimension) { - - // Use the common inner product implementation (returns raw IP, not distance) - const float ip = SQ8_SQ8_InnerProductImp(pVec1v, pVec2v, dimension); + // Exact integer dot product of the quantized bytes. Taken from the shared UINT8 helper rather + // than from the SQ8 inner-product wrapper, whose FP32 return would reintroduce the rounding + // this formulation exists to avoid. + const uint64_t q_dot = UINT8_InnerProductImp(pVec1v, pVec2v, dimension); const uint8_t *pVec1 = static_cast(pVec1v); const uint8_t *pVec2 = static_cast(pVec2v); - // Get precomputed sum of squares from both vectors - // Layout: [uint8_t values (dim)] [min_val] [delta] [sum] [sum_of_squares] - const float sum_sq_1 = - load_unaligned(pVec1 + dimension + sq8::SUM_SQUARES * sizeof(float)); - const float sum_sq_2 = - load_unaligned(pVec2 + dimension + sq8::SUM_SQUARES * sizeof(float)); - - // L2² = ||x||² + ||y||² - 2*IP(x, y) - return sum_sq_1 + sum_sq_2 - 2.0f * ip; + // Expanded difference rather than ||x||^2 + ||y||^2 - 2*IP: algebraically identical, but the + // expansion keeps any common offset in (min1 - min2) instead of letting it cancel the answer. + return static_cast(sq8::reconstructed_l2_sqr( + sq8::min_val_of(pVec1, dimension), sq8::delta_of(pVec1, dimension), + sq8::q_sum_of(pVec1, dimension), sq8::q_sum_squares_of(pVec1, dimension), + sq8::min_val_of(pVec2, dimension), sq8::delta_of(pVec2, dimension), + sq8::q_sum_of(pVec2, dimension), sq8::q_sum_squares_of(pVec2, dimension), q_dot, + dimension)); } diff --git a/src/VecSim/spaces/L2/L2_AVX512F_BW_VL_VNNI_UINT8.h b/src/VecSim/spaces/L2/L2_AVX512F_BW_VL_VNNI_UINT8.h index 350b759ea..f906e7dc4 100644 --- a/src/VecSim/spaces/L2/L2_AVX512F_BW_VL_VNNI_UINT8.h +++ b/src/VecSim/spaces/L2/L2_AVX512F_BW_VL_VNNI_UINT8.h @@ -92,5 +92,7 @@ float UINT8_L2SqrSIMD64_AVX512F_BW_VL_VNNI(const void *pVect1v, const void *pVec } while (pVect1 < pEnd1); } - return _mm512_reduce_add_epi32(sum); + // The lanes hold sums of squared byte differences, so the horizontal total is unsigned and + // reaches 255*255*dim. Reading it as a signed int wrapped it negative from dimension 33,026. + return static_cast(static_cast(_mm512_reduce_add_epi32(sum))); } diff --git a/src/VecSim/spaces/L2/L2_AVX512F_SQ8_FP16.h b/src/VecSim/spaces/L2/L2_AVX512F_SQ8_FP16.h index 9d7b1569f..0f4ed27f9 100644 --- a/src/VecSim/spaces/L2/L2_AVX512F_SQ8_FP16.h +++ b/src/VecSim/spaces/L2/L2_AVX512F_SQ8_FP16.h @@ -21,12 +21,13 @@ float SQ8_FP16_L2SqrSIMD16_AVX512F(const void *pVect1v, const void *pVect2v, siz const uint8_t *pVect1 = static_cast(pVect1v); const uint8_t *params_bytes = pVect1 + dimension; - const float x_sum_sq = load_unaligned(params_bytes + sq8::SUM_SQUARES * sizeof(float)); + const double x_sum_sq = sq8::reconstructed_sum_squares_at(params_bytes, dimension); const float16 *pVect2 = static_cast(pVect2v); const auto *query_meta_bytes = reinterpret_cast(pVect2 + dimension); const float y_sum_sq = load_unaligned(query_meta_bytes + sq8::SUM_SQUARES_QUERY * sizeof(float)); - return x_sum_sq + y_sum_sq - 2.0f * ip; + return static_cast(x_sum_sq + static_cast(y_sum_sq) - + 2.0 * static_cast(ip)); } diff --git a/src/VecSim/spaces/L2/L2_NEON_DOTPROD_SQ8_SQ8.h b/src/VecSim/spaces/L2/L2_NEON_DOTPROD_SQ8_SQ8.h index 6c55b3f03..b9d8a19dd 100644 --- a/src/VecSim/spaces/L2/L2_NEON_DOTPROD_SQ8_SQ8.h +++ b/src/VecSim/spaces/L2/L2_NEON_DOTPROD_SQ8_SQ8.h @@ -28,20 +28,19 @@ using sq8 = vecsim_types::sq8; // L2 squared distance using the common inner product implementation template // 0..63 float SQ8_SQ8_L2SqrSIMD64_NEON_DOTPROD(const void *pVec1v, const void *pVec2v, size_t dimension) { - // Use the common inner product implementation (returns raw IP, not distance) - const float ip = - SQ8_SQ8_InnerProductSIMD64_NEON_DOTPROD_IMP(pVec1v, pVec2v, dimension); + // Exact integer dot product of the quantized bytes. Taken from the shared UINT8 helper rather + // than from the SQ8 inner-product wrapper, whose FP32 return would reintroduce the rounding + // this formulation exists to avoid. + const uint64_t q_dot = UINT8_InnerProductImp(pVec1v, pVec2v, dimension); const uint8_t *pVec1 = static_cast(pVec1v); const uint8_t *pVec2 = static_cast(pVec2v); - - // Get precomputed sum of squares from both vectors - // Layout: [uint8_t values (dim)] [min_val] [delta] [sum] [sum_of_squares] - const float sum_sq_1 = - load_unaligned(pVec1 + dimension + sq8::SUM_SQUARES * sizeof(float)); - const float sum_sq_2 = - load_unaligned(pVec2 + dimension + sq8::SUM_SQUARES * sizeof(float)); - - // L2² = ||x||² + ||y||² - 2*IP(x, y) - return sum_sq_1 + sum_sq_2 - 2.0f * ip; + // Expanded difference rather than ||x||^2 + ||y||^2 - 2*IP: algebraically identical, but the + // expansion keeps any common offset in (min1 - min2) instead of letting it cancel the answer. + return static_cast(sq8::reconstructed_l2_sqr( + sq8::min_val_of(pVec1, dimension), sq8::delta_of(pVec1, dimension), + sq8::q_sum_of(pVec1, dimension), sq8::q_sum_squares_of(pVec1, dimension), + sq8::min_val_of(pVec2, dimension), sq8::delta_of(pVec2, dimension), + sq8::q_sum_of(pVec2, dimension), sq8::q_sum_squares_of(pVec2, dimension), q_dot, + dimension)); } diff --git a/src/VecSim/spaces/L2/L2_NEON_SQ8_FP16.h b/src/VecSim/spaces/L2/L2_NEON_SQ8_FP16.h index 2964c1cee..a5922361c 100644 --- a/src/VecSim/spaces/L2/L2_NEON_SQ8_FP16.h +++ b/src/VecSim/spaces/L2/L2_NEON_SQ8_FP16.h @@ -24,14 +24,15 @@ float SQ8_FP16_L2SqrSIMD16_NEON_HP(const void *pVect1v, const void *pVect2v, siz const float ip = SQ8_FP16_InnerProductSIMD16_NEON_HP_IMP(pVect1v, pVect2v, dimension); const uint8_t *params_bytes = static_cast(pVect1v) + dimension; - const float x_sum_sq = load_unaligned(params_bytes + sq8::SUM_SQUARES * sizeof(float)); + const double x_sum_sq = sq8::reconstructed_sum_squares_at(params_bytes, dimension); const uint8_t *query_meta_bytes = reinterpret_cast(static_cast(pVect2v) + dimension); const float y_sum_sq = load_unaligned(query_meta_bytes + sq8::SUM_SQUARES_QUERY * sizeof(float)); - return x_sum_sq + y_sum_sq - 2.0f * ip; + return static_cast(x_sum_sq + static_cast(y_sum_sq) - + 2.0 * static_cast(ip)); } // FMLAL (FEAT_FHM) variant — same identity, FMLAL widening-FMA IP core. @@ -41,12 +42,13 @@ float SQ8_FP16_L2SqrSIMD16_NEON_FHM(const void *pVect1v, const void *pVect2v, si SQ8_FP16_InnerProductSIMD16_NEON_HP_IMP(pVect1v, pVect2v, dimension); const uint8_t *params_bytes = static_cast(pVect1v) + dimension; - const float x_sum_sq = load_unaligned(params_bytes + sq8::SUM_SQUARES * sizeof(float)); + const double x_sum_sq = sq8::reconstructed_sum_squares_at(params_bytes, dimension); const uint8_t *query_meta_bytes = reinterpret_cast(static_cast(pVect2v) + dimension); const float y_sum_sq = load_unaligned(query_meta_bytes + sq8::SUM_SQUARES_QUERY * sizeof(float)); - return x_sum_sq + y_sum_sq - 2.0f * ip; + return static_cast(x_sum_sq + static_cast(y_sum_sq) - + 2.0 * static_cast(ip)); } diff --git a/src/VecSim/spaces/L2/L2_NEON_SQ8_FP32.h b/src/VecSim/spaces/L2/L2_NEON_SQ8_FP32.h index f6f7a6bc0..273bf4d74 100644 --- a/src/VecSim/spaces/L2/L2_NEON_SQ8_FP32.h +++ b/src/VecSim/spaces/L2/L2_NEON_SQ8_FP32.h @@ -36,12 +36,12 @@ float SQ8_FP32_L2SqrSIMD16_NEON(const void *pVect1v, const void *pVect2v, size_t // Get precomputed sum of squares from storage blob (pVect1v is SQ8 storage) const uint8_t *pVect1 = static_cast(pVect1v); - const float x_sum_sq = - load_unaligned(pVect1 + dimension + sq8::SUM_SQUARES * sizeof(float)); + const double x_sum_sq = sq8::reconstructed_sum_squares_of(pVect1, dimension); // Get precomputed sum of squares from query blob (pVect2v is FP32 query) const float y_sum_sq = static_cast(pVect2v)[dimension + sq8::SUM_SQUARES_QUERY]; // L2² = ||x||² + ||y||² - 2*IP(x, y) - return x_sum_sq + y_sum_sq - 2.0f * ip; + return static_cast(x_sum_sq + static_cast(y_sum_sq) - + 2.0 * static_cast(ip)); } diff --git a/src/VecSim/spaces/L2/L2_NEON_SQ8_SQ8.h b/src/VecSim/spaces/L2/L2_NEON_SQ8_SQ8.h index 2750997f7..bcb651db3 100644 --- a/src/VecSim/spaces/L2/L2_NEON_SQ8_SQ8.h +++ b/src/VecSim/spaces/L2/L2_NEON_SQ8_SQ8.h @@ -28,19 +28,19 @@ using sq8 = vecsim_types::sq8; // L2 squared distance using the common inner product implementation template // 0..63 float SQ8_SQ8_L2SqrSIMD64_NEON(const void *pVec1v, const void *pVec2v, size_t dimension) { - // Use the common inner product implementation (returns raw IP, not distance) - const float ip = SQ8_SQ8_InnerProductSIMD64_NEON_IMP(pVec1v, pVec2v, dimension); + // Exact integer dot product of the quantized bytes. Taken from the shared UINT8 helper rather + // than from the SQ8 inner-product wrapper, whose FP32 return would reintroduce the rounding + // this formulation exists to avoid. + const uint64_t q_dot = UINT8_InnerProductImp(pVec1v, pVec2v, dimension); const uint8_t *pVec1 = static_cast(pVec1v); const uint8_t *pVec2 = static_cast(pVec2v); - - // Get precomputed sum of squares from both vectors - // Layout: [uint8_t values (dim)] [min_val] [delta] [sum] [sum_of_squares] - const float sum_sq_1 = - load_unaligned(pVec1 + dimension + sq8::SUM_SQUARES * sizeof(float)); - const float sum_sq_2 = - load_unaligned(pVec2 + dimension + sq8::SUM_SQUARES * sizeof(float)); - - // L2² = ||x||² + ||y||² - 2*IP(x, y) - return sum_sq_1 + sum_sq_2 - 2.0f * ip; + // Expanded difference rather than ||x||^2 + ||y||^2 - 2*IP: algebraically identical, but the + // expansion keeps any common offset in (min1 - min2) instead of letting it cancel the answer. + return static_cast(sq8::reconstructed_l2_sqr( + sq8::min_val_of(pVec1, dimension), sq8::delta_of(pVec1, dimension), + sq8::q_sum_of(pVec1, dimension), sq8::q_sum_squares_of(pVec1, dimension), + sq8::min_val_of(pVec2, dimension), sq8::delta_of(pVec2, dimension), + sq8::q_sum_of(pVec2, dimension), sq8::q_sum_squares_of(pVec2, dimension), q_dot, + dimension)); } diff --git a/src/VecSim/spaces/L2/L2_NEON_UINT8.h b/src/VecSim/spaces/L2/L2_NEON_UINT8.h index aa3769867..bd280a750 100644 --- a/src/VecSim/spaces/L2/L2_NEON_UINT8.h +++ b/src/VecSim/spaces/L2/L2_NEON_UINT8.h @@ -126,7 +126,9 @@ float UINT8_L2SqrSIMD16_NEON(const void *pVect1v, const void *pVect2v, size_t di total_sum = vaddq_u32(total_sum, sum3); // Horizontal sum of the 4 elements in the combined sum register - int32_t result = vaddvq_u32(total_sum); + // Unsigned: the total is a sum of squared byte differences, reaching 255*255*dim. As a signed + // int32 this wrapped negative from dimension 33,026. + uint32_t result = vaddvq_u32(total_sum); // Return the L2 squared distance as a float return static_cast(result); diff --git a/src/VecSim/spaces/L2/L2_SSE4_SQ8_FP16.h b/src/VecSim/spaces/L2/L2_SSE4_SQ8_FP16.h index d0a0fea06..3d80a39ef 100644 --- a/src/VecSim/spaces/L2/L2_SSE4_SQ8_FP16.h +++ b/src/VecSim/spaces/L2/L2_SSE4_SQ8_FP16.h @@ -20,12 +20,13 @@ float SQ8_FP16_L2SqrSIMD16_SSE4(const void *pVect1v, const void *pVect2v, size_t const uint8_t *pVect1 = static_cast(pVect1v); const uint8_t *params_bytes = pVect1 + dimension; - const float x_sum_sq = load_unaligned(params_bytes + sq8::SUM_SQUARES * sizeof(float)); + const double x_sum_sq = sq8::reconstructed_sum_squares_at(params_bytes, dimension); const float16 *pVect2 = static_cast(pVect2v); const auto *query_meta_bytes = reinterpret_cast(pVect2 + dimension); const float y_sum_sq = load_unaligned(query_meta_bytes + sq8::SUM_SQUARES_QUERY * sizeof(float)); - return x_sum_sq + y_sum_sq - 2.0f * ip; + return static_cast(x_sum_sq + static_cast(y_sum_sq) - + 2.0 * static_cast(ip)); } diff --git a/src/VecSim/spaces/L2/L2_SSE4_SQ8_FP32.h b/src/VecSim/spaces/L2/L2_SSE4_SQ8_FP32.h index 3a3a4d12d..943ecf2b0 100644 --- a/src/VecSim/spaces/L2/L2_SSE4_SQ8_FP32.h +++ b/src/VecSim/spaces/L2/L2_SSE4_SQ8_FP32.h @@ -35,12 +35,12 @@ float SQ8_FP32_L2SqrSIMD16_SSE4(const void *pVect1v, const void *pVect2v, size_t // Get precomputed sum of squares from storage blob (pVect1v is SQ8 storage) const uint8_t *pVect1 = static_cast(pVect1v); - const float x_sum_sq = - load_unaligned(pVect1 + dimension + sq8::SUM_SQUARES * sizeof(float)); + const double x_sum_sq = sq8::reconstructed_sum_squares_of(pVect1, dimension); // Get precomputed sum of squares from query blob (pVect2v is FP32 query) const float y_sum_sq = static_cast(pVect2v)[dimension + sq8::SUM_SQUARES_QUERY]; // L2² = ||x||² + ||y||² - 2*IP(x, y) - return x_sum_sq + y_sum_sq - 2.0f * ip; + return static_cast(x_sum_sq + static_cast(y_sum_sq) - + 2.0 * static_cast(ip)); } diff --git a/src/VecSim/spaces/L2/L2_SVE2_SQ8_FP16.h b/src/VecSim/spaces/L2/L2_SVE2_SQ8_FP16.h index d9451fe2a..3c0cdf343 100644 --- a/src/VecSim/spaces/L2/L2_SVE2_SQ8_FP16.h +++ b/src/VecSim/spaces/L2/L2_SVE2_SQ8_FP16.h @@ -22,11 +22,12 @@ float SQ8_FP16_L2SqrSIMD_SVE2(const void *pVect1v, const void *pVect2v, size_t d pVect1v, pVect2v, dimension); const uint8_t *params_bytes = static_cast(pVect1v) + dimension; - const float x_sum_sq = load_unaligned(params_bytes + sq8::SUM_SQUARES * sizeof(float)); + const double x_sum_sq = sq8::reconstructed_sum_squares_at(params_bytes, dimension); const uint8_t *query_meta_bytes = reinterpret_cast(static_cast(pVect2v) + dimension); const float y_sum_sq = load_unaligned(query_meta_bytes + sq8::SUM_SQUARES_QUERY * sizeof(float)); - return x_sum_sq + y_sum_sq - 2.0f * ip; + return static_cast(x_sum_sq + static_cast(y_sum_sq) - + 2.0 * static_cast(ip)); } diff --git a/src/VecSim/spaces/L2/L2_SVE_SQ8_FP16.h b/src/VecSim/spaces/L2/L2_SVE_SQ8_FP16.h index f70ef493d..4848b5ceb 100644 --- a/src/VecSim/spaces/L2/L2_SVE_SQ8_FP16.h +++ b/src/VecSim/spaces/L2/L2_SVE_SQ8_FP16.h @@ -22,11 +22,12 @@ float SQ8_FP16_L2SqrSIMD_SVE(const void *pVect1v, const void *pVect2v, size_t di pVect1v, pVect2v, dimension); const uint8_t *params_bytes = static_cast(pVect1v) + dimension; - const float x_sum_sq = load_unaligned(params_bytes + sq8::SUM_SQUARES * sizeof(float)); + const double x_sum_sq = sq8::reconstructed_sum_squares_at(params_bytes, dimension); const uint8_t *query_meta_bytes = reinterpret_cast(static_cast(pVect2v) + dimension); const float y_sum_sq = load_unaligned(query_meta_bytes + sq8::SUM_SQUARES_QUERY * sizeof(float)); - return x_sum_sq + y_sum_sq - 2.0f * ip; + return static_cast(x_sum_sq + static_cast(y_sum_sq) - + 2.0 * static_cast(ip)); } diff --git a/src/VecSim/spaces/L2/L2_SVE_SQ8_FP32.h b/src/VecSim/spaces/L2/L2_SVE_SQ8_FP32.h index 3f95c4ae4..cd1c979dc 100644 --- a/src/VecSim/spaces/L2/L2_SVE_SQ8_FP32.h +++ b/src/VecSim/spaces/L2/L2_SVE_SQ8_FP32.h @@ -37,12 +37,12 @@ float SQ8_FP32_L2SqrSIMD_SVE(const void *pVect1v, const void *pVect2v, size_t di // Get precomputed sum of squares from storage blob (pVect1v is SQ8 storage) const uint8_t *pVect1 = static_cast(pVect1v); - const float x_sum_sq = - load_unaligned(pVect1 + dimension + sq8::SUM_SQUARES * sizeof(float)); + const double x_sum_sq = sq8::reconstructed_sum_squares_of(pVect1, dimension); // Get precomputed sum of squares from query blob (pVect2v is FP32 query) const float y_sum_sq = static_cast(pVect2v)[dimension + sq8::SUM_SQUARES_QUERY]; // L2² = ||x||² + ||y||² - 2*IP(x, y) - return x_sum_sq + y_sum_sq - 2.0f * ip; + return static_cast(x_sum_sq + static_cast(y_sum_sq) - + 2.0 * static_cast(ip)); } diff --git a/src/VecSim/spaces/L2/L2_SVE_SQ8_SQ8.h b/src/VecSim/spaces/L2/L2_SVE_SQ8_SQ8.h index 5fe194d80..a8a906b03 100644 --- a/src/VecSim/spaces/L2/L2_SVE_SQ8_SQ8.h +++ b/src/VecSim/spaces/L2/L2_SVE_SQ8_SQ8.h @@ -28,20 +28,20 @@ using sq8 = vecsim_types::sq8; // L2 squared distance using the common inner product implementation template float SQ8_SQ8_L2SqrSIMD_SVE(const void *pVec1v, const void *pVec2v, size_t dimension) { - // Use the common inner product implementation (returns raw IP, not distance) - const float ip = SQ8_SQ8_InnerProductSIMD_SVE_IMP( - pVec1v, pVec2v, dimension); + // Exact integer dot product of the quantized bytes. Taken from the shared UINT8 helper rather + // than from the SQ8 inner-product wrapper, whose FP32 return would reintroduce the rounding + // this formulation exists to avoid. + const uint64_t q_dot = + UINT8_InnerProductImp(pVec1v, pVec2v, dimension); const uint8_t *pVec1 = static_cast(pVec1v); const uint8_t *pVec2 = static_cast(pVec2v); - - // Get precomputed sum of squares from both vectors - // Layout: [uint8_t values (dim)] [min_val] [delta] [sum] [sum_of_squares] - const float sum_sq_1 = - load_unaligned(pVec1 + dimension + sq8::SUM_SQUARES * sizeof(float)); - const float sum_sq_2 = - load_unaligned(pVec2 + dimension + sq8::SUM_SQUARES * sizeof(float)); - - // L2² = ||x||² + ||y||² - 2*IP(x, y) - return sum_sq_1 + sum_sq_2 - 2.0f * ip; + // Expanded difference rather than ||x||^2 + ||y||^2 - 2*IP: algebraically identical, but the + // expansion keeps any common offset in (min1 - min2) instead of letting it cancel the answer. + return static_cast(sq8::reconstructed_l2_sqr( + sq8::min_val_of(pVec1, dimension), sq8::delta_of(pVec1, dimension), + sq8::q_sum_of(pVec1, dimension), sq8::q_sum_squares_of(pVec1, dimension), + sq8::min_val_of(pVec2, dimension), sq8::delta_of(pVec2, dimension), + sq8::q_sum_of(pVec2, dimension), sq8::q_sum_squares_of(pVec2, dimension), q_dot, + dimension)); } diff --git a/src/VecSim/spaces/computer/preprocessors.h b/src/VecSim/spaces/computer/preprocessors.h index 195be1418..532eff487 100644 --- a/src/VecSim/spaces/computer/preprocessors.h +++ b/src/VecSim/spaces/computer/preprocessors.h @@ -15,6 +15,7 @@ #include #include #include +#include #include #include #include @@ -151,10 +152,15 @@ class CosinePreprocessor : public PreprocessorInterface { * and then scaling the values to fit in the range of [0, 255]. * * Storage layout: - * | quantized_values[dim] | min_val | delta | x_sum | (x_sum_squares for L2 only) | + * | quantized_values[dim] | min_val | delta | q_sum | (q_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. + * q_sum = Σa_i: sum of the quantized bytes (uint32), + * q_sum_squares = Σa_i²: sum of their squares (uint32). + * + * The sums are over the quantized bytes, not the original values, because a stored blob only ever + * describes its reconstruction min + delta * a_i. Anything derived from the blob must describe that + * same vector; sums over the input belong to a different one. They are integers so that they stay + * exact, which the L2 formulation depends on. See sq8.h for the derivations that consume them. * * Storage metadata is always FP32 (independent of DataType) to match the asymmetric distance * kernels. The quantized blob size is: @@ -195,11 +201,14 @@ class CosinePreprocessor : public PreprocessorInterface { * * For L2: * ||x - y||² = Σx_i² - 2*Σ(x_i * y_i) + Σy_i² - * = x_sum_squares - 2 * IP(x, y) + y_sum_squares * where: - * - x_sum_squares = Σx_i² is precomputed and stored in the storage blob + * - Σx_i² describes the reconstruction and is derived from q_sum and q_sum_squares by + * sq8::reconstructed_sum_squares, in double * - IP(x, y) is computed using the formula above - * - y_sum_squares = Σy_i² is precomputed and stored in the query blob + * - y_sum_squares = Σy_i² is precomputed and stored in the query blob, in FP32, because the + * query is never quantized + * Note this identity is only safe while Σx_i² is close to the answer. It is not when the two + * vectors share a large offset, which is why the symmetric case below does not use it. * For normalized L2, x and y in this formula are the centered values x' and y'; their distance * is identical to the distance between the original vectors. * @@ -213,16 +222,20 @@ class CosinePreprocessor : public PreprocessorInterface { * = dim * min_x * min_y * + min_x * (sum_y - dim * min_y) + min_y * (sum_x - dim * min_x) * + delta_x * delta_y * Σ(qx_i * qy_i) - * = min_x * sum_y + min_y * sum_x - dim * min_x * min_y - * + 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) - * - Σqy_i = (sum_y - dim * min_y) / delta_y + * where Σqx_i and Σqy_i are read directly from each blob's q_sum, and Σ(qx_i * qy_i) is the + * integer dot product of the quantized bytes. Every term is exact until the final combination, + * which happens in double (sq8::reconstructed_ip). * * For L2: - * ||x - y||² = sum_sq_x + sum_sq_y - 2 * IP(x, y) - * where sum_sq_x, sum_sq_y are precomputed sums of the squared represented values. + * ||x - y||² = Σ((min_x - min_y) + delta_x * qx_i - delta_y * qy_i)² + * = dim * c² + delta_x² * Σqx_i² + delta_y² * Σqy_i² + * + 2c * delta_x * Σqx_i - 2c * delta_y * Σqy_i + * - 2 * delta_x * delta_y * Σ(qx_i * qy_i), with c = min_x - min_y + * Deliberately not sum_sq_x + sum_sq_y - 2 * IP(x, y), which is algebraically identical but + * subtracts two large nearly-equal numbers: for vectors sharing a large offset the answer is + * smaller than one ULP of either term and is lost entirely. Expanding the difference first keeps + * the offset in c, formed once and exactly, so the summed terms scale with the spread instead. + * See sq8::reconstructed_l2_sqr. */ // Input types accepted by QuantPreprocessor. Opt-in via std::same_as so unrelated types // (e.g. integers, double, bfloat16) are rejected at the template head with a named constraint. @@ -261,85 +274,110 @@ class QuantPreprocessor : public PreprocessorInterface { static_assert(!WithNorm || Metric != VecSimMetric_Cosine, "WithNorm does not support Cosine metric."); - // Helper function to perform quantization. This function is used by the storage preprocessing - // methods. + // Quantizes one input vector into dim bytes followed by its metadata run. Used by the storage + // preprocessing methods. + // + // Scaling is in double: [-FLT_MAX, +FLT_MAX] is valid input and max - min overflows FP32. + // Three things together keep NaN out of the conversion to a byte, and none is sufficient + // alone. find_min_max clamps its endpoints, because FP32 centering can already hold an inf and + // because min is stored as FP32 and could not carry one anyway. delta is guarded rather than + // diff, since (float)(diff / 255) underflows to zero below diff ~1.8e-43 while diff itself is + // still nonzero, leaving 1/delta = inf and scaling the minimum element, whose numerator is + // exactly zero, to 0 * inf. And to_byte bounds with fmin/fmax, which return the non-NaN + // operand where std::clamp propagates it. The first two make NaN unreachable on their own, so + // the third is a backstop, kept because this function has produced three separate + // conversion-UB defects. delta = 1 is the degenerate handling min == max already got. + // + // The sums are exact uint32 over the quantized bytes rather than the input floats; see the + // class documentation above for why, and sq8.h for the derivations that consume them. void quantize(const DataType *input, OUTPUT_TYPE *quantized) const { assert(input && quantized); float x_mean_ip = 0.0f; // only used for WithNorm auto [min_val, max_val] = find_min_max(input, x_mean_ip); - // Calculate scaling factor (typed as MetadataType because they end up as metadata). - const MetadataType diff = (max_val - min_val); - const MetadataType delta = (diff == 0.0f) ? MetadataType{1} : diff / MetadataType{255}; - const MetadataType inv_delta = MetadataType{1} / delta; + const double diff = static_cast(max_val) - static_cast(min_val); + MetadataType delta = static_cast(diff / 255.0); + if (!std::isfinite(delta) || delta <= MetadataType{0}) + delta = MetadataType{1}; - // 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{}; + const double inv_delta = 1.0 / static_cast(delta); + const double min_val_d = static_cast(min_val); - // 4 independent accumulators (sum of squares), only used for L2 - float q0{}, q1{}, q2{}, q3{}; + const auto to_byte = [inv_delta, min_val_d](float x) { + const double scaled = (static_cast(x) - min_val_d) * inv_delta; + // Bounded first, so +0.5 before truncating is std::round without the library call: + // this TU is compiled at the x86-64 baseline, where round() is an out-of-line call. + return static_cast(std::fmin(std::fmax(scaled, 0.0), 255.0) + 0.5); + }; + + // 4 independent accumulators each, mirroring the unrolled loop below. q* is L2 only. + uint32_t s0{}, s1{}, s2{}, s3{}; + uint32_t q0{}, q1{}, q2{}, q3{}; size_t i = 0; - // round dim down to the nearest multiple of 4 size_t dim_round_down = this->dim & ~size_t(3); - // Quantize the values for (; i < dim_round_down; i += 4) { - // Load once (widened to FP32 if DataType is FP16). + // transformed_value widens to FP32 if DataType is FP16, and centers if WithNorm. const float x0 = transformed_value(input, i); const float x1 = transformed_value(input, i + 1); const float x2 = transformed_value(input, i + 2); const float x3 = transformed_value(input, i + 3); - // 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] = static_cast(std::round((x0 - min_val) * inv_delta)); - quantized[i + 1] = static_cast(std::round((x1 - min_val) * inv_delta)); - quantized[i + 2] = static_cast(std::round((x2 - min_val) * inv_delta)); - quantized[i + 3] = static_cast(std::round((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 uint32_t a0 = to_byte(x0); + const uint32_t a1 = to_byte(x1); + const uint32_t a2 = to_byte(x2); + const uint32_t a3 = to_byte(x3); + quantized[i] = static_cast(a0); + quantized[i + 1] = static_cast(a1); + quantized[i + 2] = static_cast(a2); + quantized[i + 3] = static_cast(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 += a0 * a0; + q1 += a1 * a1; + q2 += a2 * a2; + q3 += 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); + uint32_t q_sum_squares = (q0 + q1) + (q2 + q3); for (; i < this->dim; ++i) { const float x = transformed_value(input, i); - quantized[i] = static_cast(std::round((x - min_val) * inv_delta)); - sum += x; + const uint32_t a = to_byte(x); + quantized[i] = static_cast(a); + q_sum += a; if constexpr (Metric == VecSimMetric_L2) { - sum_squares += x * x; + q_sum_squares += a * a; } } - // 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; - MetadataType buf[5] = {min_val, delta, sum}; - size_t n = 3; + // Metadata is a packed run of 4-byte fields: min and delta as FP32, the quantized sums as + // uint32, and x_mean_ip (FP32) last for IP with norm. memcpy throughout because the offset + // (dim * sizeof(uint8_t)) is not guaranteed to be 4-byte aligned. + unsigned char buf[5 * sizeof(uint32_t)]; + size_t n = 0; + const auto put = [&buf, &n](const auto &value) { + static_assert(sizeof(value) == sizeof(uint32_t)); + memcpy(buf + n, &value, sizeof(value)); + n += sizeof(value); + }; + put(min_val); + put(delta); + put(q_sum); if constexpr (Metric == VecSimMetric_L2) - buf[n++] = sum_squares; + put(q_sum_squares); if constexpr (WithNorm && Metric == VecSimMetric_IP) - buf[n++] = x_mean_ip; - memcpy(meta_dst, buf, n * sizeof(MetadataType)); + put(x_mean_ip); + memcpy(quantized + this->dim, buf, n); } // Computes and writes query metadata (FP32) in a single pass over the query values. @@ -585,7 +623,14 @@ class QuantPreprocessor : public PreprocessorInterface { if constexpr (Metric == VecSimMetric_IP) x_mean_ip += input_value * mean[i]; } - return {min_val, max_val}; + // Centering is an FP32 subtraction of finite floats, so it can reach +/-inf: FLT_MAX + // against mean -FLT_MAX gives 6.8e38. Clamping the endpoints is both the storage limit + // (min is stored as FP32) and the guard quantize() relies on; only the endpoints are + // clamped, so the loop above stays FP32. std::clamp is safe on these two, which are + // finite or +/-inf but never NaN, since finite minus finite cannot be NaN. + constexpr float representable = std::numeric_limits::max(); + return {std::clamp(min_val, -representable, representable), + std::clamp(max_val, -representable, representable)}; } } diff --git a/src/VecSim/types/sq8.h b/src/VecSim/types/sq8.h index c1e9c40b8..70b710640 100644 --- a/src/VecSim/types/sq8.h +++ b/src/VecSim/types/sq8.h @@ -11,6 +11,7 @@ #include #include #include "VecSim/vec_sim_common.h" +#include "VecSim/utils/alignment.h" namespace vecsim_types { @@ -18,12 +19,23 @@ namespace vecsim_types { struct sq8 { using value_type = uint8_t; - // Metadata layout indices (stored after quantized values) + // Metadata layout indices (stored after quantized values). Each slot is 4 bytes. + // + // MIN_VAL and DELTA are FP32. Q_SUM and Q_SUM_SQUARES are uint32 sums over the *quantized* + // bytes, not over the input floats. That distinction is the whole point: a stored blob only + // ever describes its reconstruction x_r[i] = min + delta * a[i], so every quantity derived + // from it has to describe x_r as well. Summing the original input instead mixes two different + // vectors into one formula and produces distances that can even be negative. + // + // Integer sums are also exact, which matters more than it looks: the L2 identity + // ||x||^2 + ||y||^2 - 2*IP subtracts large nearly-equal quantities, so any rounding in the + // inputs lands directly on the answer. Keeping the sums exact and combining them in double + // (see the helpers below) is what makes that subtraction trustworthy. enum MetadataIndex : size_t { MIN_VAL = 0, DELTA = 1, - SUM = 2, - SUM_SQUARES = 3 // Only for L2 + Q_SUM = 2, // uint32: sum(a[i]) + Q_SUM_SQUARES = 3 // uint32: sum(a[i]^2). Only for L2 }; enum QueryMetadataIndex : size_t { @@ -57,6 +69,104 @@ struct sq8 { static constexpr size_t query_mean_ip_index() { return query_metadata_count() - 1; } + + // Readers for the packed metadata run that follows the dim quantized bytes. The run is not + // guaranteed to be 4-byte aligned, so every field goes through an unaligned load. Kernels use + // these rather than open-coding the offsets, so a layout change has one place to happen. + // Taking the metadata pointer, for kernels that have already formed it. + static float min_val_at(const value_type *meta) { + return load_unaligned(meta + MIN_VAL * sizeof(float)); + } + static float delta_at(const value_type *meta) { + return load_unaligned(meta + DELTA * sizeof(float)); + } + static uint32_t q_sum_at(const value_type *meta) { + return load_unaligned(meta + Q_SUM * sizeof(uint32_t)); + } + static uint32_t q_sum_squares_at(const value_type *meta) { + return load_unaligned(meta + Q_SUM_SQUARES * sizeof(uint32_t)); + } + + // Taking the blob pointer, for kernels that have not. + static float min_val_of(const value_type *blob, size_t dim) { return min_val_at(blob + dim); } + static float delta_of(const value_type *blob, size_t dim) { return delta_at(blob + dim); } + static uint32_t q_sum_of(const value_type *blob, size_t dim) { return q_sum_at(blob + dim); } + static uint32_t q_sum_squares_of(const value_type *blob, size_t dim) { + return q_sum_squares_at(blob + dim); + } + + // Single source of truth for turning stored metadata back into properties of the + // reconstruction. Every kernel derives its terms here so the storage layout has exactly one + // interpretation. All of these take double and return double on purpose: the integer sums are + // exact, and the point is to keep them exact through the combination. + + // sum(x_r[i]) where x_r[i] = min + delta * a[i]. + [[gnu::always_inline]] + static double reconstructed_sum(double min_val, double delta, uint32_t q_sum, size_t dim) { + return static_cast(dim) * min_val + delta * static_cast(q_sum); + } + + // sum(x_r[i]^2), expanded so that no term is ever formed from a rounded input. + [[gnu::always_inline]] + static double reconstructed_sum_squares(double min_val, double delta, uint32_t q_sum, + uint32_t q_sum_squares, size_t dim) { + return static_cast(dim) * min_val * min_val + + 2.0 * min_val * delta * static_cast(q_sum) + + delta * delta * static_cast(q_sum_squares); + } + + // ||x_r||^2 for one stored blob. Two entry points because some kernels have the blob pointer + // and some have already formed the metadata pointer. + [[gnu::always_inline]] + static double reconstructed_sum_squares_at(const value_type *meta, size_t dim) { + return reconstructed_sum_squares(min_val_at(meta), delta_at(meta), q_sum_at(meta), + q_sum_squares_at(meta), dim); + } + [[gnu::always_inline]] + static double reconstructed_sum_squares_of(const value_type *blob, size_t dim) { + return reconstructed_sum_squares_at(blob + dim, dim); + } + + // IP(x_r, y_r) for two quantized blobs, from their exact integer sums and integer dot product. + [[gnu::always_inline]] + static double reconstructed_ip(double min1, double delta1, uint32_t q_sum1, double min2, + double delta2, uint32_t q_sum2, uint64_t q_dot, size_t dim) { + return static_cast(dim) * min1 * min2 + + min1 * delta2 * static_cast(q_sum2) + + min2 * delta1 * static_cast(q_sum1) + + delta1 * delta2 * static_cast(q_dot); + } + + // ||x_r - y_r||^2 for two quantized blobs. Expanding the difference before summing keeps any + // common offset in the (min1 - min2) term, where it is formed once and exactly, instead of + // letting it inflate the summed terms and cancel away the answer. That is why this is not + // written as reconstructed_sum_squares + reconstructed_sum_squares - 2 * reconstructed_ip: + // algebraically identical, numerically not. + [[gnu::always_inline]] + static double reconstructed_l2_sqr(double min1, double delta1, uint32_t q_sum1, + uint32_t q_sum_squares1, double min2, double delta2, + uint32_t q_sum2, uint32_t q_sum_squares2, uint64_t q_dot, + size_t dim) { + // The quadratic part is regrouped so that the integer combination is formed first: + // 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. Two consequences: + // for two blobs sharing a delta the answer is delta^2 times a non-negative integer and so + // cannot come out negative, and for a blob against itself both factors are exactly zero, so + // the distance is exactly zero. Written as six independent floating point terms it was not: + // the compiler contracts some of them into fused multiply-adds and not others, so the + // products stop rounding identically and stop cancelling, which left about -1.7e-12 at + // dimension 512. Note this is deliberately not fixed by clamping at zero, which would also + // have hidden the much larger negative values the old metadata produced. + const int64_t sq_diff_sum = static_cast(q_sum_squares1) + + static_cast(q_sum_squares2) - + 2 * static_cast(q_dot); + const double quadratic = delta1 * delta2 * static_cast(sq_diff_sum) + + (delta1 - delta2) * (delta1 * static_cast(q_sum_squares1) - + delta2 * static_cast(q_sum_squares2)); + const double c = min1 - min2; + return static_cast(dim) * c * c + 2.0 * c * delta1 * static_cast(q_sum1) - + 2.0 * c * delta2 * static_cast(q_sum2) + quadratic; + } }; } // namespace vecsim_types diff --git a/tests/benchmark/spaces_benchmarks/bm_spaces_uint8.cpp b/tests/benchmark/spaces_benchmarks/bm_spaces_uint8.cpp index 602fff719..33f819936 100644 --- a/tests/benchmark/spaces_benchmarks/bm_spaces_uint8.cpp +++ b/tests/benchmark/spaces_benchmarks/bm_spaces_uint8.cpp @@ -31,12 +31,15 @@ class BM_VecSimSpaces_Integers_UINT8 : public benchmark::Fixture { test_utils::populate_uint8_vec(v2, dim, 1234); // Store the norm in the extra space for cosine calculations - *(float *)(v1 + dim) = test_utils::integral_compute_norm(v1, dim); - *(float *)(v2 + dim) = test_utils::integral_compute_norm(v2, dim); + // memcpy because v1 + dim is not guaranteed to be 4-byte aligned for arbitrary dim. + const float norm1 = test_utils::integral_compute_norm(v1, dim); + const float norm2 = test_utils::integral_compute_norm(v2, dim); + memcpy(v1 + dim, &norm1, sizeof(norm1)); + memcpy(v2 + dim, &norm2, sizeof(norm2)); } void TearDown(const ::benchmark::State &state) { - delete v1; - delete v2; + delete[] v1; + delete[] v2; } }; diff --git a/tests/unit/test_components.cpp b/tests/unit/test_components.cpp index efb39615d..983e826b3 100644 --- a/tests/unit/test_components.cpp +++ b/tests/unit/test_components.cpp @@ -8,6 +8,8 @@ */ #include +#include +#include #include "gtest/gtest.h" #include "VecSim/vec_sim.h" #include "VecSim/spaces/computer/preprocessor_container.h" @@ -1200,6 +1202,139 @@ TEST(PreprocessorsTest, QuantizationAsymmetricAlignment) { } } +// Finite input whose range is not representable in FP32. max - min overflowed to inf, which made +// delta inf and inv_delta 0, and inf * 0 is NaN; converting that NaN to an integer is undefined +// behaviour, so this reproduced UB on AddVector for input the API accepts. Runs under UBSan in CI. +TEST(PreprocessorsTest, QuantizationHandlesNonRepresentableRange) { + std::shared_ptr allocator = VecSimAllocator::newVecsimAllocator(); + constexpr size_t dim = 4; + constexpr unsigned char alignment = 0; + constexpr float huge = std::numeric_limits::max(); + float original_blob[dim] = {-huge, 0.0f, huge, 1.0f}; + + auto quant_preprocessor = + new (allocator) QuantPreprocessor(allocator, dim); + + void *storage_blob = nullptr; + void *query_blob = nullptr; + size_t storage_blob_size = dim * sizeof(float); + size_t query_blob_size = dim * sizeof(float); + quant_preprocessor->preprocess(original_blob, storage_blob, query_blob, storage_blob_size, + query_blob_size, alignment, alignment); + ASSERT_NE(storage_blob, nullptr); + + // The extremes must land on the ends of the range and everything must stay in [0, 255]. The + // middle values collapse to the same byte, which is expected: one byte cannot resolve a range + // this wide, and losing resolution is the correct outcome where UB was not. + const auto *quantized = static_cast(storage_blob); + EXPECT_EQ(quantized[0], 0); + EXPECT_EQ(quantized[2], 255); + + // Metadata must be finite, since the kernels derive every term from it. + const auto *metadata = quantized + dim; + EXPECT_TRUE(std::isfinite(load_unaligned(metadata + sq8::MIN_VAL * sizeof(float)))); + EXPECT_TRUE(std::isfinite(load_unaligned(metadata + sq8::DELTA * sizeof(float)))); + + allocator->free_allocation(storage_blob); + if (query_blob != storage_blob) { + allocator->free_allocation(query_blob); + } + delete quant_preprocessor; +} + +// The other direction of the same defect: a range too narrow to derive a usable delta from, rather +// than too wide. diff is not zero, so the min == max guard does not fire, but (float)(diff / 255) +// underflows to zero because the smallest positive float is 1.4e-45. 1 / 0 is inf, and the minimum +// element, whose numerator is exactly zero, then scales to 0 * inf = NaN, whose conversion to an +// integer is undefined behaviour. Runs under UBSan in CI. +TEST(PreprocessorsTest, QuantizationHandlesSubnormalRange) { + std::shared_ptr allocator = VecSimAllocator::newVecsimAllocator(); + constexpr size_t dim = 4; + constexpr unsigned char alignment = 0; + // 1e-44 is subnormal, so diff / 255 has no representable non-zero FP32 value. + float original_blob[dim] = {0.0f, 1e-44f, 0.0f, 1e-44f}; + + auto quant_preprocessor = + new (allocator) QuantPreprocessor(allocator, dim); + + void *storage_blob = nullptr; + void *query_blob = nullptr; + size_t storage_blob_size = dim * sizeof(float); + size_t query_blob_size = dim * sizeof(float); + quant_preprocessor->preprocess(original_blob, storage_blob, query_blob, storage_blob_size, + query_blob_size, alignment, alignment); + ASSERT_NE(storage_blob, nullptr); + + // delta must be a usable positive float, not zero, because every kernel divides by it or + // multiplies the sums by it. + const auto *quantized = static_cast(storage_blob); + const auto *metadata = quantized + dim; + const float delta = load_unaligned(metadata + sq8::DELTA * sizeof(float)); + EXPECT_TRUE(std::isfinite(delta)); + EXPECT_GT(delta, 0.0f); + EXPECT_TRUE(std::isfinite(load_unaligned(metadata + sq8::MIN_VAL * sizeof(float)))); + + // The range is far below one quantization step, so collapsing to a single byte is the correct + // outcome. Which byte is not the point; not being undefined behaviour is. + for (size_t i = 0; i < dim; ++i) { + EXPECT_EQ(quantized[i], quantized[0]) << "index " << i; + } + + allocator->free_allocation(storage_blob); + if (query_blob != storage_blob) { + allocator->free_allocation(query_blob); + } + delete quant_preprocessor; +} + +// The same defect reached through centering instead of through the raw range. find_min_max and +// transformed_value compute input[i] - mean[i] in FP32, so finite input against a finite mean can +// center to +/-inf: FLT_MAX against -FLT_MAX is 6.8e38. That inf used to reach delta, where +// 1 / inf = 0 turned the per-element scaling into inf * 0 = NaN. Widening the range to double does +// not help, because the value is already lost upstream in FP32. Runs under UBSan in CI. +TEST(PreprocessorsTest, QuantizationHandlesNonRepresentableCenteredRange) { + std::shared_ptr allocator = VecSimAllocator::newVecsimAllocator(); + constexpr size_t dim = 4; + constexpr unsigned char alignment = 0; + constexpr float huge = std::numeric_limits::max(); + + float original_blob[dim] = {huge, 0.0f, -huge, 1.0f}; + vecsim_stl::vector mean_vec(allocator); + // Centers element 0 to +inf and element 2 to -inf, from entirely finite operands. + mean_vec.push_back(-huge); + mean_vec.push_back(0.0f); + mean_vec.push_back(huge); + mean_vec.push_back(0.0f); + + auto quant_preprocessor = + new (allocator) QuantPreprocessor(allocator, dim, mean_vec); + + void *storage_blob = nullptr; + size_t storage_blob_size = dim * sizeof(float); + quant_preprocessor->preprocessForStorage(original_blob, storage_blob, storage_blob_size, + alignment); + ASSERT_NE(storage_blob, nullptr); + + // Metadata must be finite, since the kernels derive every term from it, and delta must be + // positive so the reconstruction is not degenerate. + const auto *quantized = static_cast(storage_blob); + const auto *metadata = quantized + dim; + const float min_val = load_unaligned(metadata + sq8::MIN_VAL * sizeof(float)); + const float delta = load_unaligned(metadata + sq8::DELTA * sizeof(float)); + EXPECT_TRUE(std::isfinite(min_val)); + EXPECT_TRUE(std::isfinite(delta)); + EXPECT_GT(delta, 0.0f); + + // The centered extremes land on the ends of the byte range, which is where they belong. The + // two middle values collapse together, as in the non-centered case: one byte cannot resolve a + // range this wide, and losing resolution is the correct outcome where UB was not. + EXPECT_EQ(quantized[0], 255); + EXPECT_EQ(quantized[2], 0); + + allocator->free_allocation(storage_blob); + delete quant_preprocessor; +} + // Test edge case where all entries are equal TEST(PreprocessorsTest, QuantizationTestAllEntriesEqual) { std::shared_ptr allocator = VecSimAllocator::newVecsimAllocator(); @@ -1233,12 +1368,12 @@ TEST(PreprocessorsTest, QuantizationTestAllEntriesEqual) { ASSERT_FLOAT_EQ(min_val, 3.5f); ASSERT_FLOAT_EQ(delta, 1.0f); - // Verify sum and sum_squares for L2 metric - float expected_sum = 3.5f * dim; - float expected_sum_squares = 3.5f * 3.5f * dim; - ASSERT_FLOAT_EQ(load_unaligned(metadata + sq8::SUM * sizeof(float)), expected_sum); - ASSERT_FLOAT_EQ(load_unaligned(metadata + sq8::SUM_SQUARES * sizeof(float)), - expected_sum_squares); + // The sums are over the quantized bytes, and every value collapsed to 0 here, so both are 0. + // The vector itself is still recovered exactly: reconstructed_sum gives dim * min = 3.5 * dim, + // which is what any consumer of the metadata actually asks for. + ASSERT_EQ(load_unaligned(metadata + sq8::Q_SUM * sizeof(uint32_t)), 0u); + ASSERT_EQ(load_unaligned(metadata + sq8::Q_SUM_SQUARES * sizeof(uint32_t)), 0u); + ASSERT_DOUBLE_EQ(sq8::reconstructed_sum(min_val, delta, 0u, dim), 3.5 * dim); // Reconstruct and verify: min + quantized * delta = 3.5 + 0 * 1 = 3.5 for (size_t i = 0; i < dim; ++i) { @@ -1444,6 +1579,13 @@ class QuantPreprocessorFP16MetricTest : public testing::TestWithParam(base) + byte_offset, sizeof(v)); + return v; + } + template void runQuantizationTest() { const size_t expected_storage_size = getExpectedStorageSize(); @@ -1459,9 +1601,21 @@ class QuantPreprocessorFP16MetricTest : public testing::TestWithParam(allocator, dim); @@ -1496,12 +1650,13 @@ class QuantPreprocessorFP16MetricTest : public testing::TestWithParamfree_allocation(storage_blob); @@ -1533,11 +1688,11 @@ class QuantPreprocessorFP16MetricTest : public testing::TestWithParam(static_cast(blob), original_blob, dim)); ASSERT_FLOAT_EQ(load_meta(blob, query_meta_offset + sq8::SUM_QUERY * sizeof(float)), - baseline_sum); + baseline_y_sum); if constexpr (Metric == VecSimMetric_L2) { ASSERT_FLOAT_EQ( load_meta(blob, query_meta_offset + sq8::SUM_SQUARES_QUERY * sizeof(float)), - baseline_sum_sq); + baseline_y_sum_sq); } allocator->free_allocation(blob); } diff --git a/tests/unit/test_spaces.cpp b/tests/unit/test_spaces.cpp index 8e53c83d1..d6d39cdf4 100644 --- a/tests/unit/test_spaces.cpp +++ b/tests/unit/test_spaces.cpp @@ -371,6 +371,104 @@ TEST_F(SpacesTest, SQ8_FP32_l2sqr_no_optimization_func_test) { ASSERT_NEAR(dist, baseline, 0.01) << "SQ8_FP32_L2Sqr failed to match expected distance"; } +// The three cases below are the defects this formulation exists to fix. They are written against +// an FP64 reference over the reconstructed vectors, which is the quantity the kernels are trying to +// compute, rather than against numbers recorded from a previous run. +namespace { +// ||x_r - y_r||^2 in FP64, where each operand is read back from its own quantized blob. +double ReferenceL2SqrFromBlobs(const uint8_t *a, const uint8_t *b, size_t dim) { + const double min_a = sq8::min_val_of(a, dim), delta_a = sq8::delta_of(a, dim); + const double min_b = sq8::min_val_of(b, dim), delta_b = sq8::delta_of(b, dim); + double acc = 0.0; + for (size_t i = 0; i < dim; i++) { + const double d = (min_a + delta_a * a[i]) - (min_b + delta_b * b[i]); + acc += d * d; + } + return acc; +} +} // namespace + +// A large common offset used to cancel the answer away entirely. Quantization is exact for this +// input (q = [0, 255]), so a wrong result is the summation form's fault, not the 8 bits'. +TEST_F(SpacesTest, SQ8_SQ8_L2_survives_large_common_offset) { + constexpr size_t dim = 2; + const std::vector x = {100000.0f, 100008.0f}; + const std::vector y = {100000.0f, 100000.0f}; + + const size_t blob_size = + dim * sizeof(uint8_t) + sq8::storage_metadata_count() * sizeof(float); + std::vector qx(blob_size), qy(blob_size); + test_utils::quantize_float_vec_to_sq8_with_metadata(x.data(), dim, qx.data()); + test_utils::quantize_float_vec_to_sq8_with_metadata(y.data(), dim, qy.data()); + + // 64 is the distance between the inputs; the reference sees 64.0000076 because reconstructing + // through FP32 min and delta does not land exactly on 100008. That residue is quantization, + // which is expected and tiny. What matters is that the kernel agrees with the reference rather + // than returning 0. + const double expected = ReferenceL2SqrFromBlobs(qx.data(), qy.data(), dim); + ASSERT_NEAR(expected, 64.0, 1e-4) << "the reference itself should see 64"; + + // Previously returned 0: ||x||^2 and 2*IP are both about 4e10, one FP32 step there is ~4768, + // and the answer is 64. + EXPECT_NEAR(SQ8_SQ8_L2Sqr(qx.data(), qy.data(), dim), expected, 1e-3); +} + +// A query that does not land on a quantization grid point. The stored sums used to be taken over +// the original input while the inner product described the reconstruction, mixing two vectors into +// one formula, which drove the result negative and made a radius-0 range query return a match. +TEST_F(SpacesTest, SQ8_FP32_L2_is_non_negative_for_off_grid_query) { + constexpr size_t dim = 3; + const std::vector stored = {0.0f, 0.25f, 1.0f}; + + const size_t blob_size = + dim * sizeof(uint8_t) + sq8::storage_metadata_count() * sizeof(float); + std::vector qs(blob_size); + test_utils::quantize_float_vec_to_sq8_with_metadata(stored.data(), dim, qs.data()); + + std::vector query(dim + sq8::query_metadata_count()); + query[0] = 0.0f; + query[1] = 0.2501f; + query[2] = 1.0f; + test_utils::preprocess_sq8_fp32_query(query.data(), dim); + + // Reference: reconstruction of the stored blob against the exact query values. + const double min_s = sq8::min_val_of(qs.data(), dim), delta_s = sq8::delta_of(qs.data(), dim); + double expected = 0.0; + for (size_t i = 0; i < dim; i++) { + const double d = (min_s + delta_s * qs[i]) - static_cast(query[i]); + expected += d * d; + } + + const float dist = SQ8_FP32_L2Sqr(qs.data(), query.data(), dim); + EXPECT_GE(dist, 0.0f) << "a squared distance must never be negative"; + EXPECT_NEAR(dist, expected, 1e-6); +} + +// High dimension, near-identical vectors: the regime where combining the exact integer sums in FP32 +// still loses the answer, because the three large sums cancel. +TEST_F(SpacesTest, SQ8_SQ8_L2_matches_fp64_reference_at_high_dim) { + constexpr size_t dim = 16384; + std::vector x(dim), y(dim); + std::mt19937 gen(4242); + std::uniform_real_distribution dist_gen(0.0f, 1.0f); + for (size_t i = 0; i < dim; i++) { + x[i] = dist_gen(gen); + y[i] = x[i] + (i % 7 == 0 ? 0.0005f : 0.0f); + } + + const size_t blob_size = + dim * sizeof(uint8_t) + sq8::storage_metadata_count() * sizeof(float); + std::vector qx(blob_size), qy(blob_size); + test_utils::quantize_float_vec_to_sq8_with_metadata(x.data(), dim, qx.data()); + test_utils::quantize_float_vec_to_sq8_with_metadata(y.data(), dim, qy.data()); + + const double expected = ReferenceL2SqrFromBlobs(qx.data(), qy.data(), dim); + ASSERT_GT(expected, 0.0); + const double got = SQ8_SQ8_L2Sqr(qx.data(), qy.data(), dim); + EXPECT_LT(std::abs(got - expected) / expected, 1e-5) + << "got " << got << " expected " << expected; +} + TEST_F(SpacesTest, SQ8_FP32_odd_dim_unaligned_metadata_test) { for (const size_t dim : {1UL, 5UL, 7UL, 15UL}) { const size_t query_size = dim + sq8::query_metadata_count(); @@ -498,15 +596,19 @@ TEST_F(SpacesTest, SQ8_FP16_l2sqr_odd_dim_unaligned_metadata_test) { auto store_float = [](uint8_t *dst, float value) { std::memcpy(dst, &value, sizeof(value)); }; + auto store_u32 = [](uint8_t *dst, uint32_t value) { std::memcpy(dst, &value, sizeof(value)); }; + + // min 0 and delta 1 make the reconstruction equal to the quantized bytes, so the integer sums + // carry the same values the FP32 sums used to: 1+..+5 and 1+4+9+16+25. constexpr float min_val = 0.0f; constexpr float delta = 1.0f; - constexpr float storage_sum = 15.0f; - constexpr float storage_sum_squares = 55.0f; + constexpr uint32_t storage_q_sum = 15; + constexpr uint32_t storage_q_sum_squares = 55; uint8_t *storage_meta = storage.data() + dim; store_float(storage_meta + sq8::MIN_VAL * sizeof(float), min_val); store_float(storage_meta + sq8::DELTA * sizeof(float), delta); - store_float(storage_meta + sq8::SUM * sizeof(float), storage_sum); - store_float(storage_meta + sq8::SUM_SQUARES * sizeof(float), storage_sum_squares); + store_u32(storage_meta + sq8::Q_SUM * sizeof(uint32_t), storage_q_sum); + store_u32(storage_meta + sq8::Q_SUM_SQUARES * sizeof(uint32_t), storage_q_sum_squares); for (size_t i = 0; i < dim; i++) { query[i] = vecsim_types::FP32_to_FP16(static_cast(i + 2)); @@ -518,7 +620,7 @@ TEST_F(SpacesTest, SQ8_FP16_l2sqr_odd_dim_unaligned_metadata_test) { store_float(query_meta + sq8::SUM_QUERY * sizeof(float), query_sum); store_float(query_meta + sq8::SUM_SQUARES_QUERY * sizeof(float), query_sum_squares); - const auto *storage_sum_squares_addr = storage_meta + sq8::SUM_SQUARES * sizeof(float); + const auto *storage_sum_squares_addr = storage_meta + sq8::Q_SUM_SQUARES * sizeof(uint32_t); const auto *query_sum_squares_addr = query_meta + sq8::SUM_SQUARES_QUERY * sizeof(float); ASSERT_NE(reinterpret_cast(storage_sum_squares_addr) % alignof(float), 0u); ASSERT_NE(reinterpret_cast(query_sum_squares_addr) % alignof(float), 0u); @@ -4767,3 +4869,169 @@ TEST_F(SpacesTest, SQ8_SQ8_DispatcherAlignmentHints) { check("Cosine", &spaces::Cosine_SQ8_SQ8_GetDistFunc); } #endif // CPU_FEATURES_ARCH_X86_64 + +// =================================================================================== +// Regression tests for the exact-metadata contract. Every one of these fails on the +// commit before this series and passes after it. They deliberately use only the +// quantizer helper, sq8::MIN_VAL, sq8::DELTA and the public kernels, so the same +// source compiles on both revisions and the before/after comparison is meaningful. +// =================================================================================== + +#include + +namespace { + +// min_val and delta occupy metadata slots 0 and 1 as FP32 in every revision. +float sq8_meta_float(const uint8_t *blob, size_t dim, size_t slot) { + float v; + std::memcpy(&v, blob + dim + slot * sizeof(float), sizeof(v)); + return v; +} + +// FP64 dot product of the two reconstructions, derived only from the bytes, min and delta. +// This is the quantity every SQ8-to-SQ8 kernel is trying to compute. +double sq8_reference_dot(const uint8_t *a, const uint8_t *b, size_t dim) { + const double mn_a = sq8_meta_float(a, dim, sq8::MIN_VAL); + const double d_a = sq8_meta_float(a, dim, sq8::DELTA); + const double mn_b = sq8_meta_float(b, dim, sq8::MIN_VAL); + const double d_b = sq8_meta_float(b, dim, sq8::DELTA); + double acc = 0.0; + for (size_t i = 0; i < dim; i++) { + acc += (mn_a + d_a * a[i]) * (mn_b + d_b * b[i]); + } + return acc; +} + +std::vector sq8_quantize_l2(const std::vector &v) { + const size_t dim = v.size(); + std::vector blob(dim * sizeof(uint8_t) + + sq8::storage_metadata_count() * sizeof(float)); + test_utils::quantize_float_vec_to_sq8_with_metadata(v.data(), dim, blob.data()); + return blob; +} + +} // namespace + +// A blob against itself must be exactly zero, on every kernel. Both factors of the regrouped +// quadratic term are exactly zero for identical blobs, and multiplying by an exact zero is exact, +// so this holds regardless of how the compiler contracts the surrounding arithmetic. Before the +// metadata change the residue was 4 * sum((x_r[i] - min) * rounding_error[i]), nonzero and of +// arbitrary sign. +TEST_F(SpacesTest, SQ8_SQ8_L2_self_distance_is_exactly_zero) { + std::mt19937 gen(20260812); + std::uniform_real_distribution value_gen(-3.0f, 5.0f); + + for (const size_t dim : {3UL, 5UL, 8UL, 64UL, 512UL}) { + std::vector x(dim); + for (size_t i = 0; i < dim; i++) { + x[i] = value_gen(gen); + } + const std::vector q = sq8_quantize_l2(x); + + EXPECT_EQ(SQ8_SQ8_L2Sqr(q.data(), q.data(), dim), 0.0f) << "scalar kernel, dim " << dim; + + // Also exercise whatever this CPU dispatches to. + unsigned char alignment = 0; + auto dispatched = L2_SQ8_SQ8_GetDistFunc(dim, &alignment, nullptr); + EXPECT_EQ(dispatched(q.data(), q.data(), dim), 0.0f) << "dispatched kernel, dim " << dim; + } +} + +// The same property with numbers small enough to check by hand. x = [0, 100.5, 255] gives min 0 and +// delta exactly 1, so the bytes are [0, 101, 255] and the stored vector is [0, 101, 255]. Summing +// the input gave x_sum_squares 75125.25 while the dot product of the bytes gives 75226, and the old +// formula used the sum of squares twice, so it returned 75125.25 + 75125.25 - 2 * 75226 = -201.5. +TEST_F(SpacesTest, SQ8_SQ8_L2_self_distance_whiteboard_case) { + const std::vector x = {0.0f, 100.5f, 255.0f}; + const std::vector q = sq8_quantize_l2(x); + + ASSERT_EQ(q[0], 0); + ASSERT_EQ(q[1], 101) << "expected 100.5 to round up to byte 101"; + ASSERT_EQ(q[2], 255); + + const float dist = SQ8_SQ8_L2Sqr(q.data(), q.data(), 3); + EXPECT_GE(dist, 0.0f) << "a squared distance must never be negative"; + EXPECT_EQ(dist, 0.0f); +} + +// Two inputs that quantize to the same bytes are the same stored vector, so they must answer a +// query identically. The distance has to be a function of the stored blob and nothing else. +TEST_F(SpacesTest, SQ8_SQ8_distance_ignores_sub_quantum_input_changes) { + constexpr size_t dim = 4; + // Same min and max, and both interior values perturbed well under half a delta and away from a + // rounding boundary, so both inputs land on the identical byte string. min is nonzero: with + // min == 0 the old inner product formula collapses to delta_x * delta_y * q_dot and the + // input-derived sums drop out, which would make the IP half of this test inert. + const std::vector x = {2.0f, 2.4f, 2.72f, 3.0f}; + const std::vector x_perturbed = {2.0f, 2.4005f, 2.7205f, 3.0f}; + const std::vector y = {1.5f, 1.6f, 2.4f, 2.5f}; + + const std::vector qx = sq8_quantize_l2(x); + const std::vector qxp = sq8_quantize_l2(x_perturbed); + const std::vector qy = sq8_quantize_l2(y); + + ASSERT_TRUE(std::equal(qx.begin(), qx.begin() + dim, qxp.begin())) + << "inputs no longer quantize identically; pick new values"; + ASSERT_EQ(sq8_meta_float(qx.data(), dim, sq8::MIN_VAL), + sq8_meta_float(qxp.data(), dim, sq8::MIN_VAL)); + ASSERT_EQ(sq8_meta_float(qx.data(), dim, sq8::DELTA), + sq8_meta_float(qxp.data(), dim, sq8::DELTA)); + + EXPECT_EQ(SQ8_SQ8_L2Sqr(qx.data(), qy.data(), dim), SQ8_SQ8_L2Sqr(qxp.data(), qy.data(), dim)); + EXPECT_EQ(SQ8_SQ8_InnerProduct(qx.data(), qy.data(), dim), + SQ8_SQ8_InnerProduct(qxp.data(), qy.data(), dim)); +} + +// The inner product must equal the dot product of the two reconstructions. +TEST_F(SpacesTest, SQ8_SQ8_IP_matches_reconstruction_dot) { + const std::vector> inputs = { + {-1.0f, 0.3f, 2.0f}, + {-2.5f, -1.1f, 0.4f, 0.9f, 3.3f, 7.7f, 8.0f, 9.9f}, + }; + for (const auto &x : inputs) { + const size_t dim = x.size(); + const std::vector q = sq8_quantize_l2(x); + const double expected = sq8_reference_dot(q.data(), q.data(), dim); + + // SQ8_SQ8_InnerProduct returns 1 - IP. + const double got = 1.0 - static_cast(SQ8_SQ8_InnerProduct(q.data(), q.data(), dim)); + EXPECT_LT(std::abs(got - expected) / expected, 1e-5) + << "dim " << dim << ": got " << got << " expected " << expected; + } +} + +// Plain uint8, not SQ8. The accumulated total is 255 * 255 * dim, which passes INT_MAX from +// dimension 33,026: the scalar path was signed-overflow UB there and the AVX512 and NEON L2 +// reduces read their unsigned total back as a signed int and went negative. All-255 bytes are the +// worst case and make the expected value an exact integer. +TEST_F(SpacesTest, UINT8_L2Sqr_and_InnerProduct_are_exact_past_int32) { + for (const size_t dim : {33026UL, 40000UL}) { + std::vector v1(dim + sizeof(float), 255); + std::vector v2(dim + sizeof(float), 0); + + // L2 between all-255 and all-0 is 255^2 * dim. + const double expected_l2 = 255.0 * 255.0 * static_cast(dim); + const float l2 = UINT8_L2Sqr(v1.data(), v2.data(), dim); + EXPECT_GT(l2, 0.0f) << "dim " << dim << ": squared distance went negative"; + EXPECT_LT(std::abs(static_cast(l2) - expected_l2) / expected_l2, 1e-6) + << "scalar L2, dim " << dim; + + unsigned char alignment = 0; + auto dispatched_l2 = L2_UINT8_GetDistFunc(dim, &alignment, nullptr); + const float l2_simd = dispatched_l2(v1.data(), v2.data(), dim); + EXPECT_GT(l2_simd, 0.0f) << "dim " << dim << ": SIMD squared distance went negative"; + EXPECT_LT(std::abs(static_cast(l2_simd) - expected_l2) / expected_l2, 1e-6) + << "dispatched L2, dim " << dim; + + // IP between two all-255 vectors is 255^2 * dim, and the kernel returns 1 - IP. + const double expected_ip = 1.0 - 255.0 * 255.0 * static_cast(dim); + const double ip = static_cast(UINT8_InnerProduct(v1.data(), v1.data(), dim)); + EXPECT_LT(std::abs(ip - expected_ip) / std::abs(expected_ip), 1e-6) + << "scalar IP, dim " << dim; + + auto dispatched_ip = IP_UINT8_GetDistFunc(dim, &alignment, nullptr); + const double ip_simd = static_cast(dispatched_ip(v1.data(), v1.data(), dim)); + EXPECT_LT(std::abs(ip_simd - expected_ip) / std::abs(expected_ip), 1e-6) + << "dispatched IP, dim " << dim; + } +} diff --git a/tests/unit/unit_test_utils.h b/tests/unit/unit_test_utils.h index 0039dea4b..fa8312595 100644 --- a/tests/unit/unit_test_utils.h +++ b/tests/unit/unit_test_utils.h @@ -250,20 +250,24 @@ inline void ComputeSQ8Quantization(const float *original_blob, size_t dim, uint8 float diff = max_val - min_val; float delta = (diff == 0.0f) ? 1.0f : diff / 255.0f; - // Calculate quantized values, sum and sum_squares - float sum = 0.0f; - float sum_squares = 0.0f; + // Calculate quantized values, and sum them as integers: the stored sums describe the + // reconstruction rather than the input, and are exact (see sq8.h). + uint32_t q_sum = 0; + uint32_t q_sum_squares = 0; for (size_t i = 0; i < dim; i++) { float normalized = (original_blob[i] - min_val) / delta; - output[i] = static_cast(std::round(normalized)); - sum += original_blob[i]; - sum_squares += original_blob[i] * original_blob[i]; + const uint32_t q = static_cast(std::round(normalized)); + output[i] = static_cast(q); + q_sum += q; + q_sum_squares += q * q; } - // Store metadata: min_val, delta, sum, sum_squares. Use memcpy because the metadata region - // (output + dim) is not guaranteed to be 4-byte aligned for arbitrary dim values. - const float metadata[4] = {min_val, delta, sum, sum_squares}; - std::memcpy(output + dim, metadata, sizeof(metadata)); + // Store metadata: min_val and delta as FP32, then the two integer sums. Use memcpy because the + // metadata region (output + dim) is not guaranteed to be 4-byte aligned for arbitrary dim. + const float head[2] = {min_val, delta}; + const uint32_t sums[2] = {q_sum, q_sum_squares}; + std::memcpy(output + dim, head, sizeof(head)); + std::memcpy(output + dim + sizeof(head), sums, sizeof(sums)); } // TODO: Move all test_utils to this namespace diff --git a/tests/utils/tests_utils.h b/tests/utils/tests_utils.h index 67d93177c..d7fa5ff62 100644 --- a/tests/utils/tests_utils.h +++ b/tests/utils/tests_utils.h @@ -181,31 +181,32 @@ static void quantize_float_vec_to_sq8_with_metadata(const float *v, size_t dim, max_val = std::max(max_val, v[i]); } - float sum = 0.0f; - float square_sum = 0.0f; - for (size_t i = 0; i < dim; i++) { - sum += v[i]; - square_sum += v[i] * v[i]; - } - // Calculate delta float delta = (max_val - min_val) / 255.0f; if (delta == 0) delta = 1.0f; // Avoid division by zero - // Quantize each value + // Quantize each value, summing the quantized bytes as we go. The sums describe the + // reconstruction rather than the input, and they are integers, matching the layout the kernels + // read (see sq8.h). + uint32_t q_sum = 0; + uint32_t q_sum_squares = 0; for (size_t i = 0; i < dim; i++) { float normalized = (v[i] - min_val) / delta; normalized = std::max(0.0f, std::min(255.0f, normalized)); - qv[i] = static_cast(std::round(normalized)); + const uint32_t q = static_cast(std::round(normalized)); + qv[i] = static_cast(q); + q_sum += q; + q_sum_squares += q * q; } - // Store parameters: [min, delta, sum, square_sum] + // Store parameters: [min, delta, q_sum, q_sum_squares] auto *params = qv + dim; std::memcpy(params + sq8::MIN_VAL * sizeof(float), &min_val, sizeof(float)); std::memcpy(params + sq8::DELTA * sizeof(float), &delta, sizeof(float)); - std::memcpy(params + sq8::SUM * sizeof(float), &sum, sizeof(float)); - std::memcpy(params + sq8::SUM_SQUARES * sizeof(float), &square_sum, sizeof(float)); + std::memcpy(params + sq8::Q_SUM * sizeof(uint32_t), &q_sum, sizeof(q_sum)); + std::memcpy(params + sq8::Q_SUM_SQUARES * sizeof(uint32_t), &q_sum_squares, + sizeof(q_sum_squares)); } // Preprocess fp32 query for SQ8 IP/Cosine/L2 space. @@ -315,10 +316,9 @@ static float SQ8_FP16_NotOptimized_L2Sqr(const void *pVect1v, const void *pVect2 // Storage and query metadata sit at byte offsets that are not guaranteed 4-byte aligned // for odd `dimension`; use load_unaligned to avoid alignment UB. - const float min_val = load_unaligned(pVect1 + dimension + sq8::MIN_VAL * sizeof(float)); - const float delta = load_unaligned(pVect1 + dimension + sq8::DELTA * sizeof(float)); - const float x_sum_sq = - load_unaligned(pVect1 + dimension + sq8::SUM_SQUARES * sizeof(float)); + const float min_val = sq8::min_val_of(pVect1, dimension); + const float delta = sq8::delta_of(pVect1, dimension); + const double x_sum_sq = sq8::reconstructed_sum_squares_of(pVect1, dimension); const auto *query_meta_bytes = pVect2 + dimension * sizeof(vecsim_types::float16); const float y_sum_sq = load_unaligned(query_meta_bytes + sq8::SUM_SQUARES_QUERY * sizeof(float)); @@ -330,7 +330,8 @@ static float SQ8_FP16_NotOptimized_L2Sqr(const void *pVect1v, const void *pVect2 const float dequantized = pVect1[i] * delta + min_val; ip += dequantized * vecsim_types::FP16_to_FP32(vecsim_types::float16{raw}); } - return x_sum_sq + y_sum_sq - 2.0f * ip; + return static_cast(x_sum_sq + static_cast(y_sum_sq) - + 2.0 * static_cast(ip)); } // Preprocess FP16 query for SQ8 IP/Cosine/L2 space.