From 0a79024ea8075760eae9c7b1ec5c42f05ba649d9 Mon Sep 17 00:00:00 2001 From: Joonas Lipping Date: Fri, 16 Oct 2026 08:11:22 +0300 Subject: [PATCH 1/2] fix(statistics): guard against division by zero in correlation() When both variances are near zero, sqrt(variance_x * variance_y) can underflow to zero and cause a division by zero. Clamp the divisor to epsilon instead. --- dlib/statistics/statistics.h | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/dlib/statistics/statistics.h b/dlib/statistics/statistics.h index 914f5482ff..db52d3cbbc 100644 --- a/dlib/statistics/statistics.h +++ b/dlib/statistics/statistics.h @@ -363,8 +363,11 @@ namespace dlib << "\n\tyou have to add some numbers to this object first" << "\n\tthis: " << this ); - - return covariance() / std::sqrt(variance_x()*variance_y()); + // NOTE: The variance is always positive + auto const variance = std::sqrt(variance_x() * variance_y()); + return covariance() / (variance < std::numeric_limits::epsilon() + ? std::numeric_limits::epsilon() + : variance); } T variance_x ( From 0266408ef52b1271dd2ecfb4d03882bb8c4e9a35 Mon Sep 17 00:00:00 2001 From: Davis King Date: Fri, 16 Oct 2026 22:08:21 -0400 Subject: [PATCH 2/2] Fix correlation normalization without epsilon regularization --- dlib/statistics/statistics.h | 14 ++- dlib/statistics/statistics_abstract.h | 2 + dlib/test/statistics.cpp | 140 ++++++++++++++++++++++++++ 3 files changed, 151 insertions(+), 5 deletions(-) diff --git a/dlib/statistics/statistics.h b/dlib/statistics/statistics.h index db52d3cbbc..8069333754 100644 --- a/dlib/statistics/statistics.h +++ b/dlib/statistics/statistics.h @@ -363,11 +363,15 @@ namespace dlib << "\n\tyou have to add some numbers to this object first" << "\n\tthis: " << this ); - // NOTE: The variance is always positive - auto const variance = std::sqrt(variance_x() * variance_y()); - return covariance() / (variance < std::numeric_limits::epsilon() - ? std::numeric_limits::epsilon() - : variance); + const T vx = variance_x(); + const T vy = variance_y(); + if (vx == 0 || vy == 0) + return 0; + + // Avoid forming a product that can underflow or overflow. Divide by + // the smaller standard deviation first to avoid underflowing an + // intermediate quotient when the scales are very different. + return covariance() / std::sqrt(std::min(vx,vy)) / std::sqrt(std::max(vx,vy)); } T variance_x ( diff --git a/dlib/statistics/statistics_abstract.h b/dlib/statistics/statistics_abstract.h index 24432ded26..0955ac3ee7 100644 --- a/dlib/statistics/statistics_abstract.h +++ b/dlib/statistics/statistics_abstract.h @@ -47,6 +47,7 @@ namespace dlib ensures - returns the correlation coefficient between all the elements of a and b. (i.e. how correlated is a(i) with b(i)) + Returns 0 if either computed sample variance is 0. !*/ // ---------------------------------------------------------------------------------------- @@ -437,6 +438,7 @@ namespace dlib ensures - returns the correlation coefficient between all the x and y samples presented to this object via add() + - returns 0 if variance_x() == 0 or variance_y() == 0. !*/ T variance_x ( diff --git a/dlib/test/statistics.cpp b/dlib/test/statistics.cpp index f395a55bd7..15171e5cd7 100644 --- a/dlib/test/statistics.cpp +++ b/dlib/test/statistics.cpp @@ -13,6 +13,7 @@ #include #include #include +#include #include #include "tester.h" @@ -307,6 +308,141 @@ namespace } + template + void test_correlation_scales ( + const T x_scale, + const T y_scale + ) + { + running_scalar_covariance positive, negative, partial; + positive.add(-x_scale, -y_scale); + positive.add(0, 0); + positive.add(x_scale, y_scale); + negative.add(-x_scale, y_scale); + negative.add(0, 0); + negative.add(x_scale, -y_scale); + partial.add(-x_scale, 0); + partial.add(0, -y_scale); + partial.add(x_scale, y_scale); + + const T tolerance = 8*std::numeric_limits::epsilon(); + DLIB_TEST_MSG(std::abs(positive.correlation()-1) < tolerance, + "x_scale: " << x_scale << " y_scale: " << y_scale + << " correlation: " << positive.correlation()); + DLIB_TEST_MSG(std::abs(negative.correlation()+1) < tolerance, + "x_scale: " << x_scale << " y_scale: " << y_scale + << " correlation: " << negative.correlation()); + DLIB_TEST_MSG(std::abs(partial.correlation()-T(0.5)) < tolerance, + "x_scale: " << x_scale << " y_scale: " << y_scale + << " correlation: " << partial.correlation()); + } + + template + void test_running_scalar_correlation() + { + // The middle scales make the product of the variances underflow or + // overflow, while each variance and covariance remains representable. + const T scales[] = { + T(1), T(1e-4), + std::ldexp(T(1), std::numeric_limits::min_exponent/3), + std::ldexp(T(1), std::numeric_limits::max_exponent/3), + std::sqrt(std::numeric_limits::min()), + std::ldexp(T(1), (std::numeric_limits::max_exponent-4)/2) + }; + for (const T x_scale : scales) + { + for (const T y_scale : scales) + test_correlation_scales(x_scale, y_scale); + } + + if (std::numeric_limits::has_denorm == std::denorm_present) + { + // Use enough subnormal units to represent the covariance for + // the partial-correlation case as well as the variances. + const T scale = 4*std::sqrt(std::numeric_limits::denorm_min()); + test_correlation_scales(scale, scale); + + running_scalar_covariance smallest; + const T tiny = std::sqrt(std::numeric_limits::denorm_min()); + smallest.add(-tiny, -tiny); + smallest.add(0, 0); + smallest.add(tiny, tiny); + DLIB_TEST(std::abs(smallest.correlation()-1) < 8*std::numeric_limits::epsilon()); + + // The moments are 2, 1 and -1 subnormal units, respectively. + // Multiplying the standard deviations would round their product + // to 1 unit and incorrectly return -1 instead of -1/sqrt(2). + running_scalar_covariance subnormal; + subnormal.add(-2*tiny, tiny); + subnormal.add(0, tiny); + subnormal.add(0, -tiny); + subnormal.add(2*tiny, -tiny); + subnormal.add(0, 0); + DLIB_TEST(std::abs(subnormal.correlation()+1/std::sqrt(T(2))) < + 8*std::numeric_limits::epsilon()); + } + + running_scalar_covariance largest, weak; + const T large = std::nextafter(std::sqrt(std::numeric_limits::max()/2), T(0)); + largest.add(-large, -large); + largest.add(large, large); + DLIB_TEST(std::abs(largest.correlation()-1) < 8*std::numeric_limits::epsilon()); + + // Dividing by the larger standard deviation first would underflow + // the first quotient and lose this nonzero correlation. + const T expected = std::ldexp(T(1), std::numeric_limits::min_exponent+8); + const T big = scales[3]; + const T small = 1/big; + const T weak_signal = expected*big; + weak.add(big, 0); + weak.add(weak_signal, small); + weak.add(-big, 0); + weak.add(-weak_signal, -small); + DLIB_TEST(std::abs(weak.correlation()/expected-1) < 8*std::numeric_limits::epsilon()); + + running_scalar_covariance two_samples, constant_x, constant_y, constant_both; + const T scale = scales[2]; + two_samples.add(0, 0); + two_samples.add(scale, scale); + DLIB_TEST(std::abs(two_samples.correlation()-1) < 8*std::numeric_limits::epsilon()); + for (const T x : {T(-1), T(0), T(1)}) + { + constant_x.add(1, x); + constant_y.add(x, 1); + constant_both.add(1, 2); + } + DLIB_TEST(constant_x.correlation() == 0); + DLIB_TEST(constant_y.correlation() == 0); + DLIB_TEST(constant_both.correlation() == 0); + } + + void test_vector_correlation() + { + const double scales[] = {1e-9, 1e-100, 1e100}; + for (const double x_scale : scales) + { + for (const double y_scale : scales) + { + const std::vector x = {-x_scale, 0, x_scale}; + const std::vector y = {-y_scale, 0, y_scale}; + const std::vector negative = {y_scale, 0, -y_scale}; + const std::vector partial = {0, -y_scale, y_scale}; + DLIB_TEST(std::abs(correlation(x,y)-1) < 1e-14); + DLIB_TEST(std::abs(r_squared(x,y)-1) < 1e-14); + DLIB_TEST(std::abs(correlation(x,negative)+1) < 1e-14); + DLIB_TEST(std::abs(r_squared(x,negative)-1) < 1e-14); + DLIB_TEST(std::abs(correlation(x,partial)-0.5) < 1e-14); + DLIB_TEST(std::abs(r_squared(x,partial)-0.25) < 1e-14); + } + } + const std::vector constant(3, 1); + const std::vector varying = {-1, 0, 1}; + DLIB_TEST(correlation(constant,varying) == 0); + DLIB_TEST(correlation(varying,constant) == 0); + DLIB_TEST(correlation(constant,constant) == 0); + DLIB_TEST(r_squared(constant,varying) == 0); + } + void test_running_stats() { print_spinner(); @@ -943,6 +1079,10 @@ namespace test_running_cross_covariance(); test_running_cross_covariance_sparse(); test_running_stats(); + test_running_scalar_correlation(); + test_running_scalar_correlation(); + test_running_scalar_correlation(); + test_vector_correlation(); test_skewness_and_kurtosis_1(); test_skewness_and_kurtosis_2(); test_randomize_samples();