Fix: preserve feasible inputs in L1 ball projection

pyproximal/projection/L1.py

@@ -30,7 +30,12 @@ class L1BallProj:
 
     .. math::
 
-        P_{L1_{r}} (\mathbf{x}) = sign(\mathbf{x}) P_{\operatorname{Simplex}(r)}(\mathbf{x})
+        P_{L1_{r}} (\mathbf{x}) =
+        \begin{cases}
+        \mathbf{x}, & \|\mathbf{x}\|_1 \leq r \\
+        sign(\mathbf{x}) P_{\operatorname{Simplex}(r)}(|\mathbf{x}|),
+        & \text{otherwise}
+        \end{cases}
 
     Note that this is the proximal operator of the corresponding
     indicator function :math:`\mathcal{I}_{L1_{r}}`.
@@ -49,6 +54,8 @@ def __init__(
         self.simplex = SimplexProj(n, radius, maxiter, xtol)
 
     def __call__(self, x: NDArray) -> NDArray:
+        if np.sum(np.abs(x)) <= self.radius:
+            return x
         if np.iscomplexobj(x):
             return np.exp(1j * np.angle(x)) * self.simplex(np.abs(x))
         else:

pytests/test_projection.py

@@ -1,8 +1,9 @@
 import numpy as np
 import pytest
-from numpy.testing import assert_array_almost_equal
+from numpy.testing import assert_array_almost_equal, assert_array_equal
 from pylops.basicoperators import Identity
 
+from pyproximal.projection import L1BallProj
 from pyproximal.proximal import (
     AffineSet,
     Box,
@@ -157,6 +158,38 @@ def test_L1Ball(par):
     assert moreau(l1, x, tau)
 
 
+@pytest.mark.parametrize("dtype", ["float32", "float64", "complex64", "complex128"])
+@pytest.mark.parametrize("values", [[0.125, -0.25], [1.0, 0.0], [0.0, 0.0]])
+def test_L1Ball_feasible_inputs(dtype, values):
+    """The projection fixes every point already in the ball."""
+    x = np.array(values, dtype=dtype)
+    original = x.copy()
+    projection = L1BallProj(x.size, radius=1.0)
+    ball = L1Ball(x.size, radius=1.0)
+
+    assert_array_equal(projection(x), original)
+    assert_array_equal(ball.prox(x, 0.5), original)
+    assert_array_equal(ball.prox(x, 2.0), original)
+    assert_array_equal(x, original)
+
+
+@pytest.mark.parametrize("dtype", ["complex64", "complex128"])
+def test_L1Ball_complex_interior(dtype):
+    x = np.array([0.125 + 0.125j, 0.25 - 0.125j], dtype=dtype)
+    ball = L1Ball(x.size, radius=1.0)
+    assert_array_equal(ball.prox(x, 1.0), x)
+
+
+@pytest.mark.parametrize("dtype", ["float32", "float64", "complex64", "complex128"])
+@pytest.mark.parametrize("tau", [0.5, 2.0])
+def test_L1Ball_dual_inside_scaled_ball(dtype, tau):
+    """The conjugate prox is zero when x/tau lies inside the ball."""
+    x = np.array([0.125, -0.25], dtype=dtype)
+    ball = L1Ball(x.size, radius=1.0)
+    assert_array_equal(ball.proxdual(x, tau), np.zeros_like(x))
+    assert moreau(ball, x, tau)
+
+
 @pytest.mark.parametrize("par", [(par1), (par2)])
 def test_NuclBall(par):
     """Nuclear Ball projection and proximal/dual proximal of related indicator"""