diff --git a/examples/plot_norms.py b/examples/plot_norms.py index bd80d371..0c5f59fe 100755 --- a/examples/plot_norms.py +++ b/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) diff --git a/pyproximal/proximal/L1.py b/pyproximal/proximal/L1.py index 6427d555..e21818cc 100644 --- a/pyproximal/proximal/L1.py +++ b/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): diff --git a/pytests/test_l1_dual.py b/pytests/test_l1_dual.py new file mode 100644 index 00000000..4518bcc8 --- /dev/null +++ b/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)