# A flexible block-coordinate forward-backward algorithm for non-smooth and non-convex optimization.

Luis Briceño-Arias, Paulo Gonçalves, Guillaume Lauga, Nelly Pustelnik, Elisa Riccietti

► To cite this version:

Luis Briceño-Arias, Paulo Gonçalves, Guillaume Lauga, Nelly Pustelnik, Elisa Riccietti. A flexible block-coordinate forward-backward algorithm for non-smooth and non-convex optimization. 2025. hal-05337116

HAL Id: hal-05337116

https://inria.hal.science/hal-05337116v1

Preprint submitted on 29 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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 Attribution 4.0 International License
A flexible block-coordinate forward-backward
algorithm for non-smooth and non-convex
optimization.

Luis Briceño-Arias<fnref n="2" />†, Paulo Gonçalves<fnref n="1" />†, Guillaume Lauga<fnref n="1" />*,
Nelly Pustelnik<fnref n="3" />†, Elisa Riccietti<fnref n="1" />†

Abstract

Block coordinate descent (BCD) methods are prevalent in large scale optimiza-
tion problems due to the low memory and computational costs per iteration, the
predisposition to parallelization, and the ability to exploit the structure of the
problem. The theoretical and practical performance of BCD relies heavily on the
rules defining the choice of the blocks to be updated at each iteration. We pro-
pose a new deterministic BCD framework that allows for very flexible updates,
while guaranteeing state-of-the-art convergence guarantees on non-smooth non-
convex optimization problems. While encompassing several update rules from
the literature, this framework allows for priority on updates of particular blocks
and correlations in the block selection between iterations, which is not permitted
under the classical convergent stochastic framework. This flexibility is lever-
aged in the context of multilevel optimization algorithms and, in particular, in
multilevel image restoration problems, where the efficiency of the approach is
illustrated.

<footnote n="1">¹*ENS de Lyon, CNRS, Université Claude Bernard Lyon 1, Inria, LIP, UMR 5668, 69342, Lyon cedex 07, France.</footnote>

<footnote n="2">²Departamento de Matemática, Universidad Técnica Federico Santa María, Santiago, Chile.</footnote>

<footnote n="3">³Laboratoire de Physique, ENSL, CNRS UMR 5672, F-69342, Lyon,
France.</footnote>

<footnote>*Corresponding author(s). E-mail(s): guillaume.lauga@ens-lyon.fr;
Contributing authors: luis.briceno@usm.cl; paulo.goncalves@inria.fr;
nelly.pustelnik@ens-lyon.fr; elisa.riccietti@ens-lyon.fr;</footnote>

<footnote>†These authors contributed equally to this work.</footnote>
# 1 Introduction

In this paper we introduce a new flexible block-coordinate algorithm to solve the separable-structured optimization problem

$$ \hat{\mathbf{x}} \in \underset{\mathbf{x}=(x_1, \dots, x_L) \in \mathcal{H}}{\operatorname{Argmin}} \Psi(\mathbf{x}) := f(\mathbf{x}) + \sum_{\ell=1}^{L} g_{\ell}(x_{\ell}), \quad (1) $$

where $\mathcal{H}$ is the direct sum of real separable, and finite dimensional Hilbert spaces $(\mathcal{H}_\ell)_{1 \le \ell \le L}$, $f: \mathcal{H} \to (-\infty, +\infty]$ is continuously differentiable, and, for every $\ell \in \{1, \dots, L\}$, $g_\ell: \mathcal{H}_\ell \to (-\infty, +\infty]$ is proper and lower semicontinuous. Without additional assumptions, the minimization problem can be non-smooth and non-convex. In image processing, for instance, $f$ usually encodes a fidelity with respect to some observation (e.g., a corrupted image) and functions $(g_\ell)_{1 \le \ell \le L}$ encode some prior knowledge about components of the parameters to estimate (e.g., the regularity of an image). In large-scale optimization, when the dimension of $\mathcal{H}$ is high, block-coordinate (BC) methods are widely used for their low per-iteration computational and memory costs, and their ability to exploit problem separability. The main application of BC approaches is for parameter estimation (e.g., standard machine learning [1–9] but also deep learning [10–12]). The theoretical convergence of block-coordinate methods depends on the choice of (i) the update rule of blocks $x_1, \dots, x_L$ and (ii) the function to be minimized for each block. The update rule can be of two types: stochastic or deterministic. The function to be minimized for each block can either be $\Psi$ or a linearized version of $\Psi$. The first case is referred as Block-Coordinate Descent (see, e.g., [13]), which includes the Gauss-Seidel approach (see [14] and references therein), while the second case is Block-Coordinate Proximal Gradient (BC-PG), whose particular case when $g_\ell \equiv 0$ reduces to Block-Coordinate Gradient Descent.

In [14], the linear convergence of the Block-Coordinate Descent method is proved under the strong convexity of $\Psi$ w.r.t. each block. However, this approach has two drawbacks: for each block the minimization procedure can be as complex as solving (1) and without strong convexity assumption there is no theoretical guarantees of its convergence [15]. For these reasons, we focus in this paper on the BC-PG, which offers efficient update of each block due to linearization and convergence guarantees in a larger setting.

## State-of-the-art on BC-PG.

There exists a wide literature using stochastic activation of blocks in the convex and non-convex setting, see, e.g., [1, 2, 9, 16–26]. Most of the literature on stochastic approaches studies convergence or rate of convergence of the objective function values in expectation [2, 19–21]. In [17], the almost sure convergence of the random iterates
is proved by using the concept of stochastic Quasi-Féjer sequence, first introduced in [27]. This framework is powerful and can be applied to many types of block-coordinate algorithms (see for instance primal-dual ones in [16, 28, 29]). However, it is thus far only applicable if $\Psi$ is convex.

The main limitation of stochastic approaches is the random block selection, which does not offer the flexibility to prescribe the order of block activation, even though it allows for parallel activations. Moreover, a comparison between random and deterministic block activation strategies in the strongly convex setting when $g_\ell \equiv 0$ is presented in [8], where the deterministic Gauss-Southwell rule, i.e., the strategy that updates at each iteration the block with the largest partial gradient, is shown to be preferable in practice.

Deterministic block-activation algorithms appear in a wide literature that covers both convex and non convex frameworks, with convergence supported by theoretical arguments, different from those used in the stochastic setting. Two main types of deterministic rules exist: cyclic and essentially cyclic. Cyclic refers to the sequential update of one block after the other until every one of them has been updated once. These updates encompass alternated linearized optimization techniques [30]. This approach does not permit to update the same component twice in the same cycle. Essentially cyclic rules on the other hand, also called K-cyclic [31] for $K \ge 1$, update every block at least once in any K consecutive iterations. This rule allows for the activation of a single block multiple times in a same cycle of K iterations.

Under convexity assumptions, some convergence results are available for cyclic BCD in [31–33]. These results either cover a different framework than ours (i.e., (1)) or a more restricted one. More precisely, when $g_\ell$ is the indicator function of an interval and $f$ is coordinatewise strongly convex, the linear convergence of the BC-PG with cyclic updates is proved in [32]. In [31], the convergence of a convex feasibility problem with several updating rules involving cyclic and essentially cyclic approaches is obtained. In [34], linear convergence rates of the functional values and iterates are derived under strong convexity assumptions and when the separability is in the smooth component. A similar strategy is used in [33] for minimizing the average of a finite number of strongly convex functions.

BC-PG with Gauss-Southwell rule also belongs to the class of deterministic block-activation algorithms but to the best of our knowledge the results are based on convergence of the values of the function and not of the iterates [8, 35].

Guarantees of convergence of cyclic/essentially cyclic BCD have been investigated in the nonconvex setting for instance in [6, 30, 36] under Kurdyka-Łojasiewicz (KŁ) or Łojasiewicz properties/inequalities [30, 37–39]. The authors proved the convergence to a minimizer or a critical point of $\Psi$, and derived the rate of convergence of the sequence of iterates under specific assumptions on the desingularizing function. Among the existing literature, the work most closely related to ours is [36], where the authors examine convergence analysis in the essentially cyclic case. However, their study focuses on updating a single block at each iteration, which does not encompass the parallel essentially cyclic framework that allows for the simultaneous updating of multiple blocks.
### BC-PG approaches for image reconstruction.

If BC-PG approaches are widely used in hyperparameter estimation (e.g., machine learning [1–9] and deep learning [10–12]), its application is less straightforward in image reconstruction. Indeed, the image reconstruction is modeled by a minimization problem of the form (1), whose global structure is not handled well by BC-PG methods that operate through local (patch-wise or pixel-wise) updates [2]. BC-PG is more effective for images dominated by local information, such as astronomical images, which are composed mostly of a black background with sparse point sources relative to the image size [40]. For more general images, the structure of the optimization problem can be exploited to design BC-PG schemes. For instance in [41] the authors divide the image in lattices in order to benefit from a reduced size allowing for faster convergence. These lattices are inherited from the structure of the total variation regularization, which penalizes the difference of the value of a given pixel with its neighbor on the right, and its neighbor below. Therefore, grouping pixels by selecting one every other two rows and every other two columns define four independent groups of pixels that cover the entirety of the image. In [36], blocks are defined by combining wavelet coefficients locally in an overcomplete dictionary.

In image processing, block strategies have mainly been implemented to manage large volumes of data. However, they do not necessarily accelerate the solution process for a given data volume. For example, the forward-backward algorithm may remain faster than the BC-PG algorithm if all the data is accessible (see, e.g., [42]).

### Multilevel algorithms for image reconstruction.

There exists another class of algorithms that exploit the structure of the minimization problem and have demonstrated their effectiveness on imaging problems: multilevel algorithms [43–48]. Such procedures can speed-up algorithms in the convex framework [49–51] but also converge towards a better solution in the non-convex setting [52, 53]. A multilevel algorithm tackles high dimensional optimization problems by defining a hierarchy of smaller dimensional approximations of the original problem and by alternating optimization steps on this hierarchy. This approach, similarly to BC-PG, relies on a problem decomposition to perform iterations at lower complexity.

As compared to BC-PG, multilevel methods are much faster for imaging problems [50]. However, the notion of coherence between approximations is at the core of the multilevel strategy and is not used for BC-PG algorithms [50]. Finally, the convergence for multilevel methods is limited to the convex possibly non-smooth setting or the non-convex smooth setting, while the convergence of BC-PG methods is established both in the convex and the non-convex setting, without smoothness assumptions.

### Contribution: Flexible block-coordinate proximal gradient approach.

We propose a new BC-PG algorithm that relies on a parallel and essentially cyclic rule, which enables us to fully exploit the structure of the optimization problem, whenever it exists. We will refer to it as FLEX – BC – PG, which stands for *Flexible BC-PG*.

FLEX – BC – PG can trigger parallel activation of the blocks. Such parallel updates allow for instance to update blocks that the structure of the problem groups together in their contribution to the objective function (like with Total Variation in [41]) or share
a position in a hierarchy (such as the detail coefficients of a wavelet transform [54]). For instance, with respect to the existing cyclic update framework, FLEX – BC – PG, allows for varying the size of the blocks along the iterations, which is of interest since using larger block size is beneficial to speed up the optimization [7], but may be more costly. Alternating between the two thus allows to benefit from the advantages of both.

Our first contribution is to prove that the convergence of FLEX – BC – PG to a critical point of the optimization problem is guaranteed. To the best of our knowledge, no BC-PG algorithm was shown to converge to a critical point of a non-smooth and non-convex optimization problem while allowing parallel updates in a non-stochastic setting. Our FLEX – BC – PG algorithm is thus the first that has this property, while having state-of-the-art convergence guarantees (i.e., decrease of objective function values and convergence to a critical point in a non-convex setting) We display in Figure 1 some of the update schemes covered by our framework.

Our second contribution is to study the connection between FLEX – BC – PG and multilevel algorithms. We show that, for a specific image restoration problem that has a hierarchical structure (cf. (30)), the proposed update can encompass a hierarchical choice of the blocks, mimicking multilevel iterations. In particular, we prove that the first-order coherence, usually imposed in multilevel methods, is actually essential for establishing the equivalence of the two methods. As a result, just like multilevel methods, the FLEX – BC – PG algorithm we propose is faster than the classical forward-backward (FB) in practice. Moreover, thanks to this connection, the theory presented in this article can be seen as an extension of the convergence theory of multilevel FB [50] to a non-convex setting.

## Outline.

This article is organized as follows. In Section 2, we present the state-of-the art BCD strategies to tackle Problem (1) and the proposed FLEX – BC – PG algorithm. In Section 3, we prove the convergence of FLEX – BC – PG in a non-smooth and non-convex setting, where we assume that $\Psi$ satisfies the Kurdyka-Łojasiewicz inequality. In Section 4 we show, for a particular instance of problem (1), that a multilevel forward-backward algorithm can be seen as an instance of the FLEX – BC – PG algorithm with a hierarchical update rule. Finally in Section 5, we compare several instances of our FLEX – BC – PG algorithm to BC-PG strategies encountered in the literature.

## Notations.

We introduce the notations that will be used in the following. We use bold letters to consider full vectors of variables $(\mathbf{x})$ and plain letters indexed by $\ell$ to denote blocks $(x_\ell)$. We denote by $\|\cdot\|$ the Euclidean norm on $\mathcal{H}$ and by $\|\cdot\|$ the Euclidean norm on the $L$ spaces $(\mathcal{H}_\ell)_{1\le\ell\le L}$. Similarly, the scalar product on $\mathcal{H}$ will be denoted by $\langle\langle \cdot, \cdot \rangle\rangle$ and the scalar product on $\mathcal{H}_\ell$ by $\langle \cdot, \cdot \rangle$; the potential ambiguity between two spaces is cleared up as the variables on which the scalar product is applied will be indexed by $\ell$. Note that for all $\mathbf{x}$ and $\mathbf{y}$ in $\mathcal{H}$, $\langle\langle \mathbf{x}, \mathbf{y} \rangle\rangle = \sum_{\ell=1}^L \langle x_\ell, y_\ell \rangle$. For a continuously differentiable function $f$, $\nabla_\ell f(\mathbf{x})$ is the gradient of $f$ with respect to the variables in the $\ell$-th block at some $\mathbf{x} \in \mathcal{H}$. Recall that, for a function $g$ that is proper and lower semicontinuous with $\inf_{\mathcal{H}} g > -\infty$, the proximity operator of $g$ of parameter $\tau > 0$ is
<figure><img src="image_2.png" /><figcaption>Fig. 1: Examples of update rules that are covered by our proposed algorithm FLEX - BC - PG for a problem with 4 blocks. The top row displays existing rules that are covered by our framework, the bottom displays new rules that are now covered by our framework. The variable of the ℓ-th block at iteration n is indexed as $x_{\ell}^{n}$. At the top of each scheme, we display the iteration number n, and at the bottom the cycle number k. At each iteration we highlight in red the blocks that are updated. Thus, each column depicts the activation of the blocks at a given iteration. All these rules share a common feature, necessary to make the corresponding BCD method convergent: each block must have been updated at least once during a cycle. Here, each cycle contains 4 iterations, but it is not necessary for the number of iterations per cycle to be equal to the number of blocks.</figcaption></figure>

defined as

$$ (\forall \mathbf{x} \in \mathcal{H}), \quad \text{prox}_{\tau g}(\mathbf{x}) = \underset{\mathbf{u} \in \mathcal{H}}{\operatorname{Argmin}} \frac{1}{2\tau} \|\mathbf{u} - \mathbf{x}\|^2 + g(\mathbf{u}). \quad (2) $$

where $\text{Argmin}_{u \in \mathcal{H}} h(\mathbf{u})$ denotes the set of minimizers of some function $h: \mathcal{H} \to (-\infty, +\infty]$. Given $(\mathbf{x}, \mathbf{u}) \in \mathcal{H} \times \mathcal{H}$, we have [30]

$$ \mathbf{u} \in \text{prox}_g(\mathbf{x}) \implies \mathbf{x} - \mathbf{u} \in \partial g(\mathbf{u}). \quad (3) $$

We will say that a function g is proximable when $\text{prox}_{\tau g}$ is known under closed form.
## 2 Block-coordinate forward-backward algorithm

The idea of splitting an optimization problem into smaller tasks is ubiquitous in practice and has sparked in the last twenty years a lot of research to better understand its potential from a theoretical perspective. The following paragraphs describe the bulk of these studies in the context of the BC-PG algorithm, where block updates are done using proximal-gradient descent to handle the non-smoothness of $g := \sum_{\ell=1}^{L} g_{\ell}$. A complete overview of the update methods for BC-PG algorithm may be found in [2, 7].

### 2.1 State-of-the-art strategies

The most standard formulation of block-coordinate forward-backward algorithm is the following. We index the sequence of iterates by a superscript denoting the iteration number and a subscript denoting the block of variables. Thus, $x_{\ell}^n$ denotes the $\ell$-th block at the $n$-th iteration. For convenience, we write $\mathbf{x}^n$ to denote the full variable at iteration $n$, so that $\mathbf{x}^n = (x_1^n, \dots, x_L^n)$. We denote $(\varepsilon^n)_{n \in \mathbb{N}} = (\varepsilon_1^n, \dots, \varepsilon_L^n)_{n \in \mathbb{N}}$ a sequence of variables with value in $\{0, 1\}^L$ and the step-sizes $(\tau_{\ell}^n)_{1 \le \ell \le L} \in \mathbb{R}_{++}^L$, for all $n \in \mathbb{N}$. The algorithm is initialized with $\mathbf{x}^0 = (x_1^0, \dots, x_L^0) \in \text{dom } g$ and reads

$$ \text{for } n = 0, 1, \dots \\ \quad \left[ \begin{array}{l} \text{for } \ell = 1, \dots, L \\ \quad \left[ x_{\ell}^{n+1} \in x_{\ell}^n + \varepsilon_{\ell}^n \left( \text{prox}_{\tau_{\ell}^n g_{\ell}} (x_{\ell}^n - \tau_{\ell}^n \nabla_{\ell} f(\mathbf{x}^n)) - x_{\ell}^n \right). \end{array} \right. \tag{4} $$

The settings investigated in the literature are the following ones:

1. Stochastic setting: $(\varepsilon_1^n, \dots, \varepsilon_L^n) \in \{0, 1\}^L$ are chosen randomly, enabling random parallel updates, for all $n \in \mathbb{N}$.

2. Essentially cyclic setting: given $K > 0$,

$$ (\forall j \in \mathbb{N}), \quad \bigcup_{n=j}^{j+K-1} \{\ell \mid \varepsilon_{\ell}^j = 1\} = \{1, \dots, L\}. $$

There have been numerous works to study the algorithm in (4) in the stochastic setting [2, 9, 16, 17, 19–25], in the essentially cyclic one [1, 2, 6, 8, 30, 36, 55–57] and finally using greedy rules that are efficient in practice but lack the guarantees of the previous two. This list of references is not exhaustive, but it is representative of the proof techniques used to study the convergence of algorithm (4).

### 2.2 Proposed setting

As highlighted in the introduction, despite the huge literature on block coordinate methods, no algorithm allowing a deterministic control of parallel updates with convergence guarantees has yet to be proposed.

We are also motivated to draw connections between multilevel algorithms and BCD algorithms, and the method we propose can benefit in the relevant setting (see Section 5) from the efficiency of multilevel algorithms and the convergence guarantees
of the BCD formalism. Specifically, we design a convergent block-coordinate descent algorithm for non-smooth and non-convex optimization where the updates are potentially parallel, essentially cyclic, and may be randomly reshuffled at each cycle. Our framework encompasses, for instance, BCD algorithms where the size of the blocks may vary from one iteration to the other [7] or hierarchical BCD, whose update rule mimics the behavior of multilevel algorithms (i.e., some blocks, deemed more important, are updated more often than others, see Figure 1 or Section 4). This requires to allow correlated block selection, which is not permitted under the classical convergent stochastic framework [17]. Our major contribution is therefore the convergence of the method both in function values and with respect to the set of critical points.

### FLEX-BC-PG iterations

Our update rule is parallel and essentially cyclic by setting a priori the sequence $\epsilon^n = (\epsilon_1^n, \dots, \epsilon_L^n) \in \{0, 1\}^L$ for all $n \in \mathbb{N}$. To simplify the following analysis, we rewrite algorithm (4) to explicitly incorporate the cycles. Let $K \in \mathbb{N} \setminus \{0\}$ be the number of iterations to complete one cycle. Let $(\tau_\ell^n)_{1 \le \ell \le L} \in \mathbb{R}_{++}^L$ for all $n \in \mathbb{N}$, and $\bar{\mathbf{x}}^0 = \mathbf{x}^0 = (x_1^0, \dots, x_L^0) \in \text{dom } g$. Set $k=0$. The iterations read

$$ \text{for } n = 0, 1, \dots \\ \quad \text{for } \ell = 1, \dots, L \\ \quad \quad \left[ \begin{array}{l} x_{\ell}^{n+1} = x_{\ell}^{n} + \epsilon_{\ell}^{n} \left( \text{prox}_{\tau_{\ell}^{n} g_{\ell}} (x_{\ell}^{n} - \tau_{\ell}^{n} \nabla_{\ell} f(\mathbf{x}^{n})) - x_{\ell}^{n} \right). \\ \text{if } n+1 \equiv 0 \ [K] \\ \quad k = k+1 \\ \quad \bar{\mathbf{x}}^{k} = \mathbf{x}^{n+1} \end{array} \right. \tag{5} $$

The convergence analysis relies on a cyclic rule for the updates, and we assume that each cycle consists of at most $K$ iterations. We will denote $\bar{\mathbf{x}}^k$ the iterates obtained after $k$ cycles and thus $k \cdot K$ iterations to accentuate the difference with the iterates $\mathbf{x}^n$. We will prove the convergence of $(\bar{\mathbf{x}}^k)_{k \in \mathbb{N}}$ to a critical point, which will give us, as a byproduct, the convergence of $\mathbf{x}^n$ to the same critical point. Here, $\bar{\mathbf{x}}^k = \mathbf{x}^n$ when $n = k \cdot K$. These notations are illustrated in Figure 1.

## 3 Convergence of the proposed BC-PG algorithm

We now analyse the convergence of the proposed scheme using the Kurdyka-Łojasiewicz property of $\Psi$, which allows us to provide both convergence of the objective function values and of the iterates to a critical point of $\Psi$. For completeness of the argument, we prove the convergence of a stochastic version of our FLEX-BC-PG algorithm in Appendix B.

### 3.1 Preliminaries in non-convex optimization

In the non-convex setting, we will need an appropriate notion of subgradient. For the rest of the paper, we assume that $\mathcal{H} := \mathbb{R}^N$.
**Definition 1.** *Subdifferential [58]*. Let $\Psi: \mathcal{H} \to \mathbb{R}$, and let $\mathbf{x} \in \mathcal{H}$. The Fréchet subdifferential of $\Psi$ at $\mathbf{x}$ is denoted by $\hat{\partial}\Psi(\mathbf{x})$ and is given by the set

$$ \hat{\partial}\Psi(\mathbf{x}) = \left\{ \hat{\mathbf{s}} \in \mathcal{H} \mid \liminf_{\mathbf{y} \to \mathbf{x}, \mathbf{y} \neq \mathbf{x}} \frac{1}{\|\mathbf{x} - \mathbf{y}\|} (\Psi(\mathbf{y}) - \Psi(\mathbf{x}) - \langle\langle \mathbf{y} - \mathbf{x}, \hat{\mathbf{s}} \rangle\rangle) \geq 0 \right\}. $$

If $\mathbf{x} \notin \text{dom } \Psi$, then $\hat{\partial}\Psi(\mathbf{x}) = \emptyset$. The limiting subdifferential of $\Psi$ at $\mathbf{x}$ is denoted by $\partial\Psi(\mathbf{x})$ and is given by

$$ \partial\Psi(\mathbf{x}) = \{\mathbf{s} \in \mathcal{H} \mid \exists (\mathbf{x}^k, \hat{\mathbf{s}}^k) \xrightarrow{k \to \infty} (\mathbf{x}, \mathbf{s}) \\ \text{such that } \Psi(\mathbf{x}^k) \xrightarrow{k \to \infty} \Psi(\mathbf{x}) \text{ and } (\forall k \in \mathbb{N}) \hat{\mathbf{s}}^k \in \hat{\partial}\Psi(\mathbf{x}^k)\}. $$

Recall that if $\Psi$ is convex, its subdifferential is given, for all $\mathbf{x} \in \mathcal{H}$, by

$$ \partial\Psi(\mathbf{x}) = \{\mathbf{s} \in \mathcal{H} \mid \Psi(\mathbf{s}) + \langle\langle \mathbf{s}, \mathbf{y} - \mathbf{x} \rangle\rangle \leq \Psi(\mathbf{y}), \forall \mathbf{y} \in \mathcal{H}\}. $$

Both $\hat{\partial}\Psi(\mathbf{x})$ and $\partial\Psi(\mathbf{x})$ are closed sets [58, Theorem 8.6].

The limiting subdifferential benefits from the following separability property:

**Proposition 3.1.** *Subdifferentiability property [58]*. Let $\Psi$ be defined as in problem (1). Then, for all $\mathbf{x} = (x_1, \dots, x_L) \in \mathcal{H}_1 \times \dots \times \mathcal{H}_L$, we have

$$ \partial\Psi(\mathbf{x}) = (\nabla_1 f(\mathbf{x}) + \partial g_1(x_1)) \times \dots \times (\nabla_L f(\mathbf{x}) + \partial g_L(x_L)). $$

**The Kurdyka-Łojasiewicz (KŁ) property.**

A specific class of concave and continuous functions, called desingularizing functions, are of particular interest in the KL framework to handle non-convexity.

**Definition 2.** Concave and continuous functions [30]. Let $\eta \in (0, +\infty]$. We denote by $\Phi_\eta$ the class of all concave and continuous functions $\varphi: [0, \eta) \to \mathbb{R}_+$ that satisfy the following conditions:

1. $\varphi(0) = 0$,

2. $\varphi$ is $C^1$ on $(0, \eta)$ and continuous at $0$,

3. for all $s \in (0, \eta)$, $\varphi'(s) > 0$.

Now, we can introduce the definition of a KŁ function.

**Definition 3.** Kurdyka-Łojasiewicz (KŁ) property [30]. Let $\Psi: \mathcal{H} \to (-\infty, +\infty]$ be proper and lower semicontinuous.

1. The function $\Psi$ is said to have the KL property at $\tilde{\mathbf{u}} \in \text{dom } \partial\Psi$ if there exist $\eta \in (0, +\infty]$, a neighborhood $U$ of $\tilde{\mathbf{u}}$ and a function $\varphi \in \Phi_\eta$ such that for all

$$ \mathbf{u} \in U \cap [\Psi(\tilde{\mathbf{u}}) < \Psi(\mathbf{u}) < \Psi(\tilde{\mathbf{u}}) + \eta], $$

the following inequality holds

$$ \varphi'(\Psi(\mathbf{u}) - \Psi(\tilde{\mathbf{u}})) \text{dist}(0, \partial\Psi(\mathbf{u})) \geq 1. $$
Recall that $\mathrm{dist}(\mathbf{x}, \partial\psi(\mathbf{u})) = \inf_{\mathbf{s}\in\partial\psi(\mathbf{u})} \|\mathbf{s} - \mathbf{x}\|.$

2. If $\Psi$ satisfies the KL property at each point of dom $\partial\Psi$, then $\Psi$ is called a KL function.

The following lemma presents the KL property in a practical form by unifying the notion of neighborhood across its level curves.

**Lemma 3.2. Uniformized KL property [30].** Let $\Omega$ be a compact subset of $\mathcal{H}$. Let $\Psi: \mathcal{H} \to (-\infty, +\infty]$ be a proper and lower semicontinuous function, constant on $\Omega$ and satisfying the KL inequality on $\Omega$. Then there exists $\epsilon > 0, \eta > 0$, and $\varphi \in \Phi_\eta$ such that for all $\tilde{u} \in \Omega$ and all $u$ satisfying

$$ \{\mathbf{u} \in \mathcal{H} : \mathrm{dist}(u, \Omega) < \epsilon\} \text{ and } [\Psi(\tilde{u}) < \Psi(\mathbf{u}) < \Psi(\tilde{u}) + \eta] $$

one has

$$ \varphi'(\Psi(\mathbf{u}) - \Psi(\tilde{u}))\mathrm{dist}(0, \partial\Psi(\mathbf{u})) \geq 1. \quad (6) $$

**Remark 3.3.** The KL property is satisfied by numerous classes of functions, and notably by those considered in typical optimization settings. See [39] for an overview on this property.

## 3.2 Assumptions on the functions.

The convergence of FLEX - BC - PG relies on several classical assumptions that we present in the following.

### Assumption 1.

A1 $\Psi := f + g$, where $g := \sum_{\ell=1}^L g_\ell$, is coercive, and bounded below. For all $\ell \in \{1, \dots, L\}$, $g_\ell$ and $f$ are bounded from below.

A2 $\Psi$ satisfy the KL property (Definition 3 and Lemma 3.2).

### Assumption 2.

A3 For all $\ell \in \{1, \dots, L\}$, $g_\ell$ is a lower semicontinuous, proper function.

A4 $f$ is continuously differentiable and there exist constants $(\beta_{\ell,j})_{1\le\ell,j\le L}$ in $\mathbb{R}_{++}$ such that

$$ (\forall \ell, j \in \{1, \dots, L\})(\forall \mathbf{x} \in \mathcal{H})(\forall v_j \in \mathcal{H}_j) \quad \|\nabla_\ell f(\mathbf{x} + (0, \dots, 0, v_j, 0, \dots, 0)) - \nabla_\ell f(\mathbf{x})\| \le \beta_{\ell,j} \|v_j\| $$

Assumption A1 is sufficient to assert that the sequences generated by our algorithm are bounded [38].

Assumption A4 states that every partial gradient with respect to the block is Lipschitz continuous with respect to all the blocks, which is a quite stronger assumption than being Lipschitz continuous with respect only to its block. From this assumption we can derive multiple block smoothness, a common assumption in the BCD literature (e.g., [9, Assumption S1-S2-S3]). Despite this, Assumption A4 is fairly easy to verify in practice, since it is implied by the Lipschitz continuity of $\nabla f$ with constant $\beta_f$. Indeed, we can take $\beta_{\ell,j} = \beta_f$ for all $\ell, j$. Conversely, Assumption A4 implies that
$\nabla f$ is at most $\beta_f := \sqrt{\sum_{\ell,j=1}^L \beta_{\ell,j}^2}$-Lipschitz continuous. This is a consequence of the following proposition.

**Proposition 3.4.** *Multiple block smoothness. Suppose that Assumption 2 holds. For all $\boldsymbol{\varepsilon} = (\varepsilon_\ell)_{1\le\ell\le L} \in \{0,1\}^L$, there exists $\beta = \sqrt{\sum_{\ell,j=1}^L \varepsilon_j \beta_{\ell,j}^2} > 0$ such that for all $\mathbf{v} = (v_\ell)_{1\le\ell\le L} \in (\mathcal{H}_1, \dots, \mathcal{H}_L)$ we have*

$$ (\forall \mathbf{x} \in \mathcal{H}) \quad \|\nabla f(\mathbf{x} + \boldsymbol{\varepsilon} \odot \mathbf{v}) - \nabla f(\mathbf{x})\| \le \beta \|\boldsymbol{\varepsilon} \odot \mathbf{v}\|, $$

where $\odot$ is the element-wise multiplication $\boldsymbol{\varepsilon} \odot \mathbf{v} = (\varepsilon_1 v_1, \varepsilon_2 v_2, \dots, \varepsilon_L v_L)$.

*Proof.* Let $\mathbf{x}, \mathbf{v} = (v_\ell)_{\ell=1}^L \in \mathcal{H} = \bigoplus_{\ell=1}^L \mathcal{H}_\ell$. Note that

$$ \|\nabla f(\mathbf{x} + \mathbf{v}) - \nabla f(\mathbf{x})\|^2 = \sum_{\ell=1}^{L} \|\nabla_{\ell} f(\mathbf{x} + \mathbf{v}) - \nabla_{\ell} f(\mathbf{x})\|^2. $$

Now define, $\mathbf{v}_0 = \mathbf{0} \in \mathcal{H}$ and $\mathbf{v}_j = (v_1, \dots, v_j, 0, \dots, 0) \in \mathcal{H}$ for every $j \in \{1, \dots, L\}$. Note that $\mathbf{v}_L = \mathbf{v}$ and that

$$ (\forall j \in \{1, \dots, L\}) \quad \mathbf{v}_j - \mathbf{v}_{j-1} = (0, \dots, v_j, \dots, 0). $$

For every $\ell \in \{1, \dots, L\}$, the following equality holds

$$ \|\nabla_{\ell} f(\mathbf{x} + \boldsymbol{\varepsilon} \odot \mathbf{v}) - \nabla_{\ell} f(\mathbf{x})\| = \left\| \sum_{j=1}^{L} (\nabla_{\ell} f(\mathbf{x} + \boldsymbol{\varepsilon} \odot \mathbf{v}_j) - \nabla_{\ell} f(\mathbf{x} + \boldsymbol{\varepsilon} \odot \mathbf{v}_{j-1})) \right\| $$

Then, the triangular inequality and A4 imply

$$ \begin{aligned} \|\nabla_{\ell} f(\mathbf{x} + \boldsymbol{\varepsilon} \odot \mathbf{v}) - \nabla_{\ell} f(\mathbf{x})\| &\le \sum_{j=1}^{L} \|\nabla_{\ell} f(\mathbf{x} + \boldsymbol{\varepsilon} \odot \mathbf{v}_j) - \nabla_{\ell} f(\mathbf{x} + \boldsymbol{\varepsilon} \odot \mathbf{v}_{j-1})\| \\ &\le \sum_{j=1}^{L} \beta_{\ell,j} \|\varepsilon_j v_j\| \end{aligned} \quad (7) $$

and therefore, from Cauchy-Schwarz in $\mathbb{R}^L$,

$$ \|\nabla f(\mathbf{x} + \boldsymbol{\varepsilon} \odot \mathbf{v}) - \nabla f(\mathbf{x})\|^2 \le \sum_{\ell=1}^{L} \left( \sum_{j=1}^{L} \beta_{\ell,j} \|\varepsilon_j v_j\| \right)^2 \le \sum_{\ell,j=1}^{L} \beta_{\ell,j}^2 \|\boldsymbol{\varepsilon} \odot \mathbf{v}\|^2, \quad (8) $$

deducing that $\beta = \sqrt{\sum_{\ell,j=1}^L \varepsilon_j \beta_{\ell,j}^2}$ is a Lipschitz constant of $\nabla f$ with respect to the blocks selected by $\boldsymbol{\varepsilon}$. $\square$
### 3.3 Assumptions on the update rules.

We consider an essentially cyclic update scheme for the blocks in which parallel updates of different blocks may be used, paired with a potential shuffle of the updates order, as specified in the following assumption.

**Assumption 3.** A5 Let $\mathbb{I}^n$ be the set of the blocks updated at iteration $n$, i.e.,

$$ \mathbb{I}^n = \{\ell \mid \varepsilon_{\ell}^n = 1\} \subseteq \{1, \dots, L\}. \qquad (9) $$

There exists $K \in \mathbb{N} \setminus \{0\}$ such that

$$ (\forall j \in \mathbb{N}), \quad \bigcup_{n=j}^{j+K-1} \mathbb{I}^n = \{1, \dots, L\}. $$

**Remark 3.5.** Assumption A5 does not impose any constraint on the order of the updates of the blocks. For instance, it allows for sequential update of the blocks if $K = L$ and it reduces to the classical forward-backward update for $K = 1$. As an example of the flexibility of our framework, one can shuffle the order of the updates inside every cycle without breaking convergence guarantees. Hence, a random shuffle such as in [6] is compatible with Assumption A5.

Before setting the main result, we need one more technical lemma.

**Lemma 3.6.** Let $\{\bar{\mathbf{x}}^k\}_{k \in \mathbb{N}}$ be a sequence generated by Algorithm FLEX - BC - PG. Then,

$$ \|\bar{\mathbf{x}}^{k+1} - \bar{\mathbf{x}}^k\| \le \left( \sum_{n=k \cdot K}^{(k+1) \cdot K-1} \sum_{\ell \in \mathbb{I}^n} \|x_{\ell}^{n+1} - x_{\ell}^n\| \right), $$

where $\mathbb{I}^n$ is defined in (9).

*Proof.* Recall that

$$ \|\bar{\mathbf{x}}^{k+1} - \bar{\mathbf{x}}^k\| = \sqrt{\sum_{\ell=1}^{L} \|\bar{\mathbf{x}}_{\ell}^{k+1} - \bar{\mathbf{x}}_{\ell}^k\|^2} \le \sum_{\ell=1}^{L} \sqrt{\|\bar{\mathbf{x}}_{\ell}^{k+1} - \bar{\mathbf{x}}_{\ell}^k\|^2} = \sum_{\ell=1}^{L} \|\bar{\mathbf{x}}_{\ell}^{k+1} - \bar{\mathbf{x}}_{\ell}^k\|. $$

Now for all $\ell \in \{1, \dots, L\}$, the triangular inequality yields

$$ \|\bar{\mathbf{x}}_{\ell}^{k+1} - \bar{\mathbf{x}}_{\ell}^{k}\| \le \sum_{n=k \cdot K}^{(k+1) \cdot K-1} \|x_{\ell}^{n+1} - x_{\ell}^{n}\|. $$

Thus, summing up for all $\ell$, we obtain

$$ \|\bar{\mathbf{x}}^{k+1} - \bar{\mathbf{x}}^k\| \le \left( \sum_{n=k \cdot K}^{(k+1) \cdot K-1} \sum_{\ell=1}^{L} \|x_{\ell}^{n+1} - x_{\ell}^n\| \right) = \left( \sum_{n=k \cdot K}^{(k+1) \cdot K-1} \sum_{\ell \in \mathbb{I}^n} \|x_{\ell}^{n+1} - x_{\ell}^n\| \right) $$
where the last equality follows from $x_{\ell}^{n+1} = x_{\ell}^{n}$, for all $\ell \notin \mathbb{I}^{n}$ and $n \in \mathbb{N}$. $\square$

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

## 3.4 Main result: convergence to critical points of $\Psi$

The convergence to critical points of the sequence ${\bar{\mathbf{x}}^k}_{k \in \mathbb{N}}$ is a consequence of the sufficient decrease property of our algorithm, paired with the existence at each point of a subgradient bounded by the norm of the difference between iterates.

**Proposition 3.7.** Suppose that Assumptions 1, 2, and 3 hold. Let ${\mathbf{x}^n}_{n \in \mathbb{N}}$, and ${\bar{\mathbf{x}}^k}_{k \in \mathbb{N}}$ be the sequences generated by Algorithm FLEX-BC-PG. The following assertions hold.

1. Sufficient decrease property. The sequence ${\Psi(\bar{\mathbf{x}}^k)}_{k \in \mathbb{N}}$ is non-increasing. For each $n$, let $\beta_f^n = \sqrt{\sum_{j \in \mathbb{I}^n, 1 \le \ell \le L} \varepsilon_j \beta_{\ell,j}^2}$ and $0 < \tau_{\ell}^n < 1/\beta_f^n$. Then, for all $k \ge 0$,

$$ \Psi(\bar{\mathbf{x}}^{k+1}) + \left( \sum_{n=k \cdot K}^{(k+1) \cdot K-1} \sum_{\ell \in \mathbb{I}^n} \frac{1}{2} \left( \frac{1}{\tau_{\ell}^n} - \beta_f^n \right) \|x_{\ell}^n - x_{\ell}^{n+1}\|^2 \right) \le \Psi(\bar{\mathbf{x}}^k). $$

Furthermore,

$$ \sum_{n=0}^{+\infty} \left( \sum_{\ell=1}^{L} \|x_{\ell}^n - x_{\ell}^{n+1}\|^2 \right) < +\infty, $$

which implies $\lim_{n \to +\infty} \|x_{\ell}^n - x_{\ell}^{n+1}\| = 0$ for every $\ell \in \{1, \dots, L\}$, which in turn yields $\lim_{n \to +\infty} \|\mathbf{x}^n - \mathbf{x}^{n+1}\| = 0$.

2. Subgradient bound. For each $k \in \mathbb{N}$ define

$$ \bar{B}^{k+1} = \left( \frac{x_{\ell}^{k\ell} - x_{\ell}^{k\ell+1}}{\tau_{\ell,k\ell}} - \nabla_{\ell} f(\mathbf{x}^{k\ell}) + \nabla_{\ell} f(\bar{\mathbf{x}}^{k+1}) \right)_{1 \le \ell \le L}, $$

where $k_{\ell}$ is a positive integer such that $k \cdot K \le k_{\ell} \le (k+1) \cdot K - 1$, and is the last iteration of cycle $k$ at which block $\ell$ is updated. Then $\bar{B}^{k+1} \in \partial \Psi(\bar{\mathbf{x}}^{k+1})$ and there exist positive numbers $\tau_k$ such that:

$$ \|\bar{B}^{k+1}\| \le \left( \frac{1}{\tau_k} + \beta_f \right) \left( \sum_{n=k \cdot K}^{(k+1) \cdot K-1} L \sum_{\ell \in \mathbb{I}^n} \frac{1}{2} \left( \frac{1}{\tau_{\ell}^n} - \beta_f^n \right) \|x_{\ell}^n - x_{\ell}^{n+1}\| \right). \quad (10) $$

*Proof. 1.* First, for all $\ell \notin \mathbb{I}^n$, $x_{\ell}^{n+1} = x_{\ell}^{n}$. For all $\ell \in \mathbb{I}^n$, by applying the first order optimality conditions of the proximity operator (2) we obtain:

$$ g_{\ell}(x_{\ell}^{n+1}) + \frac{1}{2\tau_{\ell}^{n}} \|x_{\ell}^{n} - x_{\ell}^{n+1}\|^{2} \le g(x_{\ell}^{n}) + \langle \nabla_{\ell} f(\mathbf{x}^{n}), x_{\ell}^{n} - x_{\ell}^{n+1} \rangle, $$
(see Appendix A, Lemma A.2), which we can sum up for all $\ell \in \mathbb{I}^n$ to obtain

$$ \sum_{\ell \in \mathbb{I}^n} \left( g_{\ell}(x_{\ell}^{n+1}) + \frac{1}{2\tau_{\ell}^n} \|x_{\ell}^n - x_{\ell}^{n+1}\|^2 \right) \leq \sum_{\ell \in \mathbb{I}^n} \left( g(x_{\ell}^n) + \langle \nabla_{\ell} f(\mathbf{x}^n), x_{\ell}^n - x_{\ell}^{n+1} \rangle \right). \quad (11) $$

We now invoke A4 from Assumption 2 and by splitting the scalar product along the blocks we get

$$ f(\mathbf{x}^n + [x_1^{n+1} - x_1^n, \dots, x_L^{n+1} - x_L^n]^\top) \leq f(\mathbf{x}^n) + \sum_{\ell \in \mathbb{I}^n} \left( \langle \nabla_{\ell} f(\mathbf{x}^n), x_{\ell}^{n+1} - x_{\ell}^n \rangle + \frac{\beta_f^n}{2} \|x_{\ell}^n - x_{\ell}^{n+1}\|^2 \right). $$

Note that $f(\mathbf{x}^n + [x_1^{n+1} - x_1^n, \dots, x_L^{n+1} - x_L^n]^\top) = f(\mathbf{x}^{n+1})$, and thus

$$ \sum_{\ell \in \mathbb{I}^n} \left( \langle \nabla_{\ell} f(\mathbf{x}^n), x_{\ell}^n - x_{\ell}^{n+1} \rangle \right) \leq f(\mathbf{x}^n) - f(\mathbf{x}^{n+1}) + \frac{\beta_f^n}{2} \sum_{\ell \in \mathbb{I}^n} \|x_{\ell}^n - x_{\ell}^{n+1}\|^2. \quad (12) $$

Combining inequalities (11) and (12) we obtain

$$ \sum_{\ell \in \mathbb{I}^n} \left( g_{\ell}(x_{\ell}^{n+1}) + \frac{1}{2\tau_{\ell}^n} \|x_{\ell}^n - x_{\ell}^{n+1}\|^2 \right) \leq \sum_{\ell \in \mathbb{I}^n} \left( g(x_{\ell}^n) + f(\mathbf{x}^n) - f(\mathbf{x}^{n+1}) + \frac{\beta_f^n}{2} \sum_{\ell \in \mathbb{I}^n} \|x_{\ell}^n - x_{\ell}^{n+1}\|^2 \right). $$

We add $\sum_{\ell \notin \mathbb{I}^n} g_{\ell}(x_{\ell}^n)$ to each side of this inequality, and since $\sum_{\ell \notin \mathbb{I}^n} g_{\ell}(x_{\ell}^n) = \sum_{\ell \notin \mathbb{I}^n} g_{\ell}(x_{\ell}^{n+1})$ (these blocks are not updated), we have

$$ \Psi(\mathbf{x}^{n+1}) + \sum_{\ell \in \mathbb{I}^n} \frac{1}{2} \left( \frac{1}{\tau_{\ell}^n} - \beta_f^n \right) \|x_{\ell}^n - x_{\ell}^{n+1}\|^2 \leq \Psi(\mathbf{x}^n). \quad (13) $$

Now bundle K block iterations together to define the sequence $(\bar{\mathbf{x}}^k)_{k \in \mathbb{N}}$. Let $k \in \mathbb{N}$, summing inequality (13) on all iterations from $n = k \cdot K$ to $n = (k+1) \cdot K - 1$, and taking into account that $\mathbf{x}^{k \cdot K} = \bar{\mathbf{x}}^k$ and $\mathbf{x}^{(k+1) \cdot K} = \bar{\mathbf{x}}^{k+1}$, we obtain

$$ \Psi(\bar{\mathbf{x}}^{k+1}) + \sum_{n=k \cdot K}^{(k+1) \cdot K-1} \sum_{\ell \in \mathbb{I}^n} \frac{1}{2} \left( \frac{1}{\tau_{\ell}^n} - \beta_f^n \right) \|x_{\ell}^n - x_{\ell}^{n+1}\|^2 \leq \Psi(\bar{\mathbf{x}}^k), $$

and the result follows.
Let us now take $n_0 \in \mathbb{N} \setminus \{0\}$. Summing up inequality (13) from $n=0$ to $n_0 - 1$, we obtain

$$ \sum_{n=0}^{n_0-1} \left( \sum_{\ell \in \mathbb{I}^n} \frac{1}{2} \left( \frac{1}{\tau_\ell^n} - \beta_f^n \right) \|x_\ell^n - x_\ell^{n+1}\|^2 \right) \le \Psi(\mathbf{x}^0) - \Psi(\mathbf{x}^{n_0}) \\ \le \Psi(\mathbf{x}^0) - \inf \Psi $$

and then, by setting $C := \min_{n=0,...,n_0-1} \frac{1}{2} \left( \frac{1}{\tau_\ell^n} - \beta_f^n \right) > 0$ we get

$$ \sum_{n=0}^{n_0-1} \left( \sum_{\ell \in \mathbb{I}^n} \|x_\ell^n - x_\ell^{n+1}\|^2 \right) \le \frac{1}{C} (\Psi(\mathbf{x}^0) - \inf \Psi) < +\infty. $$

The result follows by taking $n_0 \to +\infty$.

Now, set

$$ (\forall k \in \mathbb{N}), \quad \rho_k = \min_{k \cdot K \le n \le (k+1) \cdot K-1} \left( \frac{1}{\tau_\ell^n} - \frac{\beta_f^n}{2} \right) \ge 0, $$

and set

$$ C_k = \sum_{n=k \cdot K}^{(k+1) \cdot K-1} \sum_{\ell \in \mathbb{I}^n} \|x_\ell^n - x_\ell^{n+1}\|^2. \quad (14) $$

2. For all $k \in \mathbb{N}$ and for all $\ell \in \{1, \dots, L\}$, there exists an iteration index $k_\ell \in \mathbb{N} \cap [k \cdot K, (k+1) \cdot K - 1]$ in which block $\ell$ receives its last update in cycle $k$. Thus for all $\ell \in \{1, \dots, L\}$ from (3)

$$ \frac{x_{\ell}^{k_{\ell}} - x_{\ell}^{k_{\ell+1}}}{\tau_{\ell}^{k_{\ell}}} - \nabla_{\ell} f(\mathbf{x}^{k_{\ell}}) \in \partial g_{\ell}(x_{\ell}^{k_{\ell+1}}), $$

which yields

$$ \bar{B}_{\ell}^{k+1} := \frac{x_{\ell}^{k_{\ell}} - x_{\ell}^{k_{\ell+1}}}{\tau_{\ell}^{k_{\ell}}} - \nabla_{\ell} f(\mathbf{x}^{k_{\ell}}) + \nabla_{\ell} f(\bar{\mathbf{x}}^{k+1}) \\ \in \partial g_{\ell}(x_{\ell}^{k_{\ell+1}}) + \nabla_{\ell} f(\bar{\mathbf{x}}^{k+1}) = \partial_{\ell} \Psi(\bar{\mathbf{x}}^{k+1}). $$

where the last equality follows from Proposition 3.1. Hence, by setting $\bar{B}^{k+1} = (\bar{B}_{\ell}^{k+1})_{1 \le \ell \le L}$, we have $\bar{B}^{k+1} \in \partial \Psi(\bar{\mathbf{x}}^{k+1})$. An upper bound on the norm of $\bar{B}^{k+1}$ follows from the block-smoothness property of $f$. We have for all $\ell \in \{1, \dots, L\}$

$$ \| \nabla_{\ell} f(\bar{\mathbf{x}}^{k+1}) - \nabla_{\ell} f(\mathbf{x}^{k_{\ell}}) \| \le \sum_{j=1}^{L} \beta_{\ell,j} \| \bar{x}_{j}^{k+1} - x_{j}^{k_{\ell}} \| \\ = \sum_{j=1}^{L} \beta_{\ell,j} \left\| \sum_{n=k_{\ell}}^{(k+1) \cdot K-1} x_{j}^{n+1} - x_{j}^{n} \right\| $$
$$
\begin{align*}
&\le \sum_{j=1}^{L} \beta_{\ell,j} \sum_{n=k_\ell}^{(k+1)\cdot K-1} \|x_j^{n+1} - x_j^n\| \\
&\le \left( \max_{1\le\ell,j\le L} \beta_{\ell,j} \right) \sum_{j=1}^{L} \sum_{n=k\cdot K}^{(k+1)\cdot K-1} \|x_j^{n+1} - x_j^n\| \\
&\le \left( \max_{1\le\ell,j\le L} \beta_{\ell,j} \right) \sum_{n=k\cdot K}^{(k+1)\cdot K-1} \sum_{j\in\mathbb{I}^n} \|x_j^{n+1} - x_j^n\|.
\end{align*}
$$

Now define

$$
\nu_k = L \left( \frac{1}{\min_{k \cdot K \le n \le (k+1) \cdot K-1} \tau_{\ell}^n} + \max_{1 \le \ell, j \le L} \beta_{\ell, j} \right),
$$

and

$$
D_k = \sum_{n=k \cdot K}^{(k+1) \cdot K - 1} \sum_{j \in \mathbb{I}^n} \|x_j^{n+1} - x_j^n\|.
$$

Thus, for all $k \in \mathbb{N}$, there exists an element $\bar{B}^{k+1} \in \partial\Psi(\mathbf{x}^{k+1})$ whose norm is upper bounded by:

$$
\|\bar{B}^{k+1}\| \leq \nu_k D_k.
$$

□

The proof of the convergence of the sequence requires the study of the limit points
set, defined as follows.

**Definition 4.** *Limit points set [30].* The set of all limit points of sequences generated by FLEX – BC – PG from a starting point $\mathbf{x}^0 = \bar{\mathbf{x}}^0$ will be denoted by $\mathrm{lp}(\bar{\mathbf{x}}^0)$:

$$
\mathrm{lp}(\bar{\mathbf{x}}^0) = \{ \hat{\mathbf{x}} \in \mathcal{H} \,\exists\ \text{an increasing sequence of integers } \{k_j\}_{j \in \mathbb{N}}, \\
\text{such that } \bar{\mathbf{x}}^{k_j} \to \hat{\mathbf{x}} \text{ as } j \to +\infty \}
$$

The properties of the limit points of sequences produced by some block algorithms
are investigated in [30], small modifications are required in our context.

**Lemma 3.8.** *Properties of the limit points set.* Suppose that Assumptions 1 and 2 hold. Let {x̄^k}k∈N be a sequence generated by Algorithm FLEX – BC – PG starting from x̄^0 = x^0. The following hold:

1. $\emptyset \neq \mathrm{lp}(\bar{\mathbf{x}}^0) \subset \mathrm{crit} \ \Psi$.

2. We have

$$
\lim_{k \to \infty} \mathrm{dist}(\bar{\mathbf{x}}^k, \mathrm{lp}(\bar{\mathbf{x}}^0)) = 0.
$$

3. $\mathrm{lp}(\bar{\mathbf{x}}^0)$ is a nonempty, compact and connected set.

4. The objective function $\Psi$ is finite and constant on $\mathrm{lp}(\bar{\mathbf{x}}^0)$.

Proof. 1.Let $\hat{\mathbf{x}}$ be a limit point of $\{\bar{\mathbf{x}}^k\}_{k \in \mathbb{N}}$. By Definition 4 there exists a subsequence $\{\bar{\mathbf{x}}^{k_q}\}_{q \in \mathbb{N}}$ such that $\bar{\mathbf{x}}^{k_q} \to \hat{\mathbf{x}}$. Using assumption A3, it follows that

$$
(\forall \ell \in \{1, \dots, L\}), \quad \liminf_{q \to +\infty} g_{\ell}(x_{\ell}^{k_q}) \geq g_{\ell}(\hat{x}_{\ell}).
$$
It is enough to prove that

$$ (\forall \ell \in \{1, \dots, L\}), \quad \limsup_{q \to +\infty} g_{\ell}(x_{\ell}^{k_q}) \le g_{\ell}(\hat{x}_{\ell}). \qquad (15) $$

Indeed, since $f$ is continuous, (15) implies that

$$ \lim_{q \to \infty} \Psi(\bar{\mathbf{x}}^{k_q}) = \Psi(\hat{\mathbf{x}}). $$

For all $k \in \mathbb{N}$ and for all $\ell \in \{1, \dots, L\}$, we denote by $k_\ell \in \mathbb{N}$ the iteration in which block $\ell$ has received its last update in the $k$-th cycle. We have

$$ (\forall \ell \in \{1, \dots, L\}) \quad x_{\ell}^{k+1} = \arg \min_{x_{\ell} \in \mathcal{H}_{\ell}} \left\{ \langle x_{\ell} - x_{\ell}^{k_{\ell}}, \nabla_{\ell} f(\mathbf{x}^{k_{\ell}}) \rangle + \frac{1}{2\tau_{\ell}^{k_{\ell}}} \|x_{\ell} - x_{\ell}^{k_{\ell}}\|^2 + g_{\ell}(x_{\ell}) \right\}. $$

Thus, for $x_\ell = \hat{x}_\ell$ it holds

$$ \langle x_{\ell}^{k} - x_{\ell}^{k_{\ell}}, \nabla_{\ell} f(\mathbf{x}^{k_{\ell}}) \rangle + \frac{1}{2\tau_{\ell}^{k_{\ell}}} \|x_{\ell}^{k} - x_{\ell}^{k_{\ell}}\|^2 + g_{\ell}(x_{\ell}^{k}) \le \langle \hat{x}_{\ell} - x_{\ell}^{k_{\ell}}, \nabla_{\ell} f(\mathbf{x}^{k_{\ell}}) \rangle \\ + \frac{1}{2\tau_{\ell}^{k_{\ell}}} \|\hat{x}_{\ell} - x_{\ell}^{k_{\ell}}\|^2 + g_{\ell}(\hat{x}_{\ell}). $$

The index $k_\ell$ depends implicitly on $k$. Let $q \in \mathbb{N}$. For the rest of the proof we need to extract a converging subsequence, and to note the dependence to $q$ we will write $k_{q,\ell}$ to indicate the last update received by block $\ell$ at cycle $k_q$. Taking then $k = k_q$, we obtain

$$ \langle x_{\ell}^{k_q} - x_{\ell}^{k_{q,\ell}}, \nabla_{\ell} f(\mathbf{x}^{k_{q,\ell}}) \rangle + \frac{1}{2\tau_{\ell}^{k_{q,\ell}}} \|x_{\ell}^{k_q} - x_{\ell}^{k_{q,\ell}}\|^2 + g_{\ell}(x_{\ell}^{k_q}) \\ \le \langle \hat{x}_{\ell} - x_{\ell}^{k_{q,\ell}}, \nabla_{\ell} f(\mathbf{x}^{k_{q,\ell}}) \rangle + \frac{1}{2\tau_{\ell}^{k_{q,\ell}}} \|\hat{x}_{\ell} - x_{\ell}^{k_{q,\ell}}\|^2 + g_{\ell}(\hat{x}_{\ell}) \quad (16) $$

Now we look at the limit when $q$ goes to infinity. Using the following properties:

- $\|x_{\ell}^{k_q} - x_{\ell}^{k_{q,\ell}}\|$ goes to 0 as $q$ goes to infinity (Proposition 3.7, item 1.),

- $\nabla_{\ell} f$ is Lipschitz continuous and the sequence $(\mathbf{x}^n)_{n \in \mathbb{N}}$ is bounded (Assumption A4),

- $\|\hat{x}_{\ell} - x_{\ell}^{k_{q,\ell}}\| \le \|\hat{x}_{\ell} - x_{\ell}^{k_q}\| + \|x_{\ell}^{k_q} - x_{\ell}^{k_{q,\ell}}\|$ and both terms on the right-hand side of the inequality go to 0 as $q$ goes to infinity.

Hence, (16) yields

$$ \limsup_{q \to +\infty} g_{\ell}(x_{\ell}^{k_q}) \le g_{\ell}(\hat{x}_{\ell}) $$

Now, combining points (1) and (2) of Proposition 3.7 that $\bar{B}^k \to 0$ as $k \to +\infty$. The closedness property of $\partial\Psi$ [58, Theorem 8.6] implies that $0 \in \partial\Psi(\hat{\mathbf{x}})$, and therefore that $\hat{\mathbf{x}}$ is a critical point of $\Psi$. We have from Proposition 3.7(1) and Lemma 3.6 that
$||\bar{\mathbf{x}}^{k+1} - \bar{\mathbf{x}}^k|| \to 0$ as $k \to +\infty$, and thus 2. and 3. hold (see [30, Remark 5 & Lemma 5]). 4. is derived from 1. ([30, Lemma 5]). □

Now that we have established the decrease of the objective function at each iteration and the properties of the limit point set, we are ready to state our main result.

**Theorem 3.9.** Suppose that Assumptions 1, 2, and 3 hold. Let {$\mathbf{x}^n$}$_{n \in \mathbb{N}}$, and {$\bar{\mathbf{x}}^k$}$_{k \in \mathbb{N}}$ be the sequences generated by Algorithm FLEX – BC – PG. Then,

1. The sequence {$\bar{\mathbf{x}}^k$}$_{k \in \mathbb{N}}$ has finite length, that is,

$$ \sum_{k=1}^{\infty} |||\bar{\mathbf{x}}^k - \bar{\mathbf{x}}^{k+1}||| < \infty. $$

2. The sequence {$\bar{\mathbf{x}}^k$}$_{k \in \mathbb{N}}$ converges to a critical point $\hat{\mathbf{x}}$ of $\Psi$.

*Proof. 1.* For this crucial result, we follow the path of [30]. Since the sequence $(\bar{\mathbf{x}}^k)_{k \in \mathbb{N}}$ is bounded, there exists a sub-sequence that converges to $\hat{\mathbf{x}}$.

* As {$\Psi(\bar{\mathbf{x}}^k)$}$_{k \in \mathbb{N}}$ is a non-increasing sequence, and as the limit points set $\text{lp}(\mathbf{x}^0)$ is such that $\lim_{k \to \infty} \text{dist}(\bar{\mathbf{x}}^k, \text{lp}(\mathbf{x}^0)) = 0$ (Lemma 3.8, point (ii)), there exist $k_0 \in \mathbb{N}, \varepsilon > 0, \eta > 0$ such that for all $k > k_0$, $\bar{\mathbf{x}}^k$ belongs to:

$$ \{\mathbf{x} \in \mathbb{R}^N : \text{dist}(\mathbf{x}, \text{lp}(\mathbf{x}^0)) < \varepsilon\} \cap [\Psi(\hat{\mathbf{x}}) < \Psi(\mathbf{x}) < \Psi(\hat{\mathbf{x}}) + \eta]. $$

* Using now that $\text{lp}(\mathbf{x}^0)$ is nonempty and compact, and that $\Psi$ is constant on it (Lemma 3.8 (ii), (iii) and (iv)), by applying Lemma 3.2, we obtain

$$ (\forall k > k_0), \quad \varphi'(\Psi(\bar{\mathbf{x}}^k) - \Psi(\hat{\mathbf{x}})) \text{dist}(0, \partial\Psi(\bar{\mathbf{x}}^k)) \ge 1, \quad \text{for some } \varphi \in \Phi_{\eta} \quad (17) $$

* Let $k \in \mathbb{N}$, Proposition 3.7 gives an upper bound of $\text{dist}(0, \partial\Psi(\bar{\mathbf{x}}^k))$. Set $\rho_1 = \min_{k \in \mathbb{N}} \rho_k$ and $\rho_2 = \max_{k \in \mathbb{N}} \nu_k$. We have

$$ \text{dist}(0, \partial\Psi(\bar{\mathbf{x}}^k)) \le \rho_2 D_{k-1} $$

Combined with (17), this bound implies

$$ \varphi'(\Psi(\bar{\mathbf{x}}^k) - \Psi(\hat{\mathbf{x}})) \ge \rho_2^{-1} D_{k-1}^{-1}. \quad (18) $$

* The concavity of $\varphi$ yields that:

$$ \varphi(\Psi(\bar{\mathbf{x}}^k) - \Psi(\hat{\mathbf{x}})) - \varphi(\Psi(\bar{\mathbf{x}}^{k+1}) - \Psi(\hat{\mathbf{x}})) \ge \varphi'(\Psi(\bar{\mathbf{x}}^k) - \Psi(\hat{\mathbf{x}})) (\Psi(\bar{\mathbf{x}}^k) - \Psi(\bar{\mathbf{x}}^{k+1})). \quad (19) $$

* Now recall from the proof of Proposition 3.7, item 1., (14) that

$$ \rho_1 C_k \le \Psi(\bar{\mathbf{x}}^k) - \Psi(\bar{\mathbf{x}}^{k+1}). $$
Define

$$ (\forall p, q \in \mathbb{N}), \quad \Delta_{p,q} := \varphi(\Psi(\bar{\mathbf{x}}^p) - \Psi(\bar{\mathbf{x}})) - \varphi(\Psi(\bar{\mathbf{x}}^q) - \Psi(\bar{\mathbf{x}})), $$

By setting $\rho := \rho_1\rho_2^{-1} > 0$, we have from (18) and (19)

$$ \Delta_{k,k+1} \geq \frac{\rho C_k}{D_{k-1}} $$

and then:

$$ C_k \leq \rho^{-1} \Delta_{k,k+1} D_{k-1} $$

• We rewrite the expression of both $C_k$ and $D_k$

$$ D_k = \sum_{j=1}^{m_k} a_{k,j}, \quad \text{and} \quad C_k = \sum_{j=1}^{m_k} a_{k,j}^2. $$

where $m_k = \sum_{n=k \cdot K}^{(k+1)} \text{card}(\mathbb{I}^n)$ and $a_{k,j} = \|x_{\ell}^n - x_{\ell+1}^{n+1}\|$ where $j$ browses through $n$ and $\ell$ by increasing order. Then,

$$ D_k \leq \sqrt{m_k C_k} $$

Using $2\sqrt{ab} \leq a + b$ for all $a, b \geq 0$, and writing $M = \max_k \sqrt{m_k}$ we get that

$$ 2D_k \leq M\rho^{-1}\Delta_{k,k+1} + D_{k-1} $$

We have then

$$ \begin{aligned} 2 \sum_{i=k_0+1}^{k} D_i &\le \sum_{i=k_0+1}^{k} D_{i-1} + M\rho^{-1} \sum_{i=k_0+1}^{k} \Delta_{i,i+1} \\ &= \sum_{i=k_0}^{k-1} D_i + M\rho^{-1} \sum_{i=k_0+1}^{k} \Delta_{i,i+1} \\ &\le \sum_{i=k_0+1}^{k} D_i + D_{k_0} + M\rho^{-1} \sum_{i=k_0+1}^{k} \Delta_{i,i+1} \\ &\Rightarrow \sum_{i=k_0+1}^{k} D_i \le D_{k_0} + M\rho^{-1} \sum_{i=k_0+1}^{k} \Delta_{i,i+1} \\ &\Rightarrow \sum_{i=k_0+1}^{k} D_i \le D_{k_0} + M\rho^{-1} \Delta_{k_0+1,k+1}, \end{aligned} $$
the last line coming from the fact $\Delta_{p,q} + \Delta_{q,r} = \Delta_{p,r}$ for all $p, q, r \in \mathbb{N}$ [30, Proof of theorem 3.1]. As $\varphi \ge 0$, we have that:

$$ \Delta_{k_0+1, k+1} = \varphi(\Psi(\bar{\mathbf{x}}^{k_0+1}) - \Psi(\hat{\mathbf{x}})) - \varphi(\Psi(\bar{\mathbf{x}}^{k+1}) - \Psi(\hat{\mathbf{x}})) \le \varphi(\Psi(\bar{\mathbf{x}}^{k_0+1}) - \Psi(\hat{\mathbf{x}})) $$

Then as $D_k \ge \|\bar{\mathbf{x}}^k - \bar{\mathbf{x}}^{k+1}\|$ (Lemma 3.6), this allows us then to conclude then that $\{\bar{\mathbf{x}}^k\}_{k \in \mathbb{N}}$ has finite length:

$$ \sum_{k=1}^{\infty} \|\bar{\mathbf{x}}^k - \bar{\mathbf{x}}^{k+1}\| < \infty. $$

2. The finite length of the sequence implies that it is a Cauchy sequence and hence a convergent sequence [30, Proof of Theorem 3.1(ii)]. □

We split the proof of our main result into 4 steps, which are common when studying descent algorithm on KL functions [6, 30, 37, 38, 59]. For each step we detail how the existing proofs were adapted for our approach by highlighting the difference with the existing literature.

1. **Sufficient decrease property:** at each cycle, the objective function $\Psi$ is decreased. The decrease is controlled by the squared norm of the differences between the block updates. *Difference with the literature:* we introduce the possibility of parallel updates to decrease the function, and thus to adapt the choice of step size to the smoothness of the group of blocks considered.

2. **Subgradient upper bound:** at each cycle, we can exhibit an upper bound on one element of the subgradient of $\Psi$ at the cycle iterate. This upper bound is controlled by the norm of the differences between the block updates. *Difference with the literature:* this bound is less sharp than in the literature (as the reader can see in the proof), but it is necessary to write it in this way to apply the KL property and obtain finite length.

3. **Limit points are critical points:** the set of limit points of the sequences generated by our algorithm will be a subset of the set of critical points of $\Psi$. *Difference with the literature:* due to the possible parallel block updates, we need to be a bit more cautious when looking at the converging subsequences.

4. **Finite length of the sequences:** the sequences generated by our algorithm have finite length and thus converge. This is a consequence of the KL property satisfied by $\Psi$ (see Lemma 3) and of points (i) and (ii). *Difference with the literature:* we invoke a particular instance of Cauchy-Schwartz inequality to obtain the desired result.

**Remark 3.10.** Regarding Theorem 3.9, it may not appear obvious that algorithm FLEX – BC – PG needs to update each block in an essentially cyclic manner, but this is indeed necessary, otherwise the norm of $\bar{B}^k \in \partial\Psi(\bar{\mathbf{x}}^k)$ would not go to 0 as $k$ goes to infinity.
## 3.5 Some consequences.

Our main result implies that the sequence of iterates {$\mathbf{x}^n$}$_{n \in \mathbb{N}}$ also converges, and to establish some convergence rates depending on the KL property.

### Convergence of the sequence {$\mathbf{x}^n$}$_{n \in \mathbb{N}}$.

We have established the convergence of the sequence of cycle iterates {$\bar{\mathbf{x}}^k$}$_{n \in \mathbb{N}}$ to a critical point while decreasing the objective function. A direct consequence is that the sequence {$\mathbf{x}^n$}$_{n \in \mathbb{N}}$ converges to the same critical point.

**Proposition 3.11.** Suppose that Assumptions 1, 2, and 3 hold. Let {$\mathbf{x}^n$}$_{n \in \mathbb{N}}$, {$\bar{\mathbf{x}}^k$}$_{k \in \mathbb{N}}$ be generated by Algorithm FLEX – BC – PG. Let $\hat{\mathbf{x}} \in \text{crit } \Psi$ such that $\bar{\mathbf{x}}^k \to_{k \to +\infty} \hat{\mathbf{x}}$. Then,

$$ \lim_{n \to +\infty} \mathbf{x}^n = \hat{\mathbf{x}}. $$

*Proof.* With a similar argument than in Lemma 3.6, we have

$$ (\forall k \in \mathbb{N}), \forall n_0 \in \mathbb{N} \cap [k \cdot K, (k+1) \cdot K - 1]), \quad \left\| \mathbf{x}^{n_0} - \bar{\mathbf{x}}^k \right\| \le \sum_{n=k \cdot K}^{n_0-1} \sum_{\ell \in \mathbb{I}^n} \| x_{\ell}^{n+1} - x_{\ell}^n \|. $$

Combining the sufficient decrease property (Proposition 3.7(1)) of our algorithm and the previous bound, we have

$$ (\forall n_0 \in \mathbb{N} \cap [k \cdot K, (k+1) \cdot K - 1]), \quad \lim_{k \to +\infty} \left\| \mathbf{x}^{n_0} - \bar{\mathbf{x}}^k \right\| = 0. $$

We showed that {$\bar{\mathbf{x}}^k$}$_{k \in \mathbb{N}}$ converges to $\hat{\mathbf{x}}$, hence by definition

$$ \forall \xi_1 > 0, \exists k_1 \in \mathbb{N}, \text{ such that for all } k \ge k_1, \left\| \bar{\mathbf{x}}^k - \hat{\mathbf{x}} \right\| < \xi_1. $$

Now choose $\xi_2 > 0$ such that for all $k \ge k_2 := k_1$, $\left\| \mathbf{x}^{n_0} - \bar{\mathbf{x}}^k \right\| < \xi_2$. We have

$$ \begin{aligned} \| \mathbf{x}^{n_0} - \hat{\mathbf{x}} \| &= \| \mathbf{x}^{n_0} - \bar{\mathbf{x}}^k + \bar{\mathbf{x}}^k - \hat{\mathbf{x}} \| \\ &\le \| \mathbf{x}^{n_0} - \bar{\mathbf{x}}^k \| + \| \bar{\mathbf{x}}^k - \hat{\mathbf{x}} \| \\ &\le \xi_1 + \xi_2 \end{aligned} $$

To prove that $\forall \xi > 0$, there exists $N \in \mathbb{N}$, such that for all $n \ge N$, $\| \mathbf{x}^{n_0} - \hat{\mathbf{x}} \| < \xi$, we can take $\xi_1 = \xi_2 := \xi/2$ and take $N = \max\{k_1, k_2\} \cdot K$ to conclude. $\square$

### Convergence rates.

As developed in [30, 37], the case in which all the functions involved are semi-algebraic is an interesting one for the study of the convergence. The desingularizing function $\varphi$ can be chosen to be of the form [37] $\varphi(s) = cs^{1-\theta}$, where $c > 0$ and $\theta \in [0, 1)$. Then, depending on the value of $\theta$, the following convergence rates hold: [30, 37]

1. If $\theta = 0$ the sequence {$\bar{\mathbf{x}}^k$}$_{k \in \mathbb{N}}$ converges in a finite number of steps.
2. If $\theta \in (0, 1/2]$ there exist $\omega > 0$ and $\kappa \in [0, 1)$ such that $\|\bar{\mathbf{x}}^k - \hat{\mathbf{x}}\| \le \omega\kappa^k$.

3. If $\theta \in (1/2, 1)$ there exists $\omega > 0$ such that $\|\bar{\mathbf{x}}^k - \hat{\mathbf{x}}\| \le \omega k^{-\frac{1-\theta}{2\theta-1}}$.

## 3.6 Convexity of the regularization

We can derive a slightly different sufficient decrease property of our algorithm when assuming convexity of the regularizing functions $g_\ell$, for all $\ell$. This assumption allows us to take bigger step sizes when updating the blocks.

**Lemma 3.12.** *Sufficient decrease property: convexity of $g_\ell$.* Suppose that Assumption 2 holds, and that for all $\ell \in \{1, \dots, L\}$, $g_\ell$ is convex. Let $\{\bar{\mathbf{x}}^k\}_{k \in \mathbb{N}}$, and $\{\mathbf{x}^n\}_{n \in \mathbb{N}}$ be generated by algorithm FLEX-BC-PG. For each $n$, let $\beta_f^n = \sqrt{\sum_{j \in J_n, 1 \le \ell \le L} \varepsilon_j \beta_{\ell,j}^2}$ and $0 < \tau_\ell^n < 2/\beta_f^n$. Then

$$ \Psi(\bar{\mathbf{x}}^{k+1}) + \left( \sum_{n=k \cdot K}^{(k+1) \cdot K-1} \sum_{\ell \in \mathbb{I}^n} \left( \frac{1}{\tau_\ell^n} - \frac{\beta_f^n}{2} \right) \|x_\ell^n - x_\ell^{n+1}\|^2 \right) \le \Psi(\bar{\mathbf{x}}^k). $$

*Proof.* Due to the convexity of $g_\ell$ for all $\ell$, the associated proximity operator is single-valued. For all $\ell \in \mathbb{I}^n$, the first order optimality conditions of the proximity operator (2) yield:

$$ g_{\ell}(x_{\ell}^{n+1}) + \frac{1}{\tau_{\ell}^{n}} \|x_{\ell}^{n} - x_{\ell}^{n+1}\|^{2} \leq g(x_{\ell}^{n}) + \langle \nabla_{\ell} f(\mathbf{x}^{n}), x_{\ell}^{n} - x_{\ell}^{n+1} \rangle. \quad (20) $$

The subtle difference with the non-convex case is the factor dividing $\|x_\ell^n - x_\ell^{n+1}\|^2$. With the same derivation as in the non-convex case (proof of point (i) in Theorem 3.9), we obtain finally that:

$$ \Psi(\bar{\mathbf{x}}^{k+1}) + \left( \sum_{n=k \cdot K}^{(k+1) \cdot K-1} \sum_{\ell \in \mathbb{I}^n} \left( \frac{1}{\tau_\ell^n} - \frac{\beta_f^n}{2} \right) \|x_\ell^n - x_\ell^{n+1}\|^2 \right) \le \Psi(\bar{\mathbf{x}}^k). $$

□

# 4 Multilevel forward-backward as a flexible block-coordinate forward-backward

We have seen how to construct a convergent parallel and essentially cyclic block-coordinate forward-backward algorithm, able to handle non-convexity and non-smoothness of the objective function. This algorithm is deterministic by essence, even though a random shuffle of the order of the updates is possible.

Our algorithm is compatible with several update rules. In particular, FLEX-BC-PG can use a hierarchical update rule, to update blocks given their position in a hierarchy inherited from the function to optimize. Such hierarchy arises
for instance in image restoration problems when the regularization penalizes wavelet coefficients of the image. Such hierarchical update rule has a direct connection with multilevel algorithms.

We start by recalling some key facts about multiresolution analysis to introduce the rigorous use of wavelets made in this chapter; then present the construction of the block algorithm for two blocks; and then discuss the construction of the multilevel algorithm for two levels. This presentation will allow us to highlight that the two approaches can be rigorously equivalent in this setting. This equivalence opens practical applications for BCD methods and theoretical insights for multilevel methods.

## 4.1 Multiresolution analysis in a nutshell

Let $\mathbf{u} \in L_2(\Omega)$ be an image, where $\Omega \subset \mathbb{R}^2$. We decompose $L_2(\Omega)$ into the sum of the space of approximation coefficients $V_J$, and that of detail coefficients $V_J^\perp$ at resolution $J \in \mathbb{N}$ [54, Chapter 7], i.e.,

$$L_2(\Omega) = V_J \oplus V_J^\perp.$$

We assume that $\mathbf{u}$ lives exclusively in $V_J$ in the following, i.e., $\mathbf{u}$ is a discrete image of $2^{2^J}$ pixels. $V_J$ can be decomposed into subspaces $V_{J-1}$ and $W_{J-1}$, where $V_{J-1}$ is the space of approximation coefficients at resolution $J-1$ and $W_{J-1}$ the space of detail coefficients at resolution $J-1$. More precisely,

$$\mathbf{u} = \Pi_{V_{J-1}}^* a_{J-1} + \Pi_{W_{J-1}}^* d_{J-1},$$

where $a_{J-1}$ (resp. $d_{J-1}$) are the approximation (resp. detail) coefficients at resolution $J-1$ and $\Pi_{V_{J-1}}$ (resp. $\Pi_{W_{J-1}}$) is the linear projection operator onto $V_{J-1}$ (resp. $W_{J-1}$). By definition, $\Pi_{V_{J-1}}\Pi_{V_{J-1}}^* = \mathrm{Id}_{V_{J-1}}$ and $\Pi_{W_{J-1}}\Pi_{W_{J-1}}^* = \mathrm{Id}_{W_{J-1}}$, where $\mathrm{Id}_{V_{J-1}}$ and $\mathrm{Id}_{W_{J-1}}$ are the identity operators on $V_{J-1}$ and $W_{J-1}$, respectively. Note that the block $d$ contains the three groups of detail coefficients [54].

## 4.2 Wavelet deblurring with multilevel FB: a hierarchical block algorithm

Our $2^{2^J}$ pixels image $\mathbf{u}$ is decomposed into two independent components $a_{J-1} \in V_{J-1}$, and $d_{J-1} \in W_{J-1}$. In this setting the blocks $(x_1, x_2)$ are equal to $(a_{J-1}, d_{J-1})$. We drop the index $J-1$ in the following for simplicity, so that $(x_1, x_2) = (a, d) \in V \times W$. We aim at solving the following optimization problem:

$$\underset{\mathbf{u} \in V_J}{\operatorname{Argmin}} F(\mathbf{u}) = \frac{1}{2} \|\mathbf{A}\mathbf{u} - \mathbf{z}\|_2^2 + \lambda g(\mathbf{D}\mathbf{u}), \quad (21)$$

where $A: V_J \to V_J$ is a bounded linear operator modelling a blur, $\mathbf{z}$ is a degraded version of $\mathbf{u}$, $D$ is the wavelet transform of $\mathbf{u}$ on two levels, and $\lambda$ is a multi-valued parameter to penalize differently the approximation and detail coefficients. We also assume that $g$ is separable along $a$ ($g_a$) and $d$ ($g_d$), and proximable, hence $g \circ D$ is also proximable. We can rewrite this classical wavelet penalized least-squares problem
<figure><img src="image_4.png" /><figcaption>**Fig. 2:** Decomposition on two levels of an image **u** with a wavelet transform. We regroup the detail coefficients $d_1, d_2$ and $d_3$ into one single block **d** to simplify the presentation of our two level or two block proximal gradient descent algorithms.</figcaption></figure>

using the wavelet decomposition of **u**, i.e., the synthesis formulation

$$ (\hat{a}, \hat{d}) \in \underset{a \in V, d \in W}{\operatorname{Argmin}} \Psi(a, d) = \frac{1}{2} \|A (\Pi_V^* a + \Pi_W^* d) - z\|^2 + \lambda_a g_a(a) + \lambda_d g_d(d). \quad (22) $$

To find a minimizer of $\Psi$ with respect to $(a,d)$ is equivalent to finding a minimizer of $F$ with respect to $\mathbf{u}$, as we can recover the solution $\hat{\mathbf{u}}$ of Problem (21) from the solution $(\hat{a}, \hat{d})$ of Problem (22), by setting

$$ \hat{\mathbf{u}} = \Pi_V^* \hat{a} + \Pi_W^* \hat{d}. $$

In the following, note that $D = \begin{bmatrix} \Pi_V \\ \Pi_W \end{bmatrix}$.

**Two block-coordinate proximal gradient descent.**

To minimize $\Psi$, we first consider the block-coordinate approach. Given $(a^0, d^0) \in V \times W$, and $0 < \tau < 2/\|A^*A\|$, the iterations read

$$ \begin{cases} a^{n+1} = a^n + \varepsilon_{a,n} (\operatorname{prox}_{\tau \lambda_a g_a} (a^n - \tau \Pi_V A^* [A (\Pi_V^* a^n + \Pi_W^* d^n) - \mathbf{z}] - a^n)) \\ d^{n+1} = d^n + \varepsilon_{d,n} (\operatorname{prox}_{\tau \lambda_d g_d} (d^n - \tau \Pi_W A^* [A (\Pi_V^* a^n + \Pi_W^* d^n) - \mathbf{z}] - d^n)) \end{cases} \quad (23) $$

where for every $n \in \mathbb{N}$, $(\varepsilon_{a,n}, \varepsilon_{d,n}) \in \{0, 1\}^2$. Note that if $(\varepsilon_{a,n}, \varepsilon_{d,n}) = (1, 1)$ for all $n$, (23) reduces to the standard forward-backward algorithm. A cyclic coordinate descent algorithm is obtained by setting alternatively one of $\varepsilon_{a,n}, \varepsilon_{d,n}$ to 1. This could also be set at random, provided that $P[(\varepsilon_{a,n}, \varepsilon_{d,n}) = (0, 0)] = 0$.

**Two-level proximal gradient descent.**

We present now the construction of a two-level proximal algorithm to minimize $\Psi$ (the fine level function in the multilevel terminology [50]). We will denote by $\Psi_H$ an approximation of $\Psi$, deemed the coarse level function. Given the structure of the
problem, it is natural to define $\Psi_H$ in the approximation space $V$. Consequently, and following [50], the information transfer operator $I_h^H$ (that sends information from the fine level to the coarse level) is the restriction $R_V$ onto $V$ (i.e., $R_V(a, d) = a$) and the prolongation operator $I_h^V$ (that sends information from the coarse level to the fine level) is directly $R_V^*(a, 0)$.

By setting $A_H = A\Pi_V^*$, the coarse model is chosen as:

$$ \Psi_H(a) = \frac{1}{2} \|A_H a - \Pi_V^* \Pi_V z\|_2^2 + \lambda_a g_a(a) + \langle v_H, a \rangle, \quad (24) $$

where $v_H$ enforces the first order coherence [50, Definition 2.1] between two smoothed versions $\Psi_\mu$ and $\Psi_{H,\mu}$ of $\Psi$ and $\Psi_H$ respectively [50, Definition 2.4], with parameter $\mu > 0$ (i.e., the smoothed coarse level function to be the first order Taylor approximation of the smoothed fine level function):

$$ v_H = R_V \nabla \Psi_\mu(a^n, d^n) - \nabla \Psi_{H,\mu}(a^n), $$

with $\Psi_{H,\mu}(\cdot) = \|A_H \cdot - \Pi_V^* \Pi_V z\|_2^2 + g_{a,\mu}(\cdot)$, and $\Psi_\mu(\cdot) = \|A \cdot - z\|_2^2 + g_{a,\mu}(\cdot) + g_{d,\mu}(\cdot)$, where we have denoted $g_{a,\mu}$ the $\mu > 0$-smoothed $\ell_1$-norm (according to the principles of [60]). To go from one level to the other, multilevel algorithm employ information transfer operators.

In the following, we assume that we compute only one coarse iteration before going back to the fine level, but everything holds trivially for more coarse iterations. This iteration will yield $a^{n+1/2}$ from $a^n$. The coarse level model being non-smooth, we will use a proximal gradient step to optimize it. Accordingly, the two-level proximal gradient algorithm, starting from $(a_0, d_0) \in V \times W$, is

$$ \left\{ \begin{array}{l} a^{n+1/2} = \operatorname{prox}_{\tau \lambda_a g_a} (a^n - \tau A_H^* (A_H a^n - \Pi_V^* \Pi_V z) - \tau v_H) \\ a^{n+1} = \operatorname{prox}_{\tau \lambda_a g_a} (a^{n+1/2} - \tau \Pi_V A^* (A (\Pi_V^* a^{n+1/2} + \Pi_W^* d^n) - z)) \\ d^{n+1} = \operatorname{prox}_{\tau \lambda_d g_d} (d^n - \tau \Pi_W A^* (A (\Pi_V^* a^{n+1/2} + \Pi_W^* d^n) - z)) \end{array} \right. \quad (25) $$

The fact that Algorithms (23) and (25) are the same algorithm is not obvious at first sight. We will show in the following that this is indeed the case. After summarizing our assumptions, we will compute the first order coherence term explicitly.

**Assumption 4.** We assume that:

(i) the information transfer operator is the restriction $R_V$ onto $V$;

(ii) in the definition of $v_H$, the fine and coarse models are smoothed with the same smoothing technique, with the same smoothing parameter $\mu > 0$;

(iii) $\Psi$ and $\Psi_H$ are first order coherent with respect to their smoothed versions [50, Definition 2.1].

**Lemma 4.1.** Suppose that Assumption 4 holds. The first order coherence term $v_H$ in (24) at $(a^n, d^n)$ is given by:

$$ v_H = \Pi_V A^* (A \Pi_W^* d^n - \Pi_W^* \Pi_W z). \quad (26) $$
The first order coherence sends the contribution of the detail coefficients to the gradient of the data fidelity term to the coarse level.

*Proof.* By definition of first order coherence between smoothed functions [50, Definition 2.1], we have:

$$v_H = R_V \nabla \Psi_{\mu}(a^n, d^n) - \nabla \Psi_{H,\mu}(a^n). \quad (27)$$

The second term is the simple computation of the gradient of the coarse model:

$$\begin{aligned}
\nabla \Psi_{H,\mu}(a^n) &= \nabla \left( \frac{1}{2} \|A_H a^n - \Pi_V^* \Pi_V \mathbf{z}\|_2^2 + \lambda_a g_{a,\mu}(a^n) \right) \\
&= A_H^* (A_H a^n - \Pi_V^* \Pi_V \mathbf{z}) + \lambda_a \nabla_a g_{a,\mu}(a^n) \\
&= \Pi_V A^* (A \Pi_V^* a^n - \Pi_V^* \Pi_V \mathbf{z}) + \lambda_a \nabla_a g_{a,\mu}(a^n).
\end{aligned} \quad (28)$$

On the other hand,

$$\begin{aligned}
\nabla \Psi_{\mu}(a^n, d^n) &= \nabla \left( \frac{1}{2} \|A (\Pi_V^* a^n + \Pi_W^* d^n) - z\|_2^2 + \lambda_a g_{a,\mu}(a^n) + \lambda_d g_{d,\mu}(d^n) \right) \\
&= \left[ \Pi_V A^* (A (\Pi_V^* a^n + \Pi_W^* d^n) - z) + \lambda_a \nabla_a g_{a,\mu}(a^n) \right] \\
&\quad - \left[ \Pi_W A^* (A (\Pi_V^* a^n + \Pi_W^* d^n) - z) + \lambda_d \nabla_d g_{d,\mu}(d^n) \right],
\end{aligned}$$

and thus:

$$v_H = R_V \left( \left[ \Pi_V A^* (A (\Pi_V^* a^n + \Pi_W^* d^n) - \mathbf{z}) + \lambda_a \nabla_a g_{a,\mu}(a^n) \right] \right. \\ \left. - \Pi_V A^* (A \Pi_V^* a^n - \Pi_V^* \Pi_V \mathbf{z}) + \lambda_a \nabla_a g_{a,\mu}(a^n) \right).$$

Since $V$ and $W$ are orthogonal to each other, restricting any element of $W$ on $V$ yields 0 and thus:

$$v_H = \Pi_V A^* (A \Pi_W^* d^n - \Pi_W^* \Pi_W \mathbf{z}),$$

where we used that

$$\mathbf{z} - \Pi_V^* \Pi_V \mathbf{z} = \Pi_W^* \Pi_W \mathbf{z}.$$

□

Based on Lemma 4.1, the proximal gradient step at coarse level at iteration $n$ reads

$$\begin{aligned}
a^{n+1/2} &= \operatorname{prox}_{\tau \lambda_a g_a} (a^n - \tau A_H^* (A_H a^n - \Pi_V^* \Pi_V \mathbf{z}) - \tau v_H) \\
&= \operatorname{prox}_{\tau \lambda_a g_a} (a^n - \tau \Pi_V A^* (A \Pi_V^* a^n - \mathbf{z} + A \Pi_W^* d^n)) \\
&= \operatorname{prox}_{\tau \lambda_a g_a} (a^n - \tau \Pi_V A^* (A (\Pi_V^* a^n + \Pi_W^* d^n) - \mathbf{z})),
\end{aligned} \quad (29)$$

which fits a block-coordinate update on $a$. We summarize the consequence of this result in the following proposition:

**Proposition 4.2.** The two-level algorithm defined in (25) is equivalent to the block-coordinate algorithm defined in (23) when choosing $\varepsilon_{a,2n} = 1$ and $\varepsilon_{d,2n} = 0$, then $\varepsilon_{a,2n+1} = 1$ and $\varepsilon_{d,2n+1} = 1$ for all $n \in \mathbb{N}$ and the same initial vector $(a^0, d^0) \in V \times W$.
<figure><img src="image_5.png" /><figcaption>Fig. 3: Update scheme of the two block-coordinate descent algorithm. The blocks are updated in a cyclic fashion, first with the approximation block updated alone ($a^1$ in red), then the approximation and detail blocks updated together ($a^2$ and $d^2$ in red). We represent two cycles in this figure $k=1$ and $k=2$ for a total of $n=4$ iterations.</figcaption></figure>

A generalization of this argument to an arbitrary number of levels can be found in [61].

Multilevel algorithms have been shown to accelerate the resolution of image restoration problems in the literature [43, 44, 46, 48–51]. Therefore following update rules from the multilevel literature is of great interest for block algorithms (see example (4) in Figure 1 and Figure 3).

# 5 Numerical experiments

In this section, we present numerical experiments to assess the performance of the proposed construction of a block-coordinate descent algorithm mimicking multilevel iterations. We show that by emulating the behavior of multilevel algorithms, our FLEX – BC – PG algorithm has superior practical performance with respect to the other possible choices of update rules. We also show that our algorithm outperforms the standard FB algorithm.

## Optimization problem.

Consider the optimization problem

$$ \underset{\mathbf{u} \in \mathbb{R}^N}{\operatorname{Argmin}} F(\mathbf{u}) = \frac{1}{2} \|A\mathbf{u} - \mathbf{z}\|_2^2 + \lambda g(D\mathbf{u}), \quad (30) $$

where A encodes a Gaussian blur, and D is the 2 levels wavelet decomposition, i.e., $g(D\cdot)$ will penalize the approximation coefficients, and the three blocks of detail coefficients associated with $\mathbf{u}$ (see Figure 2). We set

$$ (\forall (a, d) \in V \times W), \quad g(a, d) = \lambda_a \sum_i \log(|a_i| + \epsilon) + \lambda_d \sum_i \log(|d_i| + \epsilon), $$
where $a_i, d_i$ denote the components of the vectors and $\epsilon > 0$. The proximity operator of $g := \{\lambda_a\} \log(|\cdot| + \epsilon)$ is known explicitly [62] and given by the following:

$$ (\forall a \in \mathbb{R}), \text{prox}_{\tau g}(a) \begin{cases} 0, & \text{if } |a| < 2\sqrt{\tau\lambda_a} - \epsilon \\ \max\left(0, \text{sign}(a) \frac{|a| - \epsilon + \sqrt{(|a| + \epsilon)^2 - 4\tau\lambda_a}}{2}\right), & \text{if } |a| = 2\sqrt{\tau\lambda_a} - \epsilon, \\ \text{sign}(a) \frac{|a| - \epsilon + \sqrt{(|a| + \epsilon)^2 - 4\tau\lambda_a}}{2}, & \text{otherwise} \end{cases} $$

The same formula holds for the detail coefficients with $\lambda_a$ replaced by $\lambda_d$. This problem is not convex, therefore stochastic BC-PG algorithms are not guaranteed to converge to a solution, but for completeness of the presentation, we include them in our experiments.

## Dataset.

In this section we consider the image of the Cameraman, of size 1024 × 1024. We will apply a Gaussian blur and a Gaussian noise to obtain the degraded image. The regularization will be done with a 2-Level log sum-Haar wavelet.

## Experimental setup.

We compare our algorithm to several versions of BC-PG and to the standard forward-backward algorithm. The block methods all consider four blocks from the wavelet decomposition (see Figure 2), with the first block corresponding to the approximation block and the remaining ones to the detail blocks<fnref n="1" />. We will consider the three following algorithms as baselines:

* FB: the forward-backward algorithm. Update rule:

$$ \forall n \in \mathbb{N}, \quad (\varepsilon_0^n, \varepsilon_1^n, \varepsilon_2^n, \varepsilon_3^n) = (1, 1, 1, 1). $$

* Cyclic BC-PG: a cyclic BC-PG algorithm that updates only one block at a time in a cyclic manner. The order of the updates is chosen randomly for one cycle at initialization, and then kept identical for all cycles:

$$ \forall n \in \mathbb{N}, \forall \ell \in \{0, 1, 2, 3\}, \quad \varepsilon_{\ell}^{n} = \begin{cases} 1 & \text{if } \ell = \sigma(n \mod 4) \\ 0 & \text{otherwise} \end{cases}, $$

where $\sigma: \{0, 1, 2, 3\} \rightarrow \{0, 1, 2, 3\}$ is a permutation, and $n \mod 4$ denotes the remainder of the division of $n$ by 4.

* Random BC-PG: the BC-PG algorithm with one randomly chosen block updated at each iteration. Let $(\{0, 1\}, \mathcal{F}, \mathbb{P})$ be a probability space where $\mathcal{F} = 2^{\{0,1\}}$, and the probability measure $\mathbb{P}$ is defined by $\mathbb{P}(\{1\}) = p, \mathbb{P}(\{0\}) = 1 - p$, where $p \in [0, 1]$. Then,

$$ \forall n \in \mathbb{N}, \exists! \ell \in \{0, 1, 2, 3\}, \quad \varepsilon_{\ell}^{n} = 1. $$

The probability of activation is uniform across the blocks.

<footnote n="1">¹Differently from the previous section, the details are not grouped in a single block $d$, they rather form three separate blocks.</footnote>
Note that FB and Cyclic BC-PG are both included in our framework, while Random BC-PG is not.

We will also consider some new schemes, which where not theoretically covered until now, and that allow to exploit the structure of the problem:

* FLEX - BC - PG: the proposed Flexible Block-Coordinate Forward-Backward algorithm, which alternates $m \le 10$ coarse updates on the approximations and $10-m$ full updates. For all $n \in \mathbb{N}$,

$$ (\forall i = 1, \dots, m), (\varepsilon_1^{10n+i}, \varepsilon_2^{10n+i}, \varepsilon_3^{10n+i}, \varepsilon_4^{10n+i}) = (1, 0, 0, 0), $$

$$ (\forall i = m + 1, \dots, 10 - m), (\varepsilon_1^{10n+i}, \varepsilon_2^{10n+i}, \varepsilon_3^{10n+i}, \varepsilon_4^{10n+i}) = (1, 1, 1, 1). $$

* Alternating FLEX - BC - PG: the proposed FLEX - BC - PG algorithm, alternating between $m$ updates on the approximation coefficients and $m-10$ updates on the detail coefficients. For all $n \in \mathbb{N}$,

$$ \forall i = 1, \dots, m, (\varepsilon_1^{10n+i}, \varepsilon_2^{10n+i}, \varepsilon_3^{10n+i}, \varepsilon_4^{10n+i}) = (1, 0, 0, 0) $$

$$ \forall i = m + 1, \dots, 10 - m, (\varepsilon_1^{10n+i}, \varepsilon_2^{10n+i}, \varepsilon_3^{10n+i}, \varepsilon_4^{10n+i}) = (0, 1, 1, 1). $$

* Stochastic FLEX - BC - PG: a stochastic version of the proposed FLEX - BC - PG, with $m \le 10$ coarse iterations followed by $10-m$ full updates in expectation. The proof of convergence of this algorithm is in Appendix B. It is based on a stochastic and parallel BC-PG algorithm, with the selection rule in Lemma B.1. For all $n \in \mathbb{N}$, we impose that $\varepsilon_2^n = \varepsilon_3^n = \varepsilon_4^n := \varepsilon^n$ and

$$ P(\varepsilon_1^n = 1) = 1, $$

$$ P(\varepsilon^n = 1) > 0. $$

In the proposed deterministic hierarchical algorithms, we update all the details coefficients simultaneously. With this choice we intend to show that our block update rule, forcing the update of the approximation coefficients at each iteration, is more efficient than the random one. Note that we also tested a Random BC-PG algorithm that splits the image in four equally sized patches, but the results were worse than all the other algorithms presented here.

The results of this experiment are shown in Figure 4. We can clearly see that FLEX - BC - PG and Alternating FLEX - BC - PG vastly outperform the baselines, while the randomized version of FLEX - BC - PG has similar performance to the others BCD algorithms. This difference in performance is due to the fact that our deterministic algorithm can exploit the knowledge of the order of the updates to avoid inefficient back and forth between the blocks, which requires computing an inverse wavelet transform and a wavelet transform to take into account the contribution of the other blocks.

Furthermore, we also compare our update rule to that of the "greedy" literature. By computing the partial gradient with respect to each block at each iteration of our FLEX - BC - PG algorithm, we note that the norm of the partial gradient associated with the approximation coefficients is higher than the norm of the other partial gradients. This gap decreases along the iterations. This indicates that a Gauss-Southwell
<figure><img src="image_6.png" /><figcaption>**Fig. 4:** Comparison of the convergence of FB (red), cyclic BC - PG (green), random BC - PG (dark green), FLEX - BC - PG (black), Stochastic FLEX - BC - PG (cyan), and Alternating (PI) FLEX - BC - PG (blue) for the deconvolution problem regularized with 2-Level log sum-Haar wavelet on a 1024 × 1024 image of the Cameraman. Degradation: Gaussian noise with σ/noise = 0.01 and a Gaussian blur of size 40 × 40 and 7 standard deviation. Parameters choice: λ_a = 1 × 10⁻¹⁰, λ_d = 1 × 10⁻⁴, m = 8.</figcaption></figure>

rule or a Gauss-Southwell-Lipschitz<fnref n="2" /> rule [7] would behave as our proposed algorithm. This further validates the interest of our framework, since Gauss-Southwell BCD only guarantee the decrease of objective function value, and in the convex case.

**6 Conclusion**

In this paper we introduce a general block-coordinate forward-backward algorithm, whose convergence is guaranteed in a non-convex setting for a wide range of update rules, encompassing known ones, e.g., cyclic, essentially cyclic, but also new ones, e.g., those inspired by multilevel algorithms. We show on a high dimensional problem that several instances of our FLEX – BC – PG algorithm are competitive with respect to BCD algorithms from the literature and the standard FB algorithm. This general algorithm widens the applicability of BCD approaches, whose update rules can now fully exploit the structure of the problem at hand.

**Acknowledgements.** This work was partially supported by the Fondation Simone et Cino Del Duca - Institut de France and the MEPHISTO (ANR-24-CE23-7039-01) project of the French National Agency for Research (ANR). The work of Luis Briceño-Arias was supported by the National Agence of Research and Development (ANID) from Chile, under the grants FONDECYT 1230257, MATH-AmSud 23-MATH-17, and Centro de Modelamiento Matemático (CMM) BASAL fund FB210005 for centers of excellence.

<footnote n="2">²The Gauss-Southwell-Lipschitz rule divides the norm of each partial gradient by its associated Lipschitz constant. Here the Lipschitz constant is equal to 1 for all possible configuration of blocks.</footnote>
References

[1] Nesterov, Y.: Efficiency of coordinate descent methods on huge-scale optimization problems. SIAM Journal on Optimization **22**(2), 341–362 (2012)

[2] Wright, S.J.: Coordinate descent algorithms. Mathematical programming **151**(1), 3–34 (2015)

[3] Larsson, J., Klopfenstein, Q., Massias, M., Wallin, J.: Coordinate descent for slope. Proceedings of Machine Learning Research, vol. 206, pp. 4802–4821 (2023)

[4] Friedman, J., Hastie, T., Tibshirani, R.: Sparse inverse covariance estimation with the graphical lasso. Biostatistics **9**(3), 432–441 (2008)

[5] Xu, Y., Yin, W.: A block coordinate descent method for regularized multi-convex optimization with applications to nonnegative tensor factorization and completion. SIAM Journal on Imaging Sciences **6**(3), 1758–1789 (2013)

[6] Xu, Y., Yin, W.: A globally convergent algorithm for nonconvex optimization based on block coordinate update. J Sci Comput **72**, 700–734 (2017)

[7] Nutini, J., Laradji, I., Schmidt, M.: Let's make block coordinate descent converge faster: faster greedy rules, message-passing, active-set complexity, and superlinear convergence. Journal of Machine Learning Research **23**(131), 1–74 (2022)

[8] Nutini, J., Schmidt, M., Laradji, I., Friedlander, M., Koepke, H.: Coordinate descent converges faster with the Gauss-Southwell rule than random selection. In: Proc. ICML’15 (2015)

[9] Salzo, S., Villa, S.: Parallel random block-coordinate forward-backward algorithm: a unified convergence analysis. Mathematical Programming **193**(1), 225–269 (2022)

[10] Zeng, J., Lau, T.T.-K., Lin, S., Yao, Y.: Global convergence of block coordinate descent in deep learning. In: Proceedings of the 36th International Conference on Machine Learning. Proceedings of Machine Learning Research, vol. 97, pp. 7313–7323 (2019)

[11] Gratton, S., Mercier, V., Riccietti, E., Toint, P.L.: A block-coordinate approach of multi-level optimization with an application to physics-informed neural networks. Computational Optimization and Applications **89**(2), 385–417 (2024)

[12] Zhang, Z., Brand, M.: Convergent block coordinate descent for training tikhonov regularized deep neural networks. In: Advances in Neural Information Processing Systems (NIPS) (2017)

[13] Hong, M., Wang, X., Razaviyayn, M., Luo, Z.-Q.: Iteration complexity analysis of block coordinate descent methods. Math. Program. **163**(1-2), 85–114 (2017)
[14] Luo, Z.-Q., Tseng, P.: Error bounds and convergence analysis of feasible descent methods: A general approach. Annals of Operations Research **46-47**(1), 157–178 (1993) https://doi.org/10.1007/BF02096261

[15] Powell, M.J.D.: On search directions for minimization algorithms. Mathematical Programming **4**(1), 193–201 (1973) https://doi.org/10.1007/BF01585765

[16] Briceño-Arias, L., Deride, J., Vega, C.: Random activations in primal-dual splittings for monotone inclusions with a priori information. J Optim Theory Appl **192**, 56–81 (2022)

[17] Combettes, P.L., Pesquet, J.-C.: Stochastic Quasi-Fejér block-coordinate fixed point iterations with random sweeping. SIAM Journal on Optimization **25**(2), 1221–1248 (2015)

[18] Lin, Q., Lu, Z., Xiao, L.: An accelerated randomized proximal coordinate gradient method and its application to regularized empirical risk minimization. SIAM Journal on Optimization **25**(4), 2244–2273 (2015)

[19] Richtárik, P., Takáč, M.: Iteration complexity of randomized block-coordinate descent methods for minimizing a composite function. Mathematical Programming **144**, 1–38 (2014)

[20] Richtárik, P., Takáč, M.: Parallel coordinate descent methods for big data optimization. Mathematical Programming **156**, 433–484 (2016)

[21] Fercoq, O., Richtárik, P.: Accelerated, parallel, and proximal coordinate descent. SIAM Journal on Optimization **25**(4), 1997–2023 (2015)

[22] Cadoni, S., Chouzenoux, E., Pesquet, J.-C., Chaux, C.: A block parallel majorize-minimize memory gradient algorithm. In: 2016 IEEE International Conference on Image Processing (ICIP), pp. 3194–3198 (2016). IEEE

[23] Namkoong, H., Sinha, A., Yadlowsky, S., Duchi, J.C.: Adaptive sampling probabilities for non-smooth optimization. In: International Conference on Machine Learning, pp. 2574–2583 (2017). PMLR

[24] Lee, C.-P., Wright, S.J.: Random permutations fix a worst case for cyclic coordinate descent. IMA Journal of Numerical Analysis **39**(3), 1246–1275 (2019)

[25] Sun, R., Ye, Y.: Worst-case complexity of cyclic coordinate descent: $o(n^2)$ gap with randomized version. Mathematical Programming **185**, 487–520 (2021)

[26] Patrascu, A., Necoara, I.: Efficient random coordinate descent algorithms for large-scale structured nonconvex optimization. Journal of Global Optimization **61**(1), 19–46 (2015) https://doi.org/10.1007/s10898-014-0151-9
[27] Ermol'ev, Y.M.: On convergence of random quasi-Féjer sequences. Cybernetics 7, 655-656 (1971)

[28] Alacaoglu, A., Fercoq, O., Cevher, V.: On the convergence of stochastic primal-dual hybrid gradient. SIAM Journal on Optimization 32(3) (2022)

[29] Chambolle, A., Delplancke, C., Ehrhardt, M.J., Schönlieb, C.-B., Tang, J.: Stochastic primal-dual hybrid gradient algorithm with adaptive step-sizes. Journal of Mathematical Imaging and Vision 66(3), 294-313 (2024) https://doi.org/10.1007/s10851-024-01174-1. Early online: 16 March 2024

[30] Bolte, J., Sabach, S., Teboulle, M.: Proximal alternating linearized minimization for nonconvex and nonsmooth problems. Mathematical Programming 146(1-2), 459-494 (2014)

[31] Combettes, P.L.: Hilbertian convex feasibility problem: Convergence of projection methods. Applied Mathematics and Optimization 35(3), 311-330 (1997) https://doi.org/10.1007/BF02683333

[32] Luo, Z.-Q., Tseng, P.: On the linear convergence of descent methods for convex essentially smooth minimization. SIAM Journal on Control and Optimization 30(2), 408-425 (1992) https://doi.org/10.1137/0330022

[33] Mokhtari, A., Gürbüzbalaban, M., Ribeiro, A.: Surpassing gradient descent provably: A cyclic incremental method with linear convergence rate. SIAM Journal on Optimization 28(2), 1420-1447 (2018) https://doi.org/10.1137/16M1101702

[34] Latafat, P., Themelis, A., Patrinos, P.: Block-coordinate and incremental aggregated proximal gradient methods for nonsmooth nonconvex problems. Mathematical Programming 193, 195-224 (2022) https://doi.org/10.1007/s10107-020-01599-7

[35] Fang, H., Fan, Z., Sun, Y., Friedlander, M.: Greed meets sparsity: Understanding and improving greedy coordinate descent for sparse optimization. In: Chiappa, S., Calandra, R. (eds.) Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics. Proceedings of Machine Learning Research, vol. 108, pp. 434-444. PMLR, ??? (2020). https://proceedings.mlr.press/v108/fang20a.html

[36] Chouzenoux, E., Pesquet, J.-C., Repetti, A.: A block coordinate variable metric forward-backward algorithm. Journal of Global Optimization 66(3), 457-485 (2016)

[37] Attouch, H., Bolte, J.: On the convergence of the proximal algorithm for nonsmooth functions involving analytic features. Mathematical Programming 116(1-2), 5-16 (2009)
[38] Attouch, H., Bolte, J., Redont, P., Soubeyran, A.: Proximal alternating minimization and projection methods for nonconvex problems: An approach based on the Kurdyka-Lojasiewicz inequality. Mathematics of operations research 35(2), 438-457 (2010)

[39] Bolte, J., Daniilidis, A., Ley, O., Mazet, L.: Characterizations of Łojasiewicz inequalities: subgradient flows, talweg, convexity. Trans. Am. Math. Soc. 362(6), 3319-3363 (2010)

[40] Sun, Y., Liu, J., Kamilov, U.: Block Coordinate Regularization by Denoising. In: Advances in Neural Information Processing Systems, vol. 32. Curran Associates, Inc., ??? (2019)

[41] Pascal, B., Pustelnik, N., Abry, P., Pesquet, J.-C.: Block-Coordinate Proximal Algorithms for Scale-Free Texture Segmentation. In: 2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 1253-1257. IEEE, Calgary, AB (2018)

[42] Onose, A., Carrillo, R.E., Repetti, A., McEwen, J.D., Thiran, J.-P., Pesquet, J.-C., Wiaux, Y.: Scalable splitting algorithms for big-data interferometric imaging in the SKA era. Monthly Notices of the Royal Astronomical Society 462(4), 4314-4335 (2016)

[43] Parpas, P.: A Multilevel Proximal Gradient Algorithm for a Class of Composite Optimization Problems. SIAM Journal on Scientific Computing 39(5), 681-701 (2017)

[44] Javaherian, A., Holman, S.: A Multi-Grid Iterative Method for Photoacoustic Tomography. IEEE Transactions on Medical Imaging (3), 696-706 (2017)

[45] Hovhannisyan, V., Panagakis, Y., Parpas, P., Zafeiriou, S.: Fast multilevel algorithms for compressive principal component pursuit. SIAM Journal on Imaging Sciences 12(1), 624-649 (2019)

[46] Fung, S.W., Wendy, Z.: Multigrid Optimization for Large-Scale Ptychographic Phase Retrieval. SIAM Journal on Imaging Sciences 13(1), 214-233 (2020)

[47] Buccini, A., Donatelli, M.: A multigrid frame based method for image deblurring. Electronic Transactions on Numerical Analysis 53, 283-312 (2020)

[48] Plier, J., Savarino, F., Kočvara, M., Petra, S.: First-Order Geometric Multilevel Optimization for Discrete Tomography. In: Scale Space and Variational Methods in Computer Vision vol. 12679, pp. 191-203. Springer, Cham (2021). Series Title: Lecture Notes in Computer Science

[49] Lauga, G., Riccietti, E., Pustelnik, N., Gonçalves, P.: Multilevel Fista For Image Restoration. IEEE ICASSP, Rhodes, Greece (4-10 June 2023)
[50] Lauga, G., Riccietti, E., Pustelnik, N., Gonçalves, P.: IML FISTA: A Multilevel Framework for Inexact and Inertial Forward-Backward. Application to Image Restoration. SIAM Journal on Imaging Sciences 17(3), 1347–1376 (2024)

[51] Lauga, G., Repetti, A., Riccietti, E., Pustelnik, N., Gonçalves, P., Wiaux, Y.: A multilevel framework for accelerating uSARA in radio-interferometric imaging. In: 2024 32nd European Signal Processing Conference (EUSIPCO), pp. 2287–2291 (2024)

[52] Laurent, N., Tachella, J., Riccietti, E., Pustelnik, N.: Multilevel plug-and-play image restoration. Preprint hal-05004914 (2025)

[53] Calandra, H., Gratton, S., Riccietti, E., Vasseur, X.: On high-order multilevel optimization strategies. SIAM Journal on Optimization 31(1), 307–330 (2021)

[54] Mallat, S.: A Wavelet Tour of Signal Processing. Elsevier, ??? (1999)

[55] Ortega, J.M., Rheinboldt, W.C.: Iterative Solution of Nonlinear Equations in Several Variables. Classics in Applied Mathematics, vol. 30. SIAM, Philadelphia (2000)

[56] Tseng, P.: Convergence of a block coordinate descent method for nondifferentiable minimization. Journal of optimization theory and applications 109, 475–494 (2001)

[57] Abboud, F., Chouzenoux, É., Pesquet, J.-C., Chenot, J.-H., Laborelli, L.: An alternating proximal approach for blind video deconvolution. Signal Processing: Image Communication 70, 21–36 (2019)

[58] Rockafellar, R.T., Wets, R.J.-B.: Variational Analysis vol. 317. Springer, ??? (2009)

[59] Repetti, A., Wiaux, Y.: Variable metric forward-backward algorithm for composite minimization problems. SIAM Journal on Optimization 31(2), 1215–1241 (2021)

[60] Beck, A., Teboulle, M.: Smoothing and First Order Methods: A Unified Framework. SIAM Journal on Optimization 22(2), 557–580 (2012) https://doi.org/10.1137/100818327

[61] Lauga, G.: Multilevel proximal methods and application to image restoration. phdthesis, École Normale Supérieure de Lyon (December 2024)

[62] Prater-Bennette, A., Shen, L., Tripp, E.E.: The proximity operator of the log-sum penalty. Journal of Scientific Computing 93(3), 67 (2022)

[63] Bertsekas, D.P.: Nonlinear Programming, 2nd edn. Athena Scientific, Belmont, MA (1999)
[64] Hansen, P.C., Nagy, J.G., O'Leary, D.P.: Deblurring Images. SIAM, ??? (2006)

Appendix A Supplementary material for
convergence proofs.

In this section, we report the descent lemmas we used to assert the convergence of our
FLEX - BC - PG algorithm.

**Lemma A.1.** *Descent lemma [55, 63].* Let $f: \mathbb{R}^N \to \mathbb{R}$ be a continuously differentiable function with Lipschitz continuous gradient and Lipschitz constant $\beta_f$. Then for any $\beta \ge \beta_f$,

$$
f(x) \le f(y) + \langle x - y, \nabla f(y) \rangle + \frac{\beta_f}{2} \|x - y\|^2 \quad \text{for every } x, y \in \mathbb{R}^N. \quad (A1)
$$

**Lemma A.2. (Non-convex) proximal-gradient descent lemma [30].** Let $f: \mathbb{R}^N \to \mathbb{R}$ be a continuously differentiable function with Lipschitz continuous gradient and Lipschitz constant $\beta_f$. Let $g: \mathbb{R}^N \to \mathbb{R}$ be a proper, lower semicontinuous function with $\inf_{\mathbb{R}^N} g > -\infty$. If

$$
y \in \operatorname{prox}_{\tau g}(x - \tau \nabla f(x)), \quad (A2)
$$

then for any $0 < \tau < \frac{1}{\beta_f}$

$$
f(y) + g(y) + \frac{1}{2} \left( \frac{1}{\tau} - \beta_f \right) \|x - y\|^2 \le f(x) + g(x). \quad (A3)
$$

*Proof.* First $\operatorname{prox}_{\tau g}(\cdot)$ is well-defined by [30, Proposition 2]. Thus, for all $\mathbf{x} \in \mathbb{R}^N$, there exists $y \in \operatorname{prox}_{\tau g}(x - \tau \nabla f(x))$. This inequality comes directly from [30, Lemma 2], but for completeness of the argument we reproduce it here. By definition of the proximity operator:

$$
y \in \arg \min_{z \in \mathbb{R}^N} \langle z - x, \nabla f(x) \rangle + g(z) + \frac{1}{2\tau} \|z - x\|^2
$$

Thus taking $z = x$ we obtain

$$
\begin{align*}
\langle y - x, \nabla f(x) \rangle + g(y) + \frac{1}{2\tau} \|y - x\|^2 &\le \langle x - x, \nabla f(x) \rangle + g(x) + \frac{1}{2\tau} \|x - x\|^2 \\
&\le g(x)
\end{align*}
$$

Now invoking Lemma A.1, we have:

$$
\langle \nabla f(x), x - y \rangle \le f(x) - f(y) + \frac{\beta_f}{2} \|x - y\|^2
$$
which yields for any $0 < \tau < \frac{1}{\beta_f}$

$$f(y) + g(y) + \frac{1}{2} \left( \frac{1}{\tau} - \beta_f \right) \|x - y\|^2 \le f(x) + g(x).$$

□

# Appendix B Convergence of a stochastic FLEX - BC - PG for convex optimization

In this section, we briefly present a convergence result for a randomized version of our Flexible Block-Coordinate Forward-Backward algorithm. The convergence result in itself is a direct application of [9, Theorem 4.9]. We aim here to construct a stochastic BC FB that, in expectation, mirrors the behavior of our multilevel algorithm and is convergent. Such algorithm follows classic rules of stochastic BCD algorithms that can update the blocks in parallel.

With such algorithm we will be able to have a complete comparison of the update rules available for BC descent algorithms. Recall that the algorithm is of the following form: Let $(\varepsilon^n)_{n \in \mathbb{N}} = (\varepsilon_1^n, \dots, \varepsilon_L^n)_{n \in \mathbb{N}}$ be a sequence of variables with value in $\{0, 1\}^L$. Let $(\tau_\ell)_{1 \le \ell \le L} \in \mathbb{R}_{++}^L$ and $\mathbf{x}^0 = (x_1^0, \dots, x_L^0) \in \text{dom } g$. Iterate

$$ \begin{array}{l} \text{for } n = 0, 1, \dots \\ \quad \left[ \begin{array}{l} \text{for } \ell = 1, \dots, L \\ \quad \left[ x_{\ell}^{n+1} = x_{\ell}^n + \varepsilon_{\ell}^n (\text{prox}_{\tau_{\ell} g_{\ell}} (x_{\ell}^n - \tau_{\ell} \nabla_{\ell} f(\mathbf{u}^n)) - x_{\ell}^n) \end{array} \end{array} \right. \end{array} \tag{B4} $$

The main difference w.r.t. the paradigm proposed in the previous sections is that $\varepsilon$ and $\mathbf{x}$ are now random variables. Consider the following assumptions:

## Assumption 5.

A6 $f: \mathcal{H} \to \mathbb{R}$ is convex and continuously differentiable,

A7 for every $\ell = 1, \dots, L$, $g_\ell: \mathcal{H}_\ell \to (-\infty, +\infty]$ is proper, convex, and lower semicontinuous.

A8 $\varepsilon = (\varepsilon_1, \dots, \varepsilon_L)$ is a random variable with values in $\{0, 1\}^L$, such that for every $\ell \in \{1, \dots, L\}$, $P(\varepsilon_\ell = 1) > 0$ and $P(\varepsilon = (0, \dots, 0)) = 0$.

We can now present a way to construct update rules to mimic our multilevel algorithm that verify Assumption A8. As multilevel algorithms mostly employ V-scheme in practice [43, 44, 46, 49–51], we present an update rule for this scheme.

**Lemma B.1.** *V-scheme probabilities for stochastic FLEX-BC-PG. Suppose that $\varepsilon = (\varepsilon_1, \dots, \varepsilon_L)$ is a random variable with values in $\{0, 1\}^L$, such that $P(\varepsilon_1 = 1) = 1$ and for every $\ell \in \{1, \dots, L-1\}$,*

* $P(\varepsilon_{\ell+1} = 1 | \varepsilon_\ell = 1) > 0$,

* $P(\varepsilon_{\ell+1} = 1 | \varepsilon_\ell = 0) = 0$.

Then, for every $\ell \in \{1, \dots, L\}$, $P(\varepsilon_\ell = 1) > 0$ and $P(\varepsilon = (0, \dots, 0)) = 0$.
Proof. The second point is straightforward. For the first point, simply remark that for every $\ell \in \{2, \dots, L\}$:

$$P(\varepsilon_{\ell} = 1) = P(\varepsilon_{\ell} = 1 | \varepsilon_{\ell-1} = 1)P(\varepsilon_{\ell-1} = 1),$$

then one directly has:

$$P(\varepsilon_{\ell} = 1) = \left( \prod_{j=2}^{j=\ell} P(\varepsilon_{j} = 1 | \varepsilon_{j-1} = 1) \right) P(\varepsilon_{1} = 1).$$

which is strictly greater than 0. □

One can see that with this construction we will update the coarsest level at each iteration, and that updating "fine" levels will also force us to update coarser levels, which is typical of multilevel methods.

The sampling of $\varepsilon$ is done sequentially by increasing $\ell$ until we reach the first zero occurrence. In order to update all levels as often as possible, the value of $P(\varepsilon_{\ell+1} = 1 | \varepsilon_{\ell} = 1)$ should be close to 1 for large $\ell$.

### Choosing the right value for the conditional probabilities.

In a typical V-scheme, a multilevel algorithm would compute $m$ iterations at each coarse level, going upwards in the resolution. After that it would compute one iteration at fine level. Thus, we should adjust the conditional probabilities of activating each block so that with high probability we update $m \ge 0$ times the coarsest level alone, then $m$ times the coarsest level and the second to last coarsest level, and so on. We thus impose for all $\ell$:

$$P(\varepsilon_{\ell} = 1) = \left(\frac{1}{m}\right)^{\ell},$$

which yields:

$$P(\varepsilon_{\ell+1} = 1 | \varepsilon_{\ell} = 1) = \frac{1}{m}.$$

### Convergence of the stochastic algorithm.

We can now state the convergence result for the stochastic version of our algorithm. The proof is a direct application of [9, Theorem 4.9] and is therefore omitted. We denote by $E$ the expected value.

**Theorem B.2.** *Convergence of stochastic FLEX-BC-PG [9, Theorem 4.9].* Let $(\varepsilon_n)_{n \in \mathbb{N}} = (\varepsilon_1^n, \dots, \varepsilon_L^n)_{n \in \mathbb{N}}$ be a sequence of independent copies of $\varepsilon$. Let $(\tau_\ell)_{1 \le \ell \le L} \in \mathbb{R}_{++}^L$ and $x_0 = (x_{1,0}, \dots, x_{L,0}) = (x^0, \dots, x^0) \in \text{dom } g$ be a constant random variable. Set $\delta = \max_{1 \le \ell \le L} \tau_\ell \beta_\ell$ (the block Lipschitz constants, see Assumption 2) and $p_{\min} = \min_{1 \le \ell \le L} P(\varepsilon_\ell = 1)$.

Set $\mathbf{Id} = \bigoplus_{\ell=1}^L \frac{1}{\tau_\ell P(\varepsilon_\ell = 1)} \mathbf{Id}_\ell$ (the identity operators on $\mathcal{H}_\ell$), $\Psi_* = \inf \Psi$, and $S_* = \arg \min \Psi \subset \mathcal{H}$. Then the following hold.

$$(i) \quad E[\Psi(\mathbf{x}^n)] \to \Psi_*.$$
(ii) Suppose that $S_* \neq \emptyset$. Then $E[\Psi(\mathbf{x}^n)] - \Psi_* = o(1/n)$ and for every integer $n \ge 1$,

$$E[\Psi(\mathbf{x}^n)] - \Psi_* \le \left[ \frac{\mathrm{dist}_{\mathrm{Id}}^2(x_0, S_*)}{2} + \left( \frac{\max\{1, (2-\delta)^{-1}\}}{\mathrm{p}_{\mathrm{min}}} - 1 \right) (\Psi(\mathbf{x}^0) - \Psi_*) \right] \frac{1}{n}$$

Moreover there exists a random variable $x_*$ taking values in $S_*$ such that $\mathbf{x}^n \to x_*$.

# Appendix C Implementation details for the numerical experiments

## C.1 Efficient computation of the gradient of the approximation

In our numerical experiments, the degradation is a Gaussian blur. As this blur is symmetric, the blurring matrix $A \in \mathbb{R}^{N \times N}$ can be expressed as a Kronecker product [64]

$$A = A_r \otimes A_c$$

where $A_r$ and $A_c$ are $\sqrt{N} \times \sqrt{N}$ real matrices that decompose the action of the blur into its vertical ($c$ for columns) and horizontal ($r$ for rows) components. Using the following relationship

$$Ax = \mathrm{vec}(A_cXA_r^\top)$$

where $\mathrm{vec}$ denotes the vectorization, and $X \in \mathbb{R}^{\sqrt{N} \times \sqrt{N}}$ is our image in its matrix form, we can avoid storing $A$ and exploit a similar relationship to compute the gradient of the "coarse" function w.r.t. the approximation coefficients without having to compute the global gradient.

Indeed, under similar conditions as for the Gaussian blur [50], the projection operation $\Pi_V$ can be written using a Kronecker product. We have [50, Section 3.2] that:

$$\Pi_V = R_{q,r} \otimes R_{q,c}$$

where $R_{q,r}$, $R_{q,c} \in \mathbb{R}^{\sqrt{N} \times \sqrt{N}/2}$ are Toeplitz matrices generated from the quadrature mirror filter $q$ [50, Section 3.2] of the wavelet transform that defined $\Pi_V$. For square images, $R_{q,r} = R_{q,c} := R$.

Therefore the gradient of the data fidelity term w.r.t. the approximation coefficients can be expressed as:

$$\Pi_V A^* (A \Pi_V^* a^n - \Pi_V^* \Pi_V z) = R A_c^\top A_c R^\top a_o R A_r^\top A_r R^\top - R A_c^\top R^\top R_z R^\top R A_r R^\top \quad (C5)$$

As $R A_c^\top A_c R^\top$, $R A_r^\top A_r R^\top$, $R A_c^\top R^\top$, $R_z R^\top$, $R A_r R^\top$, all belong to $\mathbb{R}^{N_H, N_H}$ and can be pre-computed, the gradient of $1/2 \|A \cdot -z\|_2^2$ w.r.t. the approximation coefficients can be evaluated efficiently.