# Multiresolution Adaptive Block-Coordinate Forward-Backward for Image Reconstruction

Edgar Desainte-Maréville, Marion Foare, Paulo Gonçalves, Nelly Pustelnik, Elisa Riccietti

► To cite this version:

Edgar Desainte-Maréville, Marion Foare, Paulo Gonçalves, Nelly Pustelnik, Elisa Riccietti. Multiresolution Adaptive Block-Coordinate Forward-Backward for Image Reconstruction. 2026. hal-05514715

HAL Id: hal-05514715

https://hal.science/hal-05514715v1

Preprint submitted on 25 Feb 2026

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Multiresolution Adaptive Block-Coordinate
Forward–Backward for Image Reconstruction

Edgar Desainte-Maréville¹, Marion Foare¹,³, Paulo Gonçalves¹, Nelly Pustelnik², Elisa Riccietti¹

¹ENS de Lyon, CNRS, Inria, LIP, UMR 5668, Université Claude Bernard Lyon 1, Lyon, France

²CNRS, ENS de Lyon, Laboratoire de Physique, UMR 5672, Lyon, France

³CPE Lyon, Villeurbanne, France

{edgar.desainte-mareville, marion.foare, paulo.goncalves, nelly.pustelnik, elisa.riccietti}@ens-lyon.fr

**Abstract**—Classical first-order optimization methods for imaging inverse problems scale poorly with image resolution. Wavelet-based multilevel strategies can accelerate convergence under strong blur, but their fixed coarse-to-fine schedules lose effectiveness in moderate-blur or noise-dominated regimes.

In this work, we propose an adaptive multiresolution block-coordinate Forward–Backward algorithm for image restoration. Multiresolution block selection is driven by the local magnitude of the proximal update via a stochastic non-smooth Gauss–Southwell rule applied to the wavelet decomposition of the image. This adaptive selection strategy dynamically balances updates across scales, emphasizing coarse or fine blocks according to the degradation regime. As a result, the proposed method automatically adapts to varying blur and noise levels without relying on a predefined hierarchical update scheme.

**Index Terms**—Inverse problems, image restoration, mult-scale optimization, block-coordinate Forward–Backward, Gauss–Southwell rule

## I. INTRODUCTION

Image restoration aims at recovering a ground truth image from a degraded observation. Such problems are commonly formulated within a variational framework, leading to large-scale convex optimization problems of the form

$$ \min_{x \in \mathbb{R}^n} f(x) + h(x), \quad (1) $$

where $f$ is a smooth data-fidelity term and $h$ is a possibly non-smooth regularization term. First-order methods based on proximal splitting, such as the Forward–Backward (FB) scheme [1], [2] and its accelerated variants [3], are widely used due to their simplicity and robustness. These approaches remain central in model-based neural network methods, such as Plug-and-Play (PnP) frameworks, where a trained denoiser replaces the proximal operator. Among these, Forward–Backward PnP (FB-PnP) [4], [5], [6] stands out as one of the most widely used hybridizations of neural networks and variational techniques.

Despite their strong theoretical guarantees, these methods often exhibit slow convergence, which becomes limiting when applied to high-resolution images. Improving the scalability of optimization algorithms with respect to image resolution is therefore a central challenge in imaging inverse problems.

A standard strategy to address large-scale optimization problems consists in decomposing the variables into smaller subsets, referred to as blocks, and performing block-wise updates. Block-coordinate descent (BCD) methods proceed

by iteratively updating one or several blocks of variables while keeping the remaining components fixed [7].

When the objective function takes the form (1) the resulting family of methods is referred to as block-coordinate Forward–Backward (BC-FB), and each method is determined by the block selection strategy. A simple choice consists in updating all the blocks at each iteration, which recovers the standard Forward–Backward algorithm. Alternatively, stochastic BC-FB methods rely on random block activations, typically drawn independently at each iteration according to a fixed probability distribution, and enjoy theoretical convergence guarantees [8], [9]. However, such uniform sampling strategies do not exploit the structure of the problem and may lead to minor improvements. More recently, the classical deterministic Gauss–Southwell rule [10], which selects at each iteration the block associated with the largest partial gradient norm, has been extended to a non-smooth context by selecting the block associated with the largest proximal update [11].

While block-coordinate strategies are conceptually appealing, a direct extension of these updates to imaging inverse problems would lead to patch-induced block decompositions, which can complicate the formulation of the subproblems, especially in the presence of global regularization terms, due to boundary effects and interactions between blocks [12].

A more natural way to build a block decomposition for images is to exploit their intrinsic multiscale structure by leveraging wavelet-based decompositions [13]. Specifically, recent works in multilevel optimization [14], [15] proposed a Multilevel Forward–Backward (MLFB) method that chooses blocks corresponding to approximation and detail coefficients across scales, and performs BC-FB updates according to fixed coarse-to-fine schedules, where coarse scales correspond to wavelet approximation coefficients. This approach can be interpreted as a BC-FB method with a predetermined hierarchical update order, and has been shown to accelerate convergence, especially in regimes dominated by strong blur and low noise.

However, when the blur is moderate or when the noise becomes dominant, these fixed coarse-to-fine update rules tend to lose effectiveness. This observation suggests that no single fixed ordering is optimal across different degradation regimes, and motivates the design of adaptive block-selection strategies capable of dynamically identifying the components that will
contribute the most to improving the reconstruction during the optimization process.

**Contribution** – In this work, we propose a Multiresolution Adaptive Gauss-Southwell Inspired Coordinate Forward-Backward (MAGIC-FB) algorithm for image restoration. Block selection is inspired by a proximal Gauss-Southwell rule applied across resolution levels rather than spatial patches, which, to the best of our knowledge, has never been explored in imaging. Moreover, to allow for parallel block updates, we propose a stochastic variant of the Gauss-Southwell selection rule, with probability distributions driven by the local magnitude of the proximal update. Particular attention is paid to the computation of partial gradients, which leverages the linearity of the problem to enable cheap updates.

This adaptive strategy dynamically balances updates across scales, selecting coarse or fine components depending on the degradation type. As a result, the proposed method automatically adapts to varying blur and noise levels without relying on a predefined classical multilevel update schedule, thus making the proposed method a good choice in any degradation scenario.

We demonstrate indeed through numerical experiments that, across a wide range of blur and noise regimes, the proposed strategy consistently matches or outperforms the best methods, whose performance are on the contrary heavily affected by the degradation type.

**Notations** – The Euclidean norm is denoted by $|| \cdot ||$. For a convex, proper, and lower semicontinuous function $g$, the proximal operator of $g$ is defined as $\text{prox}_g(z) := \text{argmin}_{u \in \mathbb{R}^n} \{g(u) + \frac{1}{2} ||u - z||^2\}$. For a block-structured vector $w$, $w_i$ denotes its $i$-th block, while superscripts denote iterations indices. For a function $f: \mathbb{R}^n \to \mathbb{R}$, $\nabla f$ denotes its gradient when it exists and $\nabla_{w_i} f(w)$ denotes the partial gradient with respect to the block $w_i$.

## II. WAVELET-BASED BLOCK-COORDINATE FORWARD-BACKWARD

### A. Wavelet-based imaging inverse problem

We consider an imaging inverse problem where an unknown image $\bar{x} \in \mathbb{R}^n$ is observed through a linear degradation operator $A \in \mathbb{R}^{m \times n}$ and corrupted by additive Gaussian noise:

$$y = A\bar{x} + \eta, \quad (2)$$

where $\eta \sim \mathcal{N}(0, \sigma^2 \text{Id})$. In this work, we focus on variational approaches formulated as in (1), where the penalization is expressed in the wavelet domain. The wavelet-based formulation was extremely popular a few decades ago, and receives nowadays a renewed interest for its capability of improving the reconstruction performance when combined with neural networks, as recently highlighted in [16], [17], where the multiscale structure plays a crucial role in achieving high-quality restoration.

We thus focus on a minimization problem expressed as:

$$\hat{w} \in \underset{w \in \mathbb{R}^n}{\operatorname{Argmin}} \varphi(w) := \underbrace{\frac{1}{2} \|AW^\top w - y\|^2}_{=} + g(w), \quad (3)$$

where $W \in \mathbb{R}^{n \times n}$ is an orthonormal wavelet transform with $J$ decomposition levels and the reconstructed image is $\hat{x} = W^\top \hat{w}$. The function $w \mapsto \frac{1}{2} \|AW^\top w - y\|^2$ is $L$-smooth with $L = \|AW^\top\|^2 = \|A\|^2$.

The regularization term $g: \mathbb{R}^n \to \mathbb{R}$ is assumed to be convex, proper, and lower semicontinuous, and is typically non-smooth. It is well established in the literature [18] that, for an orthonormal wavelet transform, problem (3) is equivalent to (1) with $h = g(W \cdot)$.

### B. Wavelet representation and separable regularization

Wavelet transforms provide a multiresolution representation of an image, which is decomposed into low-frequency and high-frequency components.

To avoid distinguishing between approximation and detail coefficients and to emphasize the block structure induced by the wavelet transform, we introduce the unified notation

$$Wx = [a_J^\top, d_J^\top, \dots, d_1^\top]^\top = [w_0^\top, w_1^\top, \dots, w_J^\top]^\top = w$$

where $a_J = w_0 \in \mathbb{R}^{n_0}$ denotes the approximation coefficients at the coarsest scale, and $d_J = w_{J-j+1} \in \mathbb{R}^{n_{J-j+1}}$ contains the detail coefficients at scale $j$, formed by concatenating the horizontal, vertical, and diagonal sub-bands.

We assume that the regularization term is separable across wavelet blocks, i.e., that its contribution on each block is given by proper, convex, lower semicontinuous functions $g_i: \mathbb{R}^{n_i} \to (-\infty, +\infty]$ for all $i \in \{0, \dots, J\}$ such that

$$(\forall w \in \mathbb{R}^n) \quad g(w) = \sum_{i=0}^{J} g_i(w_i). \quad (4)$$

### C. General algorithm and state-of-the-art update rules

The BC-FB method designed to solve (3) exploiting a multiresolution architecture [15] is detailed in Algorithm 1.

<caption>**Algorithm 1:** BC-FB for (3)</caption>

Input: $w^{[0]} \in \text{dom} \varphi$, step $\gamma \in (0, \frac{2}{L})$, $\varepsilon^{[k]} \in \{0, 1\}^{J+1}$
for $k = 0, 1, \dots$ do
| for $i = 0, \dots, J$ do
| | if $\varepsilon_i^{[k]} = 1$ then
| | | $w_i^{[k+1]} = \text{prox}_{\gamma g_i}(w_i^{[k]} - \gamma \nabla_i f(w_i^{[k]}))$;
| | else
| | | $w_i^{[k+1]} = w_i^{[k]}$;

The sequence $(\varepsilon_i^{[k]})_{k \in \mathbb{N}}$ encodes the update rule: for each $k$, block $i$ is updated if $\varepsilon_i^{[k]} = 1$, and left unchanged otherwise. The update rules of the state-of-the-art are the following:

* Forward-Backward: All blocks are updated at each iteration, i.e., $\varepsilon_i^{[k]} = 1, \forall k \in \mathbb{N}$ and $\forall 0 \le i \le J$.
* Uniform stochastic sampling: At each iteration k and for each block i, the update variable $\epsilon_i^{[k]}$ is drawn independently according to a Bernoulli distribution with a fixed parameter p, i.e., $\epsilon_i^{[k]} \sim B(p)$.

* Multilevel Forward-Backward: The MLFB strategy follows a deterministic coarse-to-fine update pattern with $\epsilon_i^{[k]} = 1 \quad \forall i < k \mod (J + 1)$.

The next section describes our proposed adaptive block selection strategy.

### III. MAGIC-FB

#### A. Proposed multiresolution Gauss-Southwell update

We now introduce MAGIC-FB, an adaptive multiresolution block-coordinate FB algorithm based on a new stochastic Gauss-Southwell selection rule. The method fits within the general framework in Algorithm 1 and corresponds to a specific choice of the block-update sequence $(\epsilon_i^{[k]})_{k \in \mathbb{N}}$.

Given the current iterate $w^{[k]} = (w_0^{[k]}, \dots, w_J^{[k]})$ and a stepsize $\gamma > 0$, we define for each block $0 \le i \le J$ the block-wise proximal update

$$ \Delta_i^{[k]} := w_i^{[k]} - \text{prox}_{\gamma g_i} \left( w_i^{[k]} - \gamma \nabla_{w_i} f \left( w_i^{[k]} \right) \right), $$

which quantifies the variation induced by a proximal gradient step restricted to block i. At iteration k, the deterministic Gauss-Southwell rule would select the block index $i^{[k]} \in \text{argmax}_{0 \le i \le J} \|\Delta_i^{[k]}\|.$

To leverage the flexibility of randomness and to allow for parallel updates, we introduce a stochastic variant by defining activation probabilities $p_i^{[k]} = \frac{\|\Delta_i^{[k]}\|}{\|\Delta_i^{[k]}\|}$, $\|\Delta_i^{[k]}\| := (\sum_{i=0}^J \|\Delta_i^{[k]}\|^2)^{1/2}$, which satisfy $0 \le p_i^{[k]} \le 1$. Each block is then independently activated according to $\epsilon_i^{[k]} \sim B(p_i^{[k]})$, with higher probability thus assigned to blocks associated with larger proximal updates. Another advantage of this update rule is its flexibility: it does not require fixing in advance the number of blocks updated at each iteration, allowing this number to change from one iteration to the next.

#### B. Computational efficiency of MAGIC-FB

At first glance, the greedy Gauss-Southwell rule can seem too expensive to effectively reduce the convergence time, since it requires to compute $\Delta_i^{[k]}$ for all the indices i. We detail in this subsection how, in the context of linear inverse problems, the structure of the objective function allows for fast computation of these quantities. Let $w \in \mathbb{R}^n$ be an image expressed in the wavelet domain. We denote by $S_i$ the wavelet subspace associated with the block i of size $d_i$, and $\Pi_i : \mathbb{R}^n \to S_i$ the wavelet projection operator. For any image $x \in \mathbb{R}^n$, $w_i = \Pi_i x$. We have

$$ \begin{aligned} \nabla f(w) &= (AW^{\top})^{\top} (AW^{\top} w - y) \\ &= WA^{\top} A \left( \sum_{j=0}^{J} \Pi_j^{\top} w_j \right) - WA^{\top} y \\ &= \sum_{j=0}^{J} WA^{\top} A \Pi_j^{\top} w_j - WA^{\top} y. \end{aligned} $$

The partial gradient with respect to the block $0 \le i \le J$ can therefore be written as

$$ \nabla_{w_i} f(w) = \sum_{j=0}^{J} \Pi_i A^{\top} A \Pi_j^{\top} w_j - \Pi_i A^{\top} y. $$

The matrices $(\Pi_i A^{\top} A \Pi_j^{\top})_{0 \le i, j \le J}$ and the constant terms $(\Pi_i A^{\top} y)_{0 \le i \le J}$ can be pre-computed and stored.

When evaluating $\Delta_i^{[k]}$ at iteration k, all the terms $(\Pi_i A^{\top} A \Pi_j^{\top} w_j^{[k]})_{0 \le i, j \le J}$ can be stored. If block $i^{[k]}$ is chosen and updated, in order to update the value of the gradient for the next iteration, it is sufficient to recompute the matrix-vector products of the terms that depend on block $i^{[k]}$

$$ (\Pi_i A^{\top} A \Pi_{i[k]}^{\top} w_{i[k]}^{[k+1]})_{0 \le i \le J}, $$

as all the other partial products $(\Pi_i A^{\top} A \Pi_j^{\top} w_j^{[k]})_{0 \le i, j \le J}$ remain unchanged, since $w_j^{[k]} = w_j^{[k+1]}$ for $j \ne i^{[k]}$.

This allows the partial gradients $(\nabla_{w_i} f(w))_{0 \le i \le J}$ to be evaluated efficiently from the current block coefficients $(w_i)_{0 \le i \le J}$ without recomputing global matrix-vector products, thus making the block-selection step computationally inexpensive.

As a result, the proposed greedy block-selection strategy does not increase the complexity compared to standard block-coordinate proximal gradient methods. Instead, it reallocates computational effort by selecting the most relevant block before performing the costly matrix-vector products associated with the update, thereby improving overall efficiency.

### IV. EXPERIMENTS

In this section, we evaluate the performance of the proposed algorithm, which is assessed using a performance profile [19] shown on Figure 2, complemented by representative extreme cases for qualitative analysis on Figure 1. All experiments were conducted on one of the Blaise Pascal Center's machines [20], equipped with an NVIDIA RTX 3090 Ti GPU, an Intel Xeon Gold 6544Y CPU, and 125 GB of RAM. Code is available at https://github.com/EdgarDesainteMareville/MAGIC-FB

#### A. Experimental setup

a) Test dataset: Experiments are conducted on 1024 × 1024 cropped regions extracted from 100 images of the DIV2K dataset [21]. For each test instance, the blur and noise levels are randomly drawn to span a broad range of degradation conditions.

b) Degradation model: The forward operator A corresponds to a spatially invariant Gaussian blur with variance $\sigma_{blur}$ implemented as a separable convolution with circular boundary conditions. White Gaussian noise with variance $\sigma_{noise}$ is added to the observations. Blur levels are uniformly sampled with $\sigma_{blur} \in [1, 15]$ and noise levels are log-uniformly sampled with $\sigma_{noise} \in [10^{-3}, 10^{-1}]$.

c) Multiresolution representation and regularization: We use an orthonormal Daubechies-8 wavelet transform with $J = 5$ levels to decompose the images. We use the classical sparsity promoting $\ell_1$ norm of the wavelet coefficients as a regularization term, balanced by a regularization coefficient
$\lambda > 0$: for all $i \in \{1, \dots, J\}$, $g_i = \lambda \| \cdot \|_1$ and $g_0 \equiv 0$. This choice naturally satisfies (4). Note that the approximation coefficients $w_0$ are not penalized, as they are typically not sparse. The regularization parameter is shared across all scales and is selected by a grid search to maximize the peak signal-to-noise ratio (PSNR) with respect to the ground truth.

d) *Algorithms and parameters:* We compare the following optimization schemes:

FB: classical forward-backward splitting.

Stoc-BC-FB: uniform sampling stochastic multiresolution BC-FB method, where each block has a 50% chance of being updated at each iteration.

MLFB: Multilevel Forward-Backward – fixed hierarchical coarse-to-fine multiresolution BC-FB from [15].

MAGIC-FB: the proposed stochastic multiresolution adaptive proximal Gauss-Southwell BC-FB method.

All the methods are initialized with the observed image $w^{[0]} = Wy$ and use the same stepsize $\gamma = \frac{1.9}{\|A\|^2}$, where $\|A\|$ is estimated numerically. All methods are run for 200 iterations, which are sufficient for the methods to converge in all the considered settings.

## B. Numerical results

Figure 1 reports the convergence behavior of the different optimization schemes for two extreme degradation regimes, in terms of the objective function value as a function of wall-clock time. The time required to precompute the matrix products $(\Pi_i A^\top A \Pi_j)_{i,j}$ is excluded from the reported runtimes, as these quantities are shared by all methods. For the stochastic methods, results are reported as mean and standard deviation over 10 independent runs.

In the strong blur / low noise regime where $\sigma_{\text{blur}} = 7$ and $\sigma_{\text{noise}} = 0.01$, MLFB outperforms FB and Stoc-FB, highlighting the relevance of coarse-scale updates when the degradation is a low-pass filter, as the gain associated with the fitting term outweighs the decrease in the objective function. In this setting, MAGIC-FB closely matches the behavior of MLFB, while Stoc-FB is ineffective.

In the low blur / high noise regime, where $\sigma_{\text{blur}} = 1$ and $\sigma_{\text{noise}} = 0.1$, FB is the most effective method. While MLFB does not perform well, the stochastic methods Stoc-FB and MAGIC-FB achieve fast convergence, showing the importance of parallel updates and exploration allowed by randomness.

The proposed MAGIC-FB method shows to be able to adapt to all degradation regimes, in contrast to its competitors whose performance strongly vary with the scenario.

The center top panels of Figure 1 show the empirical block activation frequencies for the multiresolution methods (excluding FB) with respect to the iterations. For the stochastic update rules, activation frequencies are averaged across runs. These heatmaps provide insight into how block selection evolves under different degradation regimes.

In the strong blur / low noise regime, MAGIC-FB essentially activates coarse-scale blocks, evidencing that, at the beginning, there is little to be gained by updating fine-scale blocks. Then, MAGIC-FB gradually shifts towards the fine scale updates,

naturally reproducing the MLFB updating rule, and achieving comparable performance.

Conversely, in the low blur / high noise regime, block activations are more evenly distributed across scales, suggesting that restoring simultaneously blocks at all scales contributes to improving the objective function.

Intermediate reconstructions at a fixed time budget of $t = 1$ s, together with their corresponding PSNR values, are displayed to provide a qualitative comparison of reconstruction quality across methods.

The performance profile shown on Figure 2 provides an evaluation of the methods efficiency across the test set of 100 DIV2K images, illustrating the proportion of problems for which each method is within a certain factor of the best performance, measured in terms of loss function value obtained after 1 second of execution. The proposed MAGIC-FB strategy attains the lowest objective value within the prescribed time budget on over 65% of the test instances, while maintaining strong robustness across performance ratios compared to the other methods.

Overall, these results show that the proposed strategy dynamically reallocates computational effort across scales and recovers, in each regime, a block-selection pattern consistent with the structure of the underlying inverse problem.

## V. CONCLUSION

In this work, we introduced MAGIC-FB, an adaptive multiresolution block-coordinate Forward-Backward method for image restoration problems. By leveraging a wavelet-based block decomposition and an adaptive block-selection driven by update magnitudes via a stochastic proximal Gauss-Southwell rule, the proposed algorithm automatically balances updates across scales. This allows for an automatic adaptation of the algorithm to different degradation regimes, ranging from strongly blurred to noise-dominated images, without relying on a fixed schedule.

Experimental results on image deblurring problems highlight the limitations of fixed update rules in multiresolution strategies when blur or noise conditions vary, and demonstrate that the proposed stochastic approach matches the best methods in the two extreme regimes. In particular, the adaptive mechanism recovers coarse-to-fine behavior in strongly blurred settings, while progressively adding fine-scale updates as noise becomes dominant.

In future work we plan to investigate the theoretical properties of the proposed adaptive scheme, including convergence under weak convexity assumptions, as well as to explore its extension to Plug-and-Play priors to improve simultaneously restoration performance and execution time.

## ACKNOWLEDGMENT

This work was funded by the ANR-24-CE23-7039 MEPHISTO project, the Simone and Cino Del Duca foundation and the PEPR IA project. We thank the Blaise Pascal Center for its computational support, using the SIDUS [20] solution.
<figure><img src="image_2.png" /><figcaption>Fig. 1. Convergence behavior and block-selection patterns for the considered methods on an image deblurring problem. Left: $\sigma_{\text{blur}} = 7$ and $\sigma_{\text{noise}} = 0.01$. Right: $\sigma_{\text{blur}} = 1$ and $\sigma_{\text{noise}} = 0.1$. Top left and right panels: objective function value as a function of time. Centered top panels: block activation patterns across wavelet scales. Bottom panels: image reconstructions at $t = 1$ s for each method.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2. Performance profiles of the considered methods across the 100 test instances. The x-axis represents the factor $\beta$ by which a method's performance is compared to the best method for each problem, while the y-axis shows the proportion of problems for which the method is within this factor of the best performance (the closer to the upper left corner, the best).</figcaption></figure>

REFERENCES

[1] P. L. Combettes and V. R. Wajs, "Signal recovery by proximal forward-backward splitting," *Multiscale Model. and Simul.*, vol. 4, no. 4, pp. 1168-1200, 2005.

[2] P. L. Combettes and J.-C. Pesquet, "Proximal splitting methods in signal processing," in *Fixed-Point Algo. for Inverse Problems in Sci. and Eng.* Springer, 2011, pp. 185-212.

[3] A. Beck and M. Teboulle, "A fast iterative shrinkage-thresholding algorithm for linear inverse problems," *SIAM J. Imaging Sci.*, vol. 2, no. 1, pp. 183-202, 2009.

[4] S. V. Venkatakrishnan, C. A. Bouman, and B. Wohlberg, "Plug-and-play priors for model based reconstruction," in *2013 IEEE Global Conf. on Sig. and Inform. Process.*, 2013, pp. 945-948.

[5] S. Hurault, A. Leclaire, and N. Papadakis, "Gradient step denoiser for convergent plug-and-play," in *ICLR*, 2022.

[6] J.-C. Pesquet, A. Repetti, M. Terris, and Y. Wiaux, "Learning maximally monotone operators for image recovery," *SIAM J. Imaging Sci.*, vol. 14, no. 3, pp. 1206-1237, 2021.

[7] P. Tseng, "Convergence of a block coordinate descent method for nondifferentiable minimization," *J. Optim. Theory Appl.*, vol. 109, no. 3, pp. 475-494, 2001.

[8] Y. Nesterov, "Efficiency of coordinate descent methods on huge-scale optimization problems," *SIAM J. Opt.*, vol. 22, no. 2, pp. 341-362, 2012.

[9] S. Salzo and S. Villa, "Parallel random block-coordinate forward-backward algorithm: A unified convergence analysis," *Math. Program.*, vol. 193, no. 1, pp. 225-269, 2022.

[10] R. V. Southwell, "Relaxation methods in engineering science - a treatise on approximate computation," *The Mathematical Gazette*, vol. 25, no. 265, pp. 180-182, 1941.

[11] J. Nutini, M. Schmidt, I. Laradji, M. Friedlander, and H. Koepke, "Coordinate descent converges faster with the gauss-southwell rule than random selection," in *ICML*, PMLR, 2015.

[12] B. Pascal, N. Pustelnik, P. Abry, and J.-C. Pesquet, "Block-coordinate proximal algorithms for scale-free texture segmentation," in *Proc. Int. Conf. Acoust., Speech Signal Process.*, IEEE, 2018, pp. 1253-1257.

[13] E. Chouzenoux, J.-C. Pesquet, and A. Repetti, "A block coordinate variable metric forward-backward algorithm," *J. Global Optim.*, vol. 66, no. 3, pp. 457-485, 2016.

[14] G. Lauga, E. Riccietti, N. Pustelnik, and P. Gonçalves, "IML FISTA: A multilevel framework for inexact and inertial forward-backward. Application to image restoration," *SIAM J. Imaging Sci.*, vol. 17, no. 3, pp. 1347-1376, 2024.

[15] L. Briceño-Arias, P. Gonçalves, G. Lauga, N. Pustelnik, and E. Riccietti, "A flexible block-coordinate forward-backward algorithm for non-smooth and non-convex optimization," *arXiv preprint*, 2025.

[16] Z. Kadkhodaie, F. Guth, S. Mallat, and E. P. Simoncelli, "Learning multi-scale local conditional probability models of images," *ICLR*, 2023.

[17] N. Laurent, J. Tachella, E. Riccietti, and N. Pustelnik, "Multi-level plug-and-play image restoration," *IEEE Trans. Comput. Imaging*, 2025.

[18] L. Chaari, N. Pustelnik, C. Chaux, and J.-C. Pesquet, "Solving inverse problems with overcomplete transforms and convex optimization techniques," in *Wavelets XIII*, SPIE, vol. 7446, 2009, pp. 254-267.

[19] E. D. Dolan and J. J. Moré, "Benchmarking optimization software with performance profiles," *Math. Program.*, vol. 91, no. 2, pp. 201-213, Jan. 2002.

[20] E. Quemener and M. Corvellec, "Sidus—the solution for extreme deduplication of an operating system," *Linux J.*, vol. 2013, no. 235, 2013.

[21] E. Agustsson and R. Timofte, "NTIRE 2017 challenge on single image super-resolution: Dataset and study," in *CVPR Workshops*, 2017.