Bug: Inconsistencies in Cosine Distance Calculation for Near-Zero Length Vectors about usearch HOT 16 OPEN

Hi @tos-kamiya! Documenting the expected behavior is indeed important, I'll try to improve there. Introducing an additional runtime epsilon parameter to distance functions, is, however, infeasible. All of the metrics must have identical signature. We just need to find the smallest constant that results in a non-zero square root and use it as an epsilon.

from usearch.

@lifthrasiir, there shouldn’t be compatibility or any other issues. This class was significantly refactored in the last releases when the macro condition was broken. Prior to this, it definitely worked, and we have enough benchmarks coverage in SimSIMD to suggest improvements over autovectorized code.

As for fast-math settings, I agree, that with SimSIMD back ON, there shouldn’t be anything left to gain from that flag 🤗

from usearch.

@tos-kamiya, thanks for opening the issue! Can you please attach the outputs of print(index) - it will show details relevant for diagnostics?

I'd also recommend to use double-precision math if you want to work with such values, simply try passing "f64" in the constructor 🤗

from usearch.

Thank you for your prompt response to my previous message. Following your suggestion, I have conducted further tests and am providing the requested details below.

The output of print(index) using the original code (f32 data type) is as follows:

usearch.Index(ScalarKind.F32 x 2, MetricKind.Cos, connectivity: 16, expansion: 128 & 64, 20 vectors in 2 levels, avx512 hardware acceleration)

I modified the code to use double-precision (f64), changing the range of the vectors to 1.0e-140 to 1.0e-180. The results were as follows:

• With this higher precision, distances for vectors in the range 1.0e-140 to 1.0e-154 consistently returned 0.
• Vectors from 1.0e-163 to 1.0e-179 yielded distances of 1.
• Other vectors in between these ranges showed varying negative or positive distances close to 0.

Modified code:

import ast

import numpy as np
from usearch.index import Index, Matches

index = Index(
    ndim=2, # Define the number of dimensions in input vectors
    metric='cos', # Choose 'l2sq', 'haversine' or other metric, default = 'ip'
    dtype='f64', # Quantize to 'f16' or 'i8' if needed, default = 'f32'
)

kvs = {}
for k in range(140, 180):
    v = np.array([ast.literal_eval("1.0e-%d" % k), 0.0])
    kvs[k] = v
    index.add(k, v)

query = np.array([1.0, 0.0])

matches: Matches = index.search(query, 100)
for m in matches:
    print("distance %s to %s %s: %g" % (query, m.key, kvs[m.key], m.distance))

print(matches)

print(index)

Output:

distance [1. 0.] to 161 [1.e-161 0.e+000]: -0.00598771
distance [1. 0.] to 160 [1.e-160 0.e+000]: -5.56646e-06
distance [1. 0.] to 159 [1.e-159 0.e+000]: -6.25753e-07
distance [1. 0.] to 158 [1.e-158 0.e+000]: -8.17014e-09
distance [1. 0.] to 157 [1.e-157 0.e+000]: -1.80596e-11
distance [1. 0.] to 156 [1.e-156 0.e+000]: -7.67386e-13
distance [1. 0.] to 155 [1.e-155 0.e+000]: -1.55431e-15
distance [1. 0.] to 140 [1.e-140 0.e+000]: 0
distance [1. 0.] to 141 [1.e-141 0.e+000]: 0
...
distance [1. 0.] to 153 [1.e-153 0.e+000]: 0
distance [1. 0.] to 154 [1.e-154 0.e+000]: 0
distance [1. 0.] to 179 [1.e-179 0.e+000]: 1
distance [1. 0.] to 178 [1.e-178 0.e+000]: 1
...
distance [1. 0.] to 171 [1.e-171 0.e+000]: 1
distance [1. 0.] to 163 [1.e-163 0.e+000]: 1
usearch.Matches(40)
usearch.Index(ScalarKind.F64 x 2, MetricKind.Cos, connectivity: 16, expansion: 128 & 64, 40 vectors in 2 levels, avx512 hardware acceleration)

I hope this additional information helps in diagnosing the issue.

Thank you once again for your assistance.

from usearch.

@tos-kamiya thank you! Now the question is - how far do we want to go to achieve numerical stability.
It’s definitely possible, but may sacrifice accuracy. Coincidentally, a similar issue was raised today in SimSIMD. The issue originates here for f32 and here for f64.

Any chance you can help tune those epsilon params? That would help up to a million devices running those libraries 🤗

from usearch.

I've been reflecting on the discussion about the behavior of cosine distances with small norm vectors. Here are some thoughts:

• Inf Values in f32 Computations: It seems that the inf values are likely due to dividing by zero when normalizing small vectors. When the norm of a vector is essentially zero, the normalization process leads to division by zero, resulting in inf values.

• Handling Zero-Norm Vectors: For vectors that effectively have a zero norm, it might be more practical to directly set the cosine value (inner product) to zero instead of performing normalization that leads to inf values. This approach could prevent the generation of vectors with inf values.

• Behavior in f64 Computations: The jump from 0 to 1 in cosine distances as vectors get smaller in f64 seems consistent. As the distance increases, it implies less similarity, which aligns with the expected behavior in practical applications.

• Documentation Suggestion: I recommend documenting that for indexes with a cos metric, it's advisable not to add vectors with small norms to avoid these issues.

• Epsilon Parameter: While I'm not entirely clear about the epsilon parameter you mentioned, introducing an eps parameter in the Index constructor could be a solution. This parameter would set the cosine to zero for vectors with norms smaller than eps, ensuring stability in calculations.

I hope these observations and suggestions are helpful for addressing the issue.

from usearch.

@ashvardanian would it be a problem to make the choice of epsilon data-dependent? And use strong heuristics to determine which particular epsilon to use? (This could perhaps be determined by doing some brute-force/fuzz-testing-style tests.)

from usearch.

Yes, @turian, that might be problematic. We try to perform all the computations in a single forward pass over the vectors. Introducing the second pass to analyze the data would be too costly. I believe we only need to tune that constant a bit, reverse-engineering the minimum value, that after the rsqrt won't result in inf.

from usearch.

After some fiddling, I believe this bug is nothing to do with epsilon in SIMSIMD and actually a bigger bug in disguise: most SIMSIMD functions are not actually used because autovectorized functions can overtake them due to a wrong conditional. Therefore what's actually being called is metric_cos_gt::operator() which lacks any guard against near-zero norms. This can be also verified by comparing the reported distance with direct outputs from the simsimd package. Fixing either one---metric_cos_gt or metric_punned_t---does give a reasonable result.

from usearch.

@lifthrasiir it depends on how you compile the package. In Python, you can use index.hardware_acceleration to check if SimSIMD is used. On most modern hardware that should be the case.

Are you compiling from source? What hardware are you using? Have you passed the compilation parameter to enable SimSIMD? Check CONTRIBUTING.md for details on those.

from usearch.

My test environment was:

• Xeon Gold 6348H @ 2.38GHz
• Ubuntu 20.04.6 LTS
• Python 3.10.13

Under this environment I confirmed that index.hardware_acceleration indeed equals "avx512", which should mean USEARCH_USE_SIMSIMD has been correctly set.

Of course this doesn't actually mean that the corresponding SimSIMD code is being used---I did confirm this by slightly altering the source code and running pip install -e . as suggested. (The current SimSIMD doesn't actually compile clean on Ubuntu 20.04, so I had to manually comment some f16 functions out. This shouldn't affect outcomes for f32 or f64 though.)

from usearch.

@ashvardanian I was talking about a batch offline developer-only step, not a second pass during inference. Almost like a fuzz testing suite, but that generates appropriate float heuristics for the codebase.

from usearch.

@turian, thats worth considering. Let me know if there is a PoC solution you can contribute 🤗

from usearch.

I guess my comment was too succinct to illustrate the evident problem, so I'll elaborate:

    inline metric_punned_t(                                //
        std::size_t dimensions,                            //
        metric_kind_t metric_kind = metric_kind_t::l2sq_k, //
        scalar_kind_t scalar_kind = scalar_kind_t::f32_k) noexcept
        : raw_arg3_(dimensions), raw_arg4_(dimensions), dimensions_(dimensions), metric_kind_(metric_kind),
          scalar_kind_(scalar_kind) {

#if USEARCH_USE_SIMSIMD
        if (!configure_with_simsimd())
            configure_with_auto_vectorized();
#endif
        configure_with_auto_vectorized();

        if (scalar_kind == scalar_kind_t::b1x8_k)
            raw_arg3_ = raw_arg4_ = divide_round_up<CHAR_BIT>(dimensions_);
    }

I hope it's clear enough that configure_with_auto_vectorized can get called twice when USEARCH_USE_SIMSIMD is defined, and configure_with_auto_vectorized will always overwrite raw_ptr_.

Even when you don't have this fact beforehand (like, when I started out), distance([1, 0], [1e-10, 0]) returning 0 should have been already suspicious. If the SimSIMD routine was indeed in use, the calculation should go like:

ab = 1 * 1e-10 = 1e-10
a2 = 1 * 1 = 1
b2 = 1e-10 * 1e-10 ≈ 1e-20
rsqrt_a2 = rsqrt14(1 + 1e-9) ≈ rsqrt14(1)
rsqrt_b2 = rsqrt14(1e-20 + 1e-9) ≈ rsqrt14(1e-9)
distance ≈ 1 - 1e-10 * rsqrt14(1) * rsqrt14(1e-9) ≈ 1 - 1e-10 * 3.16e4 ≈ 1

In the other words, there is absolutely no reason for this to return 0 or infinity even with a presence of floating-point errors because the epsilon should have masked near-zero values out. I have also manually verified that the rsqrt14 approximation itself can't cause this, because its implementation is independent of the binary exponent (e.g. rsqrt14(1.5) and rsqrt14(6) would do the identical calculations). That's why I came to consider a possibility that this code was not actually used at all.

In comparison, the metric_cos_gt code would just return 1 - ab / (sqrt(a2) * sqrt(b2)) = 1 - 1e-10 / (sqrt(1) * sqrt(1e-20)) ≈ 1 - 1e-10 / 1e-10 = 0, as observed. It is even possible to sqrt(a2) * sqrt(b2) to underflow without a2 or b2 are not (yet) zero, especially when -ffast-math is set and FTZ/DAZ might be in effect, so it can also return both signs of infinities depending on the sign of ab. (In my opinion the whole zero-detection code is flawed for this reason.)

from usearch.

Oh, thank you, @lifthrasiir - I didn't realize you meant the macro-conditional, my bad 🤦‍♂️ Would you like to author a patch PR since you are the one who noticed that?

from usearch.

@ashvardanian Glad it worked this time :-p I think it will have a performance impact, which may well be negative, in addition to obvious incompatibility issues (for example, a saved database will probably no longer work I guess?). I can commentate evident bugs but can't decide what to do now---I even don't use usearch myself---, so it's probably better for you to fix this with all decision makings.

(Aside: did you check if -ffast-math is correctly localized to the usearch library proper by the way? Because glibc apparently enables FTZ and DAZ as a part of its C runtime and can have a global effect AFAIK. And my cursory reading didn't detect any complex calculation that hasn't been rewritten and can be benefited from -ffast-math.)

from usearch.