From 79b49c84e8e3f5a1529a0614e9842027690c5e61 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Fri, 25 Sep 2026 22:20:31 +0100 Subject: [PATCH 1/3] doc: added deblurring tutorial and pywavelets doc dependency --- environment-dev-arm.yml | 1 + environment-dev.yml | 1 + pyproject.toml | 1 + tutorials/deblurring.py | 166 ++++++++++++++++++++++++++++++++++++++++ uv.lock | 7 ++ 5 files changed, 176 insertions(+) create mode 100644 tutorials/deblurring.py diff --git a/environment-dev-arm.yml b/environment-dev-arm.yml index 5198b2e..24a3b2e 100644 --- a/environment-dev-arm.yml +++ b/environment-dev-arm.yml @@ -9,6 +9,7 @@ dependencies: - scipy>=1.14.0 - pylops>=2.4.0 - scikit-image + - pywavelets - matplotlib - typing_extensions>=4.0.0 - ipython diff --git a/environment-dev.yml b/environment-dev.yml index 0eb2569..5848f4e 100644 --- a/environment-dev.yml +++ b/environment-dev.yml @@ -9,6 +9,7 @@ dependencies: - scipy>=1.14.0 - pylops>=2.4.0 - scikit-image + - pywavelets - matplotlib - typing_extensions>=4.0.0 - ipython diff --git a/pyproject.toml b/pyproject.toml index 9053a8f..f2bb575 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -60,6 +60,7 @@ doc = [ "scikit-image", "bm3d", "bm4d", + "pywavelets", "pooch", "sphinx", "sphinx-design", diff --git a/tutorials/deblurring.py b/tutorials/deblurring.py new file mode 100644 index 0000000..7e60089 --- /dev/null +++ b/tutorials/deblurring.py @@ -0,0 +1,166 @@ +r""" +Deblurring with multiple wavelet priors +======================================= +This tutorial considers the problem of deblurring an image that has been +convolved with a 2D Gaussian kernel and contaminated by random noise: + +.. math:: + \mathbf{y} = \mathbf{B} \mathbf{x} + \mathbf{n} + +where :math:`\mathbf{B}` is a 2D convolution operator and +:math:`\mathbf{n} \sim \mathcal{N}(0, \sigma^2)`. + +Natural images are rarely well represented by a single sparsifying transform: +piecewise-constant regions with sharp edges are compactly described by the +Haar wavelet, whilst smoother, textured regions are better captured by +higher-order wavelets such as the Daubechies-4 (db4). We will therefore +promote sparsity in both domains at the same time by solving: + +.. math:: + \argmin_\mathbf{x} \frac{1}{2}\|\mathbf{B}\mathbf{x}-\mathbf{y}\|_2^2 + + \lambda_1 \|\mathbf{W}_1 \mathbf{x}\|_1 + + \lambda_2 \|\mathbf{W}_2 \mathbf{x}\|_1 + +where :math:`\mathbf{W}_1` and :math:`\mathbf{W}_2` are the Haar and db4 +2D wavelet transforms, respectively. + +By stacking the two wavelet transforms into a single operator +:math:`\mathbf{K} = [\mathbf{W}_1^T, \mathbf{W}_2^T]^T`, the regularization +term becomes a separable function of :math:`\mathbf{K}\mathbf{x}`, whose +proximal operator is simply the stack of the proximal operators of the two +L1 norms. This problem can therefore be solved with the +:func:`pyproximal.optimization.primaldual.PrimalDual` solver +(Chambolle-Pock algorithm). + +""" + +import matplotlib.pyplot as plt +import numpy as np +import pylops +from pylops.utils.metrics import psnr + +import pyproximal + +plt.close("all") +np.random.seed(0) + +############################################################################### +# Let's start by loading an image, converting it to grayscale, and cropping +# it to a square of size power of two (so that the wavelet transforms do not +# require any padding). +im = plt.imread("../testdata/butterfly.png") +im = np.mean(im[..., :3], axis=-1) + +n = 256 +h, w = im.shape +im = im[(h - n) // 2 : (h - n) // 2 + n, (w - n) // 2 : (w - n) // 2 + n] +ny, nx = im.shape + +############################################################################### +# We now create a Gaussian blurring kernel and the associated +# :class:`pylops.signalprocessing.Convolve2D` operator. This is used to +# create the blurred data, to which we add white Gaussian noise with a +# signal-to-noise ratio of 30 dB. +nh = 11 +sigh = 2.0 +t = np.arange(-(nh // 2), nh // 2 + 1) +hy, hx = np.meshgrid(t, t, indexing="ij") +h = np.exp(-(hx**2 + hy**2) / (2 * sigh**2)) +h /= h.sum() + +Bop = pylops.signalprocessing.Convolve2D( + dims=(ny, nx), h=h, offset=(nh // 2, nh // 2), dtype="float64" +) + +yclean = Bop @ im.ravel() +snr_db = 30.0 +sigman = np.linalg.norm(yclean) / np.sqrt(yclean.size) * 10 ** (-snr_db / 20) +y = yclean + sigman * np.random.randn(yclean.size) + +fig, axs = plt.subplots(1, 3, figsize=(12, 4)) +axs[0].imshow(im, cmap="gray", vmin=0, vmax=1) +axs[0].set_title("True image") +axs[0].axis("off") +axs[1].imshow(y.reshape(ny, nx), cmap="gray", vmin=0, vmax=1) +axs[1].set_title(f"Blurred+noisy (PSNR={psnr(im, y.reshape(ny, nx), 1.0):.1f} dB)") +axs[1].axis("off") +axs[2].imshow(h, cmap="viridis") +axs[2].set_title(rf"Blurring kernel ($\sigma$={sigh})") +axs[2].axis("off") +plt.tight_layout() + +############################################################################### +# Next, we define the two wavelet transforms and stack them vertically into +# the operator :math:`\mathbf{K}`. +W1op = pylops.signalprocessing.DWT2D(dims=(ny, nx), wavelet="haar", level=3) +W2op = pylops.signalprocessing.DWT2D(dims=(ny, nx), wavelet="db4", level=3) +Kop = pylops.VStack([W1op, W2op]) +n1, n2 = W1op.shape[0], W2op.shape[0] + +############################################################################### +# We are now ready to define the proximal operators. The data term is a +# :class:`pyproximal.L2` norm with the blurring operator, whilst the +# regularization term is a :class:`pyproximal.VStack` of two +# :class:`pyproximal.L1` norms, each acting on the portion of the vector +# :math:`\mathbf{K}\mathbf{x}` associated with one of the two wavelet +# transforms. +lam1 = 1e-3 +lam2 = 1e-3 + +f = pyproximal.L2(Op=Bop, b=y) +g = pyproximal.VStack( + [pyproximal.L1(sigma=lam1), pyproximal.L1(sigma=lam2)], nn=[n1, n2] +) + +############################################################################### +# The Primal-Dual algorithm converges provided that +# :math:`\tau \mu \|\mathbf{K}\|_2^2 < 1`. Since each wavelet transform is +# (close to) orthogonal, we have :math:`\|\mathbf{K}\|_2^2 \le 2`, and choose +# the step lengths accordingly (with some safety margin). To appreciate the +# benefit of combining the two priors, we also solve the problem using a +# single wavelet transform at a time (for which +# :math:`\|\mathbf{W}_i\|_2^2 \le 1`). +tau = 0.4 +mu = 1.0 / (tau * 2.5) +mu1 = 1.0 / (tau * 1.25) +niter = 400 +x0 = np.zeros(ny * nx) + +# Haar only +xhaar = pyproximal.optimization.primaldual.PrimalDual( + f, pyproximal.L1(sigma=lam1), W1op, tau=tau, mu=mu1, x0=x0, niter=niter +) +xhaar = xhaar.reshape(ny, nx) + +# db4 only +xdb4 = pyproximal.optimization.primaldual.PrimalDual( + f, pyproximal.L1(sigma=lam2), W2op, tau=tau, mu=mu1, x0=x0, niter=niter +) +xdb4 = xdb4.reshape(ny, nx) + +# Haar + db4 +xdual = pyproximal.optimization.primaldual.PrimalDual( + f, g, Kop, tau=tau, mu=mu, x0=x0, niter=niter, show=True +) +xdual = xdual.reshape(ny, nx) + +############################################################################### +# Finally, we compare the three reconstructions. +fig, axs = plt.subplots(2, 2, figsize=(9, 9)) +for ax, img, title in zip( + axs.ravel(), + [y.reshape(ny, nx), xhaar, xdb4, xdual], + ["Blurred+noisy", "Haar prior", "db4 prior", "Haar + db4 priors"], + strict=True, +): + ax.imshow(img, cmap="gray", vmin=0, vmax=1) + ax.set_title(f"{title} (PSNR={psnr(im, img, 1.0):.1f} dB)") + ax.axis("off") +plt.tight_layout() + +############################################################################### +# Both single-prior reconstructions already remove most of the blur. The +# Haar prior favours piecewise-constant solutions and is therefore prone +# to blocky artefacts, whilst the db4 prior produces smoother images. By +# combining the two priors, we obtain a reconstruction that benefits from +# the strengths of both transforms. diff --git a/uv.lock b/uv.lock index b74d368..212e221 100644 --- a/uv.lock +++ b/uv.lock @@ -638,18 +638,23 @@ wheels = [ { url = "https://files.pythonhosted.org/packages/8d/3c/4716e285f8459a5707ef39d84e9a7547bf16f87c583ac4789c0763d6d986/cupy_cuda13x-14.2.0-cp311-cp311-manylinux2014_aarch64.whl", hash = "sha256:a22e2f17c97940cd566f351f186d8fe04285280d052ca1512ff223febbe9979c", size = 74074375, upload-time = "2026-08-20T02:41:04.327Z" }, { url = "https://files.pythonhosted.org/packages/ac/cb/41b6d8f1880e6f3128f8084b7622b1a3ca2d615ab77af9ef313b4aafed3d/cupy_cuda13x-14.2.0-cp311-cp311-manylinux2014_x86_64.whl", hash = "sha256:f94843e54500a9c8c169135f79572bfaba946efe4f7178e834c227280c1b0c12", size = 69969199, upload-time = "2026-08-20T02:41:07.842Z" }, { url = "https://files.pythonhosted.org/packages/f8/eb/25b8c1a415ade6479ed684c5cbff1e6f5002a92cbf26ef826df56a40cf69/cupy_cuda13x-14.2.0-cp311-cp311-win_amd64.whl", hash = "sha256:3893acd8d4dc78d113736045eaa306e496ee7bd8253bd959513ba28c91ed396d", size = 36082966, upload-time = "2026-08-20T02:41:10.723Z" }, + { url = "https://files.pythonhosted.org/packages/10/72/ad885cc21a21208b94044b34f50a840394b2e9dbb6111ed43959366aa6eb/cupy_cuda13x-14.2.0-cp311-cp311-win_arm64.whl", hash = "sha256:e0ad42a1ef35549282e2087a7c1e6ec8b7aa624a95e8fc0a60da777ca87e71de", size = 35136070, upload-time = "2026-09-01T06:41:38.284Z" }, { url = "https://files.pythonhosted.org/packages/dd/c8/68e07ff959f6a9f93169e37a64a4c24223487102b8625191e4a566cab465/cupy_cuda13x-14.2.0-cp312-cp312-manylinux2014_aarch64.whl", hash = "sha256:0a9b6538079c5151dcd99cf44b94cc3f1242315b7de0498b4619cef9b2a7c2c8", size = 73382974, upload-time = "2026-08-20T02:41:13.648Z" }, { url = "https://files.pythonhosted.org/packages/7e/bf/c4c25e97552320451afdcd459d4b6e0b715459594e998b23051c06db64c3/cupy_cuda13x-14.2.0-cp312-cp312-manylinux2014_x86_64.whl", hash = "sha256:9aafadbe29db3d4285026044aafeb28f2a3de4bd0df70718ec6892b7d27a4ded", size = 69439247, upload-time = "2026-08-20T02:41:17.375Z" }, { url = "https://files.pythonhosted.org/packages/5e/88/0a4d2d075485c4a82bfcd4725300ff87d17de36dce27b36f720d286f0f70/cupy_cuda13x-14.2.0-cp312-cp312-win_amd64.whl", hash = "sha256:9e8bc8a7693c5c212cd1450d728f47e1f05aa07ae003de30415893d0daf7d243", size = 35986921, upload-time = "2026-08-20T02:41:20.253Z" }, + { url = "https://files.pythonhosted.org/packages/a4/ab/8a756f553397b8509906d4a1d8083d92d3d21cbde115504e7a5db97de61f/cupy_cuda13x-14.2.0-cp312-cp312-win_arm64.whl", hash = "sha256:0dae6d788752c987c1cd87adb903640620820feb895339e8bc51a01fba0fb804", size = 35049803, upload-time = "2026-09-01T06:41:43.442Z" }, { url = "https://files.pythonhosted.org/packages/3c/02/3cb39774c71a16275c50b5d1e92ddcefa948ff1e69f464249b70a79b0abe/cupy_cuda13x-14.2.0-cp313-cp313-manylinux2014_aarch64.whl", hash = "sha256:449af7f608dc70ce64e8c45cf5a6ea0bbef51623e7c806b1a43445ae3c76e1e0", size = 72969322, upload-time = "2026-08-20T02:41:23.311Z" }, { url = "https://files.pythonhosted.org/packages/4e/27/216370ae8c39fa9835645be3c27513ab72f52cf4f51faf58d799fba981b6/cupy_cuda13x-14.2.0-cp313-cp313-manylinux2014_x86_64.whl", hash = "sha256:5b79a606f6639d9bed74ef7072f36756a95a804c86d7adb367b39307ac771600", size = 69048238, upload-time = "2026-08-20T02:41:26.831Z" }, { url = "https://files.pythonhosted.org/packages/5d/75/9df3c545baad4ff2e384db1ff6c8005fa4040492c4ee5672ca8b8711d29d/cupy_cuda13x-14.2.0-cp313-cp313-win_amd64.whl", hash = "sha256:f9a2a55abf889f6df68b5f62cd5f3ec4a2f305ac57ed09becd52e16201edcd86", size = 35965490, upload-time = "2026-08-20T02:41:29.943Z" }, + { url = "https://files.pythonhosted.org/packages/5c/3d/ee85ce39c2f32e8561c81c507df458a914b45d9a0eb1a54e5396c2cfd632/cupy_cuda13x-14.2.0-cp313-cp313-win_arm64.whl", hash = "sha256:5ddaf38291f60ec20e8c7889b28e339bf17ad1b4df768b8bb5e6c9f9e445efbb", size = 35029087, upload-time = "2026-09-01T06:41:48.593Z" }, { url = "https://files.pythonhosted.org/packages/e5/70/8c80c9e7010193afed8ccb7110d6173aa23bf3a6ce2fc2450244b0fa50a1/cupy_cuda13x-14.2.0-cp314-cp314-manylinux2014_aarch64.whl", hash = "sha256:7ac57c4e16c62f265b5ecc2b45d6bbafe9dc70a99b849ed6cfb8d16f32204589", size = 72840193, upload-time = "2026-08-20T02:41:33.146Z" }, { url = "https://files.pythonhosted.org/packages/c2/b8/4f4c4f34fc31ab8d136ed965a919505297974629d9852c50e15bfd616281/cupy_cuda13x-14.2.0-cp314-cp314-manylinux2014_x86_64.whl", hash = "sha256:ff0bdebd1b43c0c6db53095784c787c4e4eae671356cf521c4b7482ed78a1a7e", size = 68422828, upload-time = "2026-08-20T02:41:36.527Z" }, { url = "https://files.pythonhosted.org/packages/7d/5c/58941ece5ed0b2cc8304f8ba81b5e8334759d409306475d53c52ac99809b/cupy_cuda13x-14.2.0-cp314-cp314-win_amd64.whl", hash = "sha256:b99a8bf9d5391954c0c08675a8cb39219923fe2ee249faa3546ebbf75783a360", size = 36110882, upload-time = "2026-08-20T02:41:39.837Z" }, + { url = "https://files.pythonhosted.org/packages/e2/c5/cc7d1a9b764b3f828e09b48a73d0fdeaad15ecd2c7e50507f6d5fcf5fd1a/cupy_cuda13x-14.2.0-cp314-cp314-win_arm64.whl", hash = "sha256:c61c59ba665948c0a855926997ee96026420b628816b2e406998ec797caddb22", size = 35184806, upload-time = "2026-09-01T06:41:58.317Z" }, { url = "https://files.pythonhosted.org/packages/6f/54/a582ac70e9e75d3fad7ae073051bcd50dd128a1f66dd275291dc5f48d997/cupy_cuda13x-14.2.0-cp314-cp314t-manylinux2014_aarch64.whl", hash = "sha256:e8bb453551deced08ec5c4576da61ebf1e2b4a45aeff48e09f47ec332f7cad3d", size = 73336756, upload-time = "2026-08-20T02:41:42.718Z" }, { url = "https://files.pythonhosted.org/packages/f9/aa/222151ba5061d7a01efe4009efad9f9a7164e15bb2a494d9aa86efc97d4f/cupy_cuda13x-14.2.0-cp314-cp314t-manylinux2014_x86_64.whl", hash = "sha256:313500c5d415f65d388e84a9d2eb7f07f97c23c85dc5785f986e75d46d9c9693", size = 68792663, upload-time = "2026-08-20T02:41:46.079Z" }, { url = "https://files.pythonhosted.org/packages/d4/20/13dd37dbeaec879c8348fde2bfca67d96daa56a6e63e4a5083ff61660e53/cupy_cuda13x-14.2.0-cp314-cp314t-win_amd64.whl", hash = "sha256:c517fc55502a7e3d1a8b9310abad95cec61dded5a38d328fe0220bcf7b036268", size = 37004845, upload-time = "2026-08-20T02:41:48.963Z" }, + { url = "https://files.pythonhosted.org/packages/d1/5b/b1367f79abd90a259dbb7f0603dd162c03c3a76ec1d02c2527b0c79d46ba/cupy_cuda13x-14.2.0-cp314-cp314t-win_arm64.whl", hash = "sha256:83885db963eb0b95a0757b4907a37115a1d2aebfa9cf43b99ebac8db35713d37", size = 35485280, upload-time = "2026-09-01T06:41:53.618Z" }, ] [[package]] @@ -1927,6 +1932,7 @@ doc = [ { name = "nbsphinx" }, { name = "numpydoc" }, { name = "pooch" }, + { name = "pywavelets" }, { name = "scikit-image" }, { name = "shibuya" }, { name = "sphinx", version = "9.0.4", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version < '3.12' or (extra == 'extra-10-pyproximal-gpu-cu12' and extra == 'extra-10-pyproximal-gpu-cu13')" }, @@ -1969,6 +1975,7 @@ doc = [ { name = "nbsphinx" }, { name = "numpydoc" }, { name = "pooch" }, + { name = "pywavelets" }, { name = "scikit-image" }, { name = "shibuya" }, { name = "sphinx" }, From 1d0249f35acfa8685c032d3ffb7fa581674d4b04 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Fri, 25 Sep 2026 22:20:57 +0100 Subject: [PATCH 2/3] feat: added support for vector sigma in L1 --- pyproximal/proximal/L1.py | 27 ++++++++------ pytests/test_norms.py | 77 +++++++++++++++++++++++++++++---------- 2 files changed, 73 insertions(+), 31 deletions(-) diff --git a/pyproximal/proximal/L1.py b/pyproximal/proximal/L1.py index 42129ae..6427d55 100644 --- a/pyproximal/proximal/L1.py +++ b/pyproximal/proximal/L1.py @@ -54,7 +54,8 @@ class L1(ProxOperator): Proximal operator of the :math:`\ell_1` norm: - :math:`\sigma\|\mathbf{x} - \mathbf{g}\|_1 = \sigma \sum_i |x_i - g_i|` for - 1-dimensional arrays; + 1-dimensional arrays (or :math:`\sum_i \sigma_i |x_i - g_i|` when + :math:`\sigma` is a vector of the same size as :math:`\mathbf{x}`); - :math:`\sum_j \sigma_j \|\mathbf{X}_j - \mathbf{g}\|_1 = \sum_j \sigma_j \sum_i |x_{i,j} - g_i|` for 2-dimensional arrays. @@ -62,10 +63,10 @@ class L1(ProxOperator): ---------- sigma : :obj:`float` or :obj:`numpy.ndarray` or :obj:`func`, optional Multiplicative coefficient of L1 norm. This can be a constant number, a list - of values (for 2-dimensional inputs, acting on the second dimension) or - a function that is called passing a counter which keeps track of how many - times the ``prox`` method has been invoked before and returns a scalar (or a list of) - ``sigma`` to be used. + of values (for 1-dimensional inputs, acting element-wise; for 2-dimensional + inputs, acting on the second dimension) or a function that is called passing + a counter which keeps track of how many times the ``prox`` method has been + invoked before and returns a scalar (or a list of) ``sigma`` to be used. g : :obj:`numpy.ndarray`, optional Vector to be subtracted @@ -98,6 +99,12 @@ class L1(ProxOperator): \sigma, & x_i > \sigma\\ \end{cases} + When :math:`\mathbf{g} \neq \mathbf{0}`, the convex conjugate becomes + :math:`f^*(\mathbf{x}) = \mathbf{x}^T\mathbf{g} + + \mathcal{I}_{\|\cdot\|_{\infty} <=\sigma}(\mathbf{x})`, and the dual + operator is obtained by applying the same projection to + :math:`\mathbf{x} - \tau \mathbf{g}`. + .. [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. @@ -122,8 +129,8 @@ def __init__( def __call__(self, x: NDArray) -> float: sigma = _current_sigma(self.sigma, self.count) if self.g is None: - return float(np.sum(sigma * np.sum(np.abs(x), axis=0))) - return float(np.sum(sigma * np.sum(np.abs(x - self.g), axis=0))) + return float(np.sum(sigma * np.abs(x))) + return float(np.sum(sigma * np.abs(x - self.g))) def _increment_count(func: Callable[..., Any]) -> Callable[..., Any]: """Increment counter""" @@ -147,11 +154,7 @@ def prox(self, x: NDArray, tau: float) -> NDArray: @_check_tau def proxdual(self, x: NDArray, tau: float) -> NDArray: - if not isinstance(self.gdual, np.ndarray): - x = self.box(x) - else: - x = self._proxdual_moreau(x, tau) - return x + return self.box(x - tau * self.gdual) class L1Ball(ProxOperator): diff --git a/pytests/test_norms.py b/pytests/test_norms.py index 65c2ebb..1a3ace1 100644 --- a/pytests/test_norms.py +++ b/pytests/test_norms.py @@ -314,23 +314,48 @@ def test_L2Convolve(par): @pytest.mark.parametrize("par", [(par1), (par2)]) -def test_L1(par): - """L1 norm and grad and proximal/dual proximal""" +def test_L1_scalar(par): + """L1 norm (with and without g) and grad and proximal/dual proximal + with sigma as scalar""" np.random.seed(10) - l1 = L1(sigma=par["sigma"]) + g = np.random.normal(0.0, 1.0, par["nx"]).astype(par["dtype"]) + x = np.random.normal(0.0, 1.0, par["nx"]).astype(par["dtype"]) + tau = 2.0 + for gg in [None, g]: + l1 = L1(sigma=par["sigma"], g=gg) + xg = x if gg is None else x - gg + + # norm + assert l1(x) == pytest.approx(par["sigma"] * np.sum(np.abs(xg))) + + # grad (note that since this norm is non-smooth, the + # gradient method is not implemented; as such the gradient + # of the Moreau envelope is called instead + _ = l1.grad(x) + + # prox / dualprox + assert moreau(l1, x, tau) + + +@pytest.mark.parametrize("par", [(par1), (par2)]) +def test_L1_vector(par): + """L1 norm and proximal/dual proximal with sigma as vector""" + np.random.seed(10) + + sigma = np.random.uniform(0.1, 2.0, par["nx"]).astype(par["dtype"]) + l1 = L1(sigma=sigma) # norm x = np.random.normal(0.0, 1.0, par["nx"]).astype(par["dtype"]) - assert l1(x) == par["sigma"] * np.sum(np.abs(x)) + assert l1(x) == pytest.approx(np.sum(sigma * np.abs(x))) - # grad (note that since this norm is non-smooth, the - # gradient method is not implemented; as such the gradient - # of the Moreau envelope is called instead - _ = l1.grad(x) + # prox + tau = 2.0 + xprox = np.sign(x) * np.maximum(np.abs(x) - tau * sigma, 0.0) + assert_array_almost_equal(l1.prox(x, tau), xprox) # prox / dualprox - tau = 2.0 assert moreau(l1, x, tau) @@ -351,20 +376,34 @@ def test_L1_func(par): @pytest.mark.parametrize("par", [(par1), (par2)]) -def test_L1_diff(par): - """L1 norm of difference and proximal/dual proximal""" +def test_L1_2d(par): + """L1 norm and proximal/dual proximal of 2-dimensional input with + sigma as vector acting on the second dimension""" np.random.seed(10) - g = np.random.normal(0.0, 1.0, par["nx"]).astype(par["dtype"]) - l1 = L1(sigma=par["sigma"], g=g) - - # norm - x = np.random.normal(0.0, 1.0, par["nx"]).astype(par["dtype"]) - assert l1(x) == par["sigma"] * np.sum(np.abs(x - g)) + ncol = 3 + sigma = np.random.uniform(0.1, 2.0, ncol).astype(par["dtype"]) + g = np.random.normal(0.0, 1.0, (par["nx"], 1)).astype(par["dtype"]) - # prox / dualprox + x = np.random.normal(0.0, 1.0, (par["nx"], ncol)).astype(par["dtype"]) tau = 2.0 - assert moreau(l1, x, tau) + for gg in [None, g]: + l1 = L1(sigma=sigma, g=gg) + xg = x if gg is None else x - gg + + # norm + norm = sum(sigma[j] * np.sum(np.abs(xg[:, j])) for j in range(ncol)) + assert l1(x) == pytest.approx(norm) + + # prox + xprox = np.sign(xg) * np.maximum(np.abs(xg) - tau * sigma, 0.0) + if gg is not None: + xprox += gg + assert_array_almost_equal(l1.prox(x, tau), xprox) + + # dualprox (conjugate of sigma||x - g||_1 is + i_{|x|<=sigma}) + xdual = x if gg is None else x - tau * gg + assert_array_almost_equal(l1.proxdual(x, tau), np.clip(xdual, -sigma, sigma)) @pytest.mark.parametrize("par", [(par1), (par2)]) From c9322debbc76f18ae57330180e4ff56a69ca93da Mon Sep 17 00:00:00 2001 From: mrava87 Date: Fri, 25 Sep 2026 22:34:54 +0100 Subject: [PATCH 3/3] minor: updated steps in deblurring tutorial --- tutorials/deblurring.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/tutorials/deblurring.py b/tutorials/deblurring.py index 7e60089..ac97281 100644 --- a/tutorials/deblurring.py +++ b/tutorials/deblurring.py @@ -120,21 +120,25 @@ # benefit of combining the two priors, we also solve the problem using a # single wavelet transform at a time (for which # :math:`\|\mathbf{W}_i\|_2^2 \le 1`). + +Lw = 1.0 # maxeig(W_i^H W_i) +L = 2.0 # maxeig(K^H K) + tau = 0.4 -mu = 1.0 / (tau * 2.5) -mu1 = 1.0 / (tau * 1.25) +muw = 0.8 / (tau * Lw) +mu = 0.8 / (tau * L) niter = 400 x0 = np.zeros(ny * nx) # Haar only xhaar = pyproximal.optimization.primaldual.PrimalDual( - f, pyproximal.L1(sigma=lam1), W1op, tau=tau, mu=mu1, x0=x0, niter=niter + f, pyproximal.L1(sigma=lam1), W1op, tau=tau, mu=muw, x0=x0, niter=niter ) xhaar = xhaar.reshape(ny, nx) # db4 only xdb4 = pyproximal.optimization.primaldual.PrimalDual( - f, pyproximal.L1(sigma=lam2), W2op, tau=tau, mu=mu1, x0=x0, niter=niter + f, pyproximal.L1(sigma=lam2), W2op, tau=tau, mu=muw, x0=x0, niter=niter ) xdb4 = xdb4.reshape(ny, nx)