Preserve complex magnitudes and current weights in the L1 dual proximal

examples/plot_norms.py

@@ -126,6 +126,13 @@
 plt.legend()
 plt.tight_layout()
 
+###############################################################################
+# For complex data, the dual L1 proximal projects onto magnitude disks and
+# preserves phase, rather than clipping real and imaginary components.
+complex_l1 = pyproximal.L1(sigma=1.0)
+complex_x = np.array([0j, 2j, 3 + 4j])
+print("Complex L1 dual projection: ", complex_l1.proxdual(complex_x, 0.5))
+
 ###############################################################################
 # We consider now the TV norm.
 tv = pyproximal.TV(dims=(nx,), sigma=1.0)

pyproximal/proximal/L1.py

@@ -105,6 +105,10 @@ class L1(ProxOperator):
     operator is obtained by applying the same projection to
     :math:`\mathbf{x} - \tau \mathbf{g}`.
 
+    For complex inputs, this is a radial projection onto the disks
+    :math:`|x_i| \leq \sigma_i`, preserving the phase of each shifted input.
+    Callable weights use the current primal-call counter in both operators.
+
     .. [1] Chambolle, and A., Pock, "A first-order primal-dual algorithm for
         convex problems with applications to imaging", Journal of Mathematical
         Imaging and Vision, 40, 8pp. 120–145. 2011.
@@ -154,7 +158,14 @@ def prox(self, x: NDArray, tau: float) -> NDArray:
 
     @_check_tau
     def proxdual(self, x: NDArray, tau: float) -> NDArray:
-        return self.box(x - tau * self.gdual)
+        sigma = _current_sigma(self.sigma, self.count)
+        self.box = BoxProj(-sigma, sigma)
+        x = x - tau * self.gdual
+        if np.iscomplexobj(x):
+            denominator = np.maximum(np.abs(x), sigma)
+            denominator = np.where(denominator == 0, 1, denominator)
+            return x * (sigma / denominator)
+        return self.box(x)
 
 
 class L1Ball(ProxOperator):

pytests/test_l1_dual.py

@@ -0,0 +1,39 @@
+import numpy as np
+import pytest
+
+from pyproximal import L1
+from pyproximal.utils import moreau
+
+
+@pytest.mark.parametrize("dtype", [np.complex64, np.complex128])
+@pytest.mark.parametrize("weighted", [False, True])
+@pytest.mark.parametrize("shifted", [False, True])
+def test_l1_complex_dual_projection(dtype, weighted, shifted):
+    """Complex dual projection preserves phase and satisfies Moreau's identity."""
+    x = np.array([0, 2j, 3 + 4j, -4 + 3j, 0.1 - 0.1j], dtype=dtype)
+    sigma = np.array([0, 0.2, 0.7, 2, 0.1]) if weighted else 0.7
+    g = np.array([0, 0.1j, 0.2, -0.3j, 0.1], dtype=dtype) if shifted else None
+    tau = 0.3
+    z = x if g is None else x - tau * g
+    expected = np.array(
+        [
+            value if abs(value) <= radius else radius * value / abs(value)
+            for value, radius in zip(z, np.broadcast_to(sigma, x.shape), strict=True)
+        ]
+    )
+    operator = L1(sigma=sigma, g=g)
+    np.testing.assert_allclose(operator.proxdual(x, tau), expected, atol=2e-7)
+    assert moreau(operator, x, tau, tol=2e-6)
+
+
+@pytest.mark.parametrize("complex_input", [False, True])
+def test_l1_dual_current_weights(complex_input):
+    """Dual bounds use the same callable weight as the latest primal call."""
+    x = np.array([3.0, -4.0])
+    if complex_input:
+        x = x.astype(complex) * 1j
+    operator = L1(sigma=lambda count: 1.0 + count)
+    operator.prox(x, 0.5)
+    np.testing.assert_allclose(operator.proxdual(x, 0.5), 2 * x / np.abs(x))
+    assert operator.count == 1
+    assert moreau(operator, x, 0.5)