Avoid overflow and underflow in cosine similarity for finite vectors

algorithms/math/cosine_similarity.py

@@ -15,27 +15,27 @@
 from __future__ import annotations
 
 import math
+from collections.abc import Iterable
 
 
-def _l2_distance(vec: list[float]) -> float:
+def _l2_distance(vec: Iterable[float]) -> float:
     """Calculate the L2 (Euclidean) norm of a vector.
 
     Args:
-        vec: Input vector as a list of numbers.
+        vec: Input vector as an iterable of numbers.
 
     Returns:
         The L2 norm of the vector.
     """
-    norm = 0.0
-    for element in vec:
-        norm += element * element
-    norm = math.sqrt(norm)
-    return norm
+    return math.sqrt(math.fsum(element * element for element in vec))
 
 
 def cosine_similarity(vec1: list[float], vec2: list[float]) -> float:
     """Calculate cosine similarity between two vectors.
 
+    Scale each vector before multiplying components to avoid overflow and
+    underflow for very large or very small finite inputs.
+
     Args:
         vec1: First vector.
         vec2: Second vector (must be same length as vec1).
@@ -58,14 +58,12 @@ def cosine_similarity(vec1: list[float], vec2: list[float]) -> float:
             + str(len(vec2))
         )
 
-    norm_a = _l2_distance(vec1)
-    norm_b = _l2_distance(vec2)
-
-    similarity = 0.0
-
-    for vec1_element, vec2_element in zip(vec1, vec2, strict=False):
-        similarity += vec1_element * vec2_element
-
-    similarity /= norm_a * norm_b
-
-    return similarity
+    scale_a = max((abs(element) for element in vec1), default=1.0)
+    scale_b = max((abs(element) for element in vec2), default=1.0)
+    norm_a = _l2_distance(element / scale_a for element in vec1)
+    norm_b = _l2_distance(element / scale_b for element in vec2)
+    similarity = math.fsum(
+        (left / scale_a) * (right / scale_b)
+        for left, right in zip(vec1, vec2, strict=False)
+    )
+    return similarity / (norm_a * norm_b)

tests/test_cosine_scaling.py

@@ -0,0 +1,48 @@
+"""Cosine similarity is invariant under independent positive scaling."""
+
+import math
+import sys
+
+import pytest
+
+from algorithms.math.cosine_similarity import cosine_similarity
+
+
+@pytest.mark.parametrize("left_scale", [1e-200, 1.0, 1e200])
+@pytest.mark.parametrize("right_scale", [1e-200, 1.0, 1e200])
+def test_independent_vector_scales(left_scale, right_scale):
+    """Same, opposite, perpendicular and general directions retain their angle."""
+    cases = [
+        ([1, 2, 3], [1, 2, 3], 1.0),
+        ([1, 2, 3], [-1, -2, -3], -1.0),
+        ([1, 2, 3], [2, -1, 0], 0.0),
+        ([1, 2, 3], [1, -1, 2], 5 / math.sqrt(84)),
+    ]
+    for left, right, expected in cases:
+        scaled_left = [value * left_scale for value in left]
+        scaled_right = [value * right_scale for value in right]
+        result = cosine_similarity(scaled_left, scaled_right)
+        assert math.isfinite(result)
+        assert result == pytest.approx(expected, abs=1e-15)
+
+
+@pytest.mark.parametrize("magnitude", [sys.float_info.max, math.ulp(0.0)])
+def test_extreme_representable_components(magnitude):
+    """Even a nonrepresentable raw norm must not spoil a representable cosine."""
+    same = cosine_similarity([magnitude, magnitude], [magnitude, magnitude])
+    orthogonal = cosine_similarity([magnitude, magnitude], [magnitude, -magnitude])
+    assert same == pytest.approx(1)
+    assert orthogonal == pytest.approx(0)
+
+
+@pytest.mark.parametrize("left,right", [([], []), ([0, 0], [1, 2]), ([1, 2], [0, 0])])
+def test_zero_vector_exception_is_preserved(left, right):
+    """Zero vectors still have undefined similarity and raise the existing error."""
+    with pytest.raises(ZeroDivisionError):
+        cosine_similarity(left, right)
+
+
+def test_dimension_validation_is_preserved():
+    """Mismatched vector lengths keep their explicit validation error."""
+    with pytest.raises(ValueError, match="same length"):
+        cosine_similarity([1], [1, 2])