From 20140d34d6c207755af73e9c10ba2253c7bb776d Mon Sep 17 00:00:00 2001 From: Ahmed Eldeeb <62363199+deeb01@users.noreply.github.com> Date: Sun, 20 Sep 2026 14:49:39 -0700 Subject: [PATCH 1/3] Warn when the convolutional barycenter kernel underflows (closes #458) The 2D convolutional barycenters build K = exp(-(x-y)**2 / reg) on a unit grid. For small reg the distant entries underflow, so mass can no longer be moved across the image and the barycenter collapses towards the arithmetic mean of the inputs. That reads as an over-diffuse result rather than a failure, which is what #458 reports. Measured on two Gaussians of equal width separated by 0.15, comparing the recovered width against the true one: the default solver is 236% off at reg=5e-4, 182% at 1e-3, 20% at 2e-3 and correct from 4e-3 up, while method='sinkhorn_log' is exact throughout. It is not a convergence problem, 200000 iterations at stopThr=1e-12 return the same wrong answer, and it is not the debiasing, which reproduces the barycenter of {mu, mu} to 4e-4. Sinkhorn forms products and ratios of kernel entries, so only about half the exponent range is usable. The warning triggers when the smallest kernel exponent falls below log(tiny)/2, which lands where the measured error stops being negligible. test_convolutional_barycenter_non_square used reg=1e-3, which underflows on any grid. Its uniform input needs no transport so the assertion held, but the kernel was degenerate; it now uses 1e-2 and exercises the same code path. --- RELEASES.md | 1 + ot/bregman/_convolutional.py | 30 +++++++++++++++++ test/test_bregman.py | 65 ++++++++++++++++++++++++++++++++++-- 3 files changed, 94 insertions(+), 2 deletions(-) diff --git a/RELEASES.md b/RELEASES.md index 063701229..43e8f2d7f 100644 --- a/RELEASES.md +++ b/RELEASES.md @@ -15,6 +15,7 @@ #### Closed issues +- Warn when the convolution kernel in `ot.bregman.convolutional_barycenter2d` and `convolutional_barycenter2d_debiased` underflows at small `reg`. Past that point mass can no longer cross the image and the barycenter collapses towards the arithmetic mean of the inputs, which looks like an over-diffuse result rather than an error; `method='sinkhorn_log'` is exact in that regime (PR #867, Issue #458) - Remove a leftover debug `print` from `ot.utils.projection_sparse_simplex` with `axis=1`, and make the `ot.datasets.make_gauss_hd` docstring a raw string so importing `ot` no longer emits a `SyntaxWarning` (PR #860) - Fix `ot.dist` ignoring the weights `w` for `metric="cityblock"`, which returned the unweighted distance although the weights are documented for this metric (PR #859) - Fix swapped arguments to `div_to_product` in `ot.gromov.fused_unbalanced_across_spaces_cost`: with `reg_type="independent"` (UCOOT) the entropic terms used the plan marginals as the reference measures and vice versa (PR #855, Issue #854) diff --git a/ot/bregman/_convolutional.py b/ot/bregman/_convolutional.py index 9a8253240..3c4842c42 100644 --- a/ot/bregman/_convolutional.py +++ b/ot/bregman/_convolutional.py @@ -10,6 +10,8 @@ import warnings +import numpy as np + from ..backend import get_backend from ..utils import list_to_array @@ -20,6 +22,33 @@ ) +def _warn_if_kernel_underflows(nx, M1, M2, reg): + """Warn when exp(M) is about to lose the transport to underflow. + + The Sinkhorn iterations form products and ratios of kernel entries, so the + usable exponent range is roughly half that of the float type. Past that the + convolution can no longer move mass across the image and the barycenter + degenerates towards the arithmetic mean of the inputs, which looks like an + over-diffuse result rather than an error. + """ + min_exponent = float(min(nx.min(M1), nx.min(M2))) + try: + tiny = np.finfo(nx.to_numpy(M1).dtype).tiny + except (TypeError, ValueError): # pragma: no cover - exotic dtypes + return + # half the exponent range, i.e. the exponent of sqrt(tiny) + safe_exponent = np.log(tiny) / 2 + if min_exponent < safe_exponent: + warnings.warn( + f"reg={reg:g} is small enough that the convolution kernel " + f"underflows: its smallest exponent is {min_exponent:.0f} against a " + f"usable limit of {safe_exponent:.0f}. The result will be too " + "diffuse, and more iterations will not help. Use " + "method='sinkhorn_log' for this regularization.", + stacklevel=3, + ) + + def _get_convol_img_fn(nx, width, height, reg, type_as, log_domain=False): """Return the convolution operator for 2D images. @@ -34,6 +63,7 @@ def _get_convol_img_fn(nx, width, height, reg, type_as, log_domain=False): # If normal domain is selected, we can use M1 and M2 to compute the convolution if not log_domain: + _warn_if_kernel_underflows(nx, M1, M2, reg) K1, K2 = nx.exp(M1), nx.exp(M2) def convol_imgs(imgs): diff --git a/test/test_bregman.py b/test/test_bregman.py index 17b400306..8f8165ba4 100644 --- a/test/test_bregman.py +++ b/test/test_bregman.py @@ -1434,13 +1434,74 @@ def test_screenkhorn(nx): np.testing.assert_allclose(G_sink.sum(1), G_screen.sum(1), atol=1e-02) +def test_convolutional_barycenter_kernel_underflow_warns(): + """Small reg silently loses the transport (issue #458). + + exp(-(x-y)**2 / reg) underflows for distant pixels, mass can no longer + cross the image, and the barycenter collapses towards the arithmetic mean + of the inputs, which reads as an over-diffuse result rather than an error. + """ + rng = np.random.RandomState(0) + n, sigma, sep = 32, 0.06, 0.15 + t = np.linspace(0, 1, n) + X, Y = np.meshgrid(t, t, indexing="ij") + + def gauss(cx): + g = np.exp(-((X - cx) ** 2 + (Y - 0.5) ** 2) / (2 * sigma**2)) + return g / g.sum() + + A = np.stack([gauss(0.5 - sep), gauss(0.5 + sep)]) + + with pytest.warns(UserWarning, match="underflow"): + ot.bregman.convolutional_barycenter2d_debiased(A, 1e-04) + + with pytest.warns(UserWarning, match="underflow"): + ot.bregman.convolutional_barycenter2d(A, 1e-04) + + # a usable kernel must stay silent, and so must the log-domain solver + with warnings.catch_warnings(): + warnings.simplefilter("error", UserWarning) + ot.bregman.convolutional_barycenter2d_debiased(A, 1e-02) + ot.bregman.convolutional_barycenter2d_debiased(A, 1e-04, method="sinkhorn_log") + + +def test_convolutional_barycenter_debiased_preserves_width(): + """The debiased barycenter of two equal-width Gaussians keeps that width. + + Janati et al. 2020. The log-domain solver gets this right at every reg; + the default one only where its kernel has not underflowed. + """ + n, sigma, sep = 32, 0.06, 0.15 + t = np.linspace(0, 1, n) + X, Y = np.meshgrid(t, t, indexing="ij") + + def gauss(cx): + g = np.exp(-((X - cx) ** 2 + (Y - 0.5) ** 2) / (2 * sigma**2)) + return g / g.sum() + + A = np.stack([gauss(0.5 - sep), gauss(0.5 + sep)]) + + def width(img): + px = img.sum(axis=1) + mx = (px * t).sum() + return np.sqrt(((t - mx) ** 2 * px).sum()) + + bar = ot.bregman.convolutional_barycenter2d_debiased( + A, 1e-03, method="sinkhorn_log" + ) + np.testing.assert_allclose(width(bar), width(A[0]), rtol=0.05) + + def test_convolutional_barycenter_non_square(nx): # test for image with height not equal width A = np.ones((2, 2, 3)) / (2 * 3) A_nx = nx.from_numpy(A) - b_np = ot.bregman.convolutional_barycenter2d(A, 1e-03) - b = nx.to_numpy(ot.bregman.convolutional_barycenter2d(A_nx, 1e-03)) + # reg=1e-3 underflows the convolution kernel on a unit grid, which does not + # affect a uniform image but does emit a warning; 1e-2 exercises the same + # non-square code path with a usable kernel + b_np = ot.bregman.convolutional_barycenter2d(A, 1e-02) + b = nx.to_numpy(ot.bregman.convolutional_barycenter2d(A_nx, 1e-02)) np.testing.assert_allclose(np.ones((2, 3)) / (2 * 3), b, atol=1e-02) np.testing.assert_allclose(np.ones((2, 3)) / (2 * 3), b, atol=1e-02) From 2bea3bc2d10c4d49cca28343a68e44fc946955e4 Mon Sep 17 00:00:00 2001 From: Ahmed Eldeeb <62363199+deeb01@users.noreply.github.com> Date: Tue, 22 Sep 2026 16:52:52 -0700 Subject: [PATCH 2/3] Correct PR number in RELEASES.md --- RELEASES.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/RELEASES.md b/RELEASES.md index 43e8f2d7f..a73635696 100644 --- a/RELEASES.md +++ b/RELEASES.md @@ -15,7 +15,7 @@ #### Closed issues -- Warn when the convolution kernel in `ot.bregman.convolutional_barycenter2d` and `convolutional_barycenter2d_debiased` underflows at small `reg`. Past that point mass can no longer cross the image and the barycenter collapses towards the arithmetic mean of the inputs, which looks like an over-diffuse result rather than an error; `method='sinkhorn_log'` is exact in that regime (PR #867, Issue #458) +- Warn when the convolution kernel in `ot.bregman.convolutional_barycenter2d` and `convolutional_barycenter2d_debiased` underflows at small `reg`. Past that point mass can no longer cross the image and the barycenter collapses towards the arithmetic mean of the inputs, which looks like an over-diffuse result rather than an error; `method='sinkhorn_log'` is exact in that regime (PR #873, Issue #458) - Remove a leftover debug `print` from `ot.utils.projection_sparse_simplex` with `axis=1`, and make the `ot.datasets.make_gauss_hd` docstring a raw string so importing `ot` no longer emits a `SyntaxWarning` (PR #860) - Fix `ot.dist` ignoring the weights `w` for `metric="cityblock"`, which returned the unweighted distance although the weights are documented for this metric (PR #859) - Fix swapped arguments to `div_to_product` in `ot.gromov.fused_unbalanced_across_spaces_cost`: with `reg_type="independent"` (UCOOT) the entropic terms used the plan marginals as the reference measures and vice versa (PR #855, Issue #854) From 752c789c81cce1e9fd999a8fabcc463920867f5f Mon Sep 17 00:00:00 2001 From: Ahmed Eldeeb <62363199+deeb01@users.noreply.github.com> Date: Tue, 22 Sep 2026 16:58:54 -0700 Subject: [PATCH 3/3] Simplify the underflow check and point the warning at the caller Review of the previous commit: - it converted the whole width x width exponent matrix with nx.to_numpy purely to read a dtype, which copies the array and forces a device sync on GPU backends, once per barycenter call; - it compared two backend scalars with the builtin min, relying on __lt__ returning a Python bool across all five backends; - the reduction was unnecessary. The grid is always linspace(0, 1, n), so the most negative exponent is -1 / reg whatever the image size. The check now derives the exponent from reg alone and reads the dtype from a single-element array. stacklevel was 3, which attributed the warning to ot/bregman/_convolutional.py rather than to the caller. Measured against the real call stack: 4 lands on the user's call site. --- ot/bregman/_convolutional.py | 17 +++++++++++------ 1 file changed, 11 insertions(+), 6 deletions(-) diff --git a/ot/bregman/_convolutional.py b/ot/bregman/_convolutional.py index 3c4842c42..143758dcd 100644 --- a/ot/bregman/_convolutional.py +++ b/ot/bregman/_convolutional.py @@ -22,7 +22,7 @@ ) -def _warn_if_kernel_underflows(nx, M1, M2, reg): +def _warn_if_kernel_underflows(nx, reg, type_as, stacklevel): """Warn when exp(M) is about to lose the transport to underflow. The Sinkhorn iterations form products and ratios of kernel entries, so the @@ -30,12 +30,17 @@ def _warn_if_kernel_underflows(nx, M1, M2, reg): convolution can no longer move mass across the image and the barycenter degenerates towards the arithmetic mean of the inputs, which looks like an over-diffuse result rather than an error. + + The grid is always ``linspace(0, 1, n)``, so the most negative exponent is + ``-1 / reg`` whatever the image size, and no reduction over the kernel is + needed. """ - min_exponent = float(min(nx.min(M1), nx.min(M2))) try: - tiny = np.finfo(nx.to_numpy(M1).dtype).tiny + dtype = nx.to_numpy(nx.zeros((1,), type_as=type_as)).dtype + tiny = np.finfo(dtype).tiny except (TypeError, ValueError): # pragma: no cover - exotic dtypes return + min_exponent = -1.0 / reg # half the exponent range, i.e. the exponent of sqrt(tiny) safe_exponent = np.log(tiny) / 2 if min_exponent < safe_exponent: @@ -45,11 +50,11 @@ def _warn_if_kernel_underflows(nx, M1, M2, reg): f"usable limit of {safe_exponent:.0f}. The result will be too " "diffuse, and more iterations will not help. Use " "method='sinkhorn_log' for this regularization.", - stacklevel=3, + stacklevel=stacklevel, ) -def _get_convol_img_fn(nx, width, height, reg, type_as, log_domain=False): +def _get_convol_img_fn(nx, width, height, reg, type_as, log_domain=False, stacklevel=4): """Return the convolution operator for 2D images. The function constructed is equivalent to blurring on horizontal then vertical directions.""" @@ -63,7 +68,7 @@ def _get_convol_img_fn(nx, width, height, reg, type_as, log_domain=False): # If normal domain is selected, we can use M1 and M2 to compute the convolution if not log_domain: - _warn_if_kernel_underflows(nx, M1, M2, reg) + _warn_if_kernel_underflows(nx, reg, type_as, stacklevel=stacklevel + 1) K1, K2 = nx.exp(M1), nx.exp(M2) def convol_imgs(imgs):