Mila Nikolova
email: nikolova@cmla.ens-cachan.fr

Pauline Tan
email: pauline.tan@cmla.ens-cachan.fr

Alternating structure-adapted proximal gradient descent for nonconvex block-regularised problems

Keywords: Alternating minimization, Block coordinate descent, Global convergence, Kurdyka-Lojasiewicz property, Nonconvex-nonsmooth optimization, Proximal forward-backward, Proximal gradient descent, Subanalytic functions

There has been an increasing interest in constrained nonconvex regularized block biconvex / multiconvex optimization problems. We introduce an approach that effectively exploits the biconvex / multiconvex structure of the coupling term and enables complex application-dependent regularization terms to be used. The proposed Alternating Structure-Adapted Proximal gradient descent algorithm enjoys simple well defined updates. Global convergence of the algorithm to a critical point is proved using the so-called Kurdyka-Lojasiewicz property for subanalytic functions. Moreover, we prove that a large class of useful objective functions obeying our assumptions are subanalytic and thus satisfy the Kurdyka-Lojasiewicz property.

Introduction

Recently, there has been an increasing interest in the design and the analysis of regularized block biconvex (multiconvex) optimization problems. In the first part of this work we consider problems of the form: minimize x,y J(x, y) = F (x) + G(y) + H(x, y),

where x and y belong to real finite-dimensional real spaces. Such an objective is also known as a two block optimization model with blocks x and y. The functions F and G can be nonconvex and are supposed continuously differentiable inside the domain of J (for example, smooth approximations of nonsmooth functions). It worths noting that J is generally nonconvex even if F and G are convex.

Then we consider the extension of (1) to N blocks living in finite-dimensional real spaces {U i } N i=1 , in which case J :

U 1 × • • • × U N → R ∪ {+∞} reads as minimize z J z (1) , . . . , z (N ) := N i=1 F i z (i) + H(z).
(

) 2
The coupling function H is smooth enough, and the possibly nonconvex functions F i : U i → R ∪ {+∞} are continuously differentiable on the domain of J. In many practical cases, H is multiconvex, i.e., z (i) → H(z) is convex for any i.

Optimization problems of the form (1) -( 2) are widely used in engineering, science and finance. They are rich enough to cover various practical applications, such as blind source separation, blind deconvolution, nonnegative matrix factorization, structured total least squares, multi-modal learning for image classification [START_REF] Zhang | Classify social image by integrating multi-modal content[END_REF], patch-based methods for inpainting [START_REF] Arias | Analysis of a variational framework for exemplar-based image inpainting[END_REF], to list a few.

In order to simplify the reading of the paper, most of the presentation is on the two block problem [START_REF] Absil | Convergence of the iterates of descent methods for analytic cost functions[END_REF]; the results on the multi block problem [START_REF] Aguerrebere | A hyperprior bayesian approach for solving image inverse problems[END_REF] are given later on.

General alternating minimization schemes

The most intuitive way to solve problems of the form given in (1) is to use alternating minimization which which generates a sequence of iterates {x k , y k } k∈N defined by Choose algorithm ALG to solve the x-and y-updates; for each k = 1, 2, . . . , compute

x k ∈ arg min x J(x, y k-1 ) using ALG, y k ∈ arg min y J(x k , y) using ALG.

(3)

This classic scheme is known as block coordinate Gauss-Seidel method and as block coordinate descent (BCD). It was introduced in [START_REF] Hildreth | A quadratic programming procedure[END_REF] (1957) and further developed for various problems, see e.g., [START_REF] Tseng | Convergence of a block coordinate descent method for nondifferentiable minimization[END_REF][START_REF] Gorski | Biconvex sets and optimization with biconvex functions: a survey and extensions[END_REF][START_REF] Arias | Analysis of a variational framework for exemplar-based image inpainting[END_REF][START_REF] Beck | On the convergence of block coordinate descent type methods[END_REF][START_REF] Aguerrebere | A hyperprior bayesian approach for solving image inverse problems[END_REF]. If J is continuously differentiable and if the minimum in each step is uniquely attained, convergence to a critical point holds [START_REF] Bertsekas | Nonlinear programming[END_REF]Prop. 2.7.1]. A general convergence result on descent methods for real-analytic (possibly nonconvex) objectives was obtained by Absil, Mahony, and Andrews in [START_REF] Absil | Convergence of the iterates of descent methods for analytic cost functions[END_REF].

A way to relax the requirements for convergence of the BCD in ( 3) is to consider the proximal BCD scheme:

Choose algorithm ALG to solve the x-and y-updates; Select step-sizes τ > 0 and σ > 0;

for each k = 1, 2, . . . , compute

x k ∈ arg min x J(x, y k-1 ) + 1 2τ x -x k-1 2 using ALG, y k ∈ arg min y J(x k , y) + 1 2σ y -y k-1 2 using ALG.

(

) 4
This approach was introduced for convex functions J by Auslender [7, sec. 4]. An extension to J strictly quasi-convex with respect to its blocks was proposed in [25, sec. 7] and to nonsmooth objectives in [START_REF] Tseng | Convergence of a block coordinate descent method for nondifferentiable minimization[END_REF].

Convergence facts on (4) for other nonconvex nonsmooth objectives were found in [START_REF] Xu | A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion[END_REF].

Note that the BCD and the proximal BCD schemes in [START_REF] Arias | Analysis of a variational framework for exemplar-based image inpainting[END_REF] and in [START_REF] Attouch | On the convergence of the proximal algorithm for nonsmooth functions involving analytic features[END_REF], respectively, employ a minimization algorithm ALG which heavily determines the numerical issue.

Review of related literature

An efficient approach to solve the difficulties arising with the proximal BCD for nonconvex and nonsmooth objectives J was proposed in 2013 by Xu and Yin [START_REF] Xu | A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion[END_REF] for block multiconvex differentiable H and by Bolte, Sabach, and Teboulle [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF] for two block continuously differentiable H. Using the smoothness of H, the idea was to apply a proximal linearized BCD for (x k , y k ):

for each k = 1, 2, . . . , compute find τ k > 0 according to Lip ∇ x H •, y k-1

x k ∈ arg min x x, ∇ x H x k-1 , y k-1 + F (x) + 1 2τ k x -x k-1 2 , find σ k > 0 according to Lip ∇ y H x k , • y k ∈ arg min y y, ∇ y H x k , y k-1 + G(y) + 1 2σ k y -y k-1 2 , (5)
where the Lipschitz constant is computed (or estimated) at each step. The scheme needs regularizers (F, G) that are "simple" in the sense that their proximity operator

arg min x F (x) + 1 2τ x -z 2 (6)
has a closed-form expression. Since then, this approach became very popular. It was further successfully used and improved in numerous works, see e.g., [START_REF] Hesse | Proximal heterogeneous block implicit-explicit method and application to blind ptychographic diffraction imaging[END_REF][START_REF] Xu | Alternating proximal gradient method for sparse nonnegative Tucker decomposition[END_REF][START_REF] Beck | An alternating semiproximal method for nonconvex regularized structured total least squares problems[END_REF][START_REF] Pock | Inertial proximal alternating linearized minimization (iPALM) for nonconvex and nonsmooth problems[END_REF][START_REF] Chouzenoux | A block coordinate variable metric forwardbackward algorithm[END_REF][START_REF]A globally convergent algorithm for nonconvex optimization based on block coordinate update[END_REF]. We recall that a proximal linearized BCD is equivalent to an alternating proximal gradient descent and to an alternating proximal forward-backward. The advantages of this approach compared to the schemes in subsection 1.1 are tremendous. All the three schemes -the BCD (3), the proximal BCD (4) and the proximal linearized BDC (5) -were analyzed and compared in [START_REF] Xu | A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion[END_REF]. The conclusion [38, p. 1795] is that in general the three schemes give different solutions, and that the proximal linearized BDC (5) needs less computation and offers a better decrease of the objective function than the other algorithms. One however remarks the cost to compute or to estimate the step-sizes (τ k , σ k ) at each iteration and the possible variability of {τ k } k≥0 and {σ k } k≥0 during the iterations. An attempt to fix this issue is the preconditioning proposed by Chouzenoux, Pesquet, and Repetti [START_REF] Chouzenoux | A block coordinate variable metric forwardbackward algorithm[END_REF]. One can also remark that if F or G is truly nonconvex, the solution in (6) can be composed out of non-connected sets, for an example see [17, p. 487].

Remark 1. [Choices for H in applications]

In nearly all applications solved using a scheme of the form (5), the coupling term H is biquadratic (multiquadratic), sometimes combined with a bilinear term. When restricted to two blocks, one finds H(x, y) = L 0 (x • y) -w 2 + L 1 (x), L 2 (x, y) where "•" denotes a product, w is a known matrix, • stands for Frobenius norm and L i are linear forms. Such a term is used in [START_REF] Beck | An alternating semiproximal method for nonconvex regularized structured total least squares problems[END_REF] and with L 0 = 1 and L 1 = L 2 = 0 in [START_REF] Xu | A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion[END_REF][START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF][START_REF] Pock | Inertial proximal alternating linearized minimization (iPALM) for nonconvex and nonsmooth problems[END_REF][START_REF]A globally convergent algorithm for nonconvex optimization based on block coordinate update[END_REF]. In all these applications, (∇ x H, ∇ y H) are only locally Lipschitz.

Remark 2. [Choices for (F, G) in applications] The most typical form of F (resp. G) is the non-negativity constraint x ≥ 0 (resp., y ≥ 0) -see [START_REF] Xu | A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion[END_REF][START_REF] Hesse | Proximal heterogeneous block implicit-explicit method and application to blind ptychographic diffraction imaging[END_REF][START_REF] Pock | Inertial proximal alternating linearized minimization (iPALM) for nonconvex and nonsmooth problems[END_REF][START_REF]A globally convergent algorithm for nonconvex optimization based on block coordinate update[END_REF] and [17, sec. 4.2, p. 486]. We note that in those cases F and G are smooth inside their domains. Other forms for F or G are the 2 norm (possibly squared) [START_REF] Pock | Inertial proximal alternating linearized minimization (iPALM) for nonconvex and nonsmooth problems[END_REF][START_REF] Beck | An alternating semiproximal method for nonconvex regularized structured total least squares problems[END_REF], the counting function 0 [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF] and p , 0 ≤ p < 1 [START_REF]A globally convergent algorithm for nonconvex optimization based on block coordinate update[END_REF].

A unified approach for proving the convergence of proximal splitting methods for nonconvex and nonsmooth problems was developed by Attouch, Bolte, and Svaiter in their seminal work [START_REF] Attouch | Convergence of descent methods for semi-algebraic and tame problems: proximal algorithms, forward-backward splitting, and regularized Gauss-Seidel methods[END_REF]. A central assumption in order to prove global convergence of the iterates to a critical point is that the objective function J satisfies the so-called Kurdyka-Lojasiewicz (K L) property [START_REF] Bolte | The lojasiewicz inequality for nonsmooth subanalytic functions with applications to subgradient dynamical systems[END_REF][START_REF] Bolte | Clarke subgradients of stratifiable functions[END_REF]. In several articles [START_REF] Hesse | Proximal heterogeneous block implicit-explicit method and application to blind ptychographic diffraction imaging[END_REF][START_REF] Beck | An alternating semiproximal method for nonconvex regularized structured total least squares problems[END_REF][START_REF] Pock | Inertial proximal alternating linearized minimization (iPALM) for nonconvex and nonsmooth problems[END_REF][START_REF] Beck | An alternating semiproximal method for nonconvex regularized structured total least squares problems[END_REF] convergence is proven using the methodology proposed in [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF].

Motivation and proposed algorithm ASAP

The coupling term H typically has a good structure in the sense of being "simple" with respect to each variable x and y. In particular, the biconvexity (block multiconvexity) of H is a very strong structural property. The proximal linearized gradient approach in (5) (subsection 1.2) does not benefit from these features. Our motivation is to build an algorithm that exploits the biconvexity (block multiconvexity) of H. This means reversing the splitting used in (5), i.e., using the proximity operators with respect to H(•, y) and H(x, •) -"simple" in practice, see Remark 1 -instead of those of F and G. The argument behind the splitting in ( 5) is to use F and G nonsmooth. Alternatively, we attach to H the convex constraints on F and G and assume that F and G are differentiable on their domains.

Remark 3. Smooth approximations to solve nonsmooth optimization problems are recommended in a unified framework by Beck and Teboulle [START_REF] Beck | Smoothing and first order methods: a unified framework[END_REF]. Very numerous methods with nonsmooth regularizers use their smooth approximations in the numerical algorithms [START_REF] Chen | Smoothing methods for nonsmooth, nonconvex minimization[END_REF][START_REF] Hintermüller | Nonconvex TV q -models in image restoration: Analysis and a trustregion regularization-based superlinearly convergent solver[END_REF][START_REF]Penalty methods for a class of non-Lipschitz optimization problems[END_REF][START_REF] Pock | Inertial proximal alternating linearized minimization (iPALM) for nonconvex and nonsmooth problems[END_REF]. Smooth (stiff) functions are customarily used for sparse recovery [START_REF] Chen | Smoothing nonlinear conjugate gradient method for image restoration using nonsmooth nonconvex minimization[END_REF][START_REF] Lai | An unconstrained -q minimization with 0 < q ≤ 1 for sparse solution of underdetermined linear systems[END_REF][START_REF]Penalty methods for a class of non-Lipschitz optimization problems[END_REF]. A study by Chen of smoothing methods giving rise to efficient calculations in numerical schemes for nonsmooth optimization can be found in [START_REF] Chen | Smoothing methods for nonsmooth, nonconvex minimization[END_REF]. In other cases the estimates of x, y must not be sparse and hence the regularizers must be smooth; a practical application on fringe separation is presented in subsection 7.2.

For simplicity, here we present the case of biconvex coupling terms H. We reformulate the problem so that H contains the biconvex constraints and let F and G be continuously differentiable. Thus we propose the simple Alternating Structure-Adapted Proximal gradient descent (ASAP) algorithm1 sketched below:

ASAP (2 block) Q = 2 if H biconvex, Q = 1 otherwise; set τ ∈ (0, Q/Lip(∇F )) and σ ∈ (0, Q/Lip(∇G));
for each k = 1, 2, . . . , compute:

x k = arg min x x, ∇F x k-1 + H(x, y k-1 ) + 1 2τ x -x k-1 2 ; y k = arg min y y, ∇G y k-1 + H(x k , y) + 1 2σ y -y k-1 2 . ( 7
)
Let us provide some details about the proposed algorithm.

-The step-sizes τ and σ in the proposed ASAP algorithm depend only on F and Gand eventually on the form of the coupling term via Q ∈ {1, 2}, so they are fixed in advance. On can always use Q = 1 but in case when H is biconvex Q = 2 and thus twice larger step-sizes are allowed.

-The iterates x k and y k are always well defined. When H is biconvex, they are uniquely defined (even if F or G are nonconvex). Furthermore, Remark 1 shows that in many important applications (x k , y k ) are the explicit solutions of well-posed quadratic problems. In general, data-fidelity terms are defined using functions H(x, •) and H(•, y) that are "simple" in the sense of [START_REF] Attouch | Convergence of descent methods for semi-algebraic and tame problems: proximal algorithms, forward-backward splitting, and regularized Gauss-Seidel methods[END_REF].

-Regularizers F and G can have complex structure in order to capture various application-dependent features; one example is presented in subsection 7.2.

Our assumptions are very simple and easy to verify. Besides, a series of additional results help understanding the algorithm. We exhibit important classes of optimization problems that can be solved using the proposed ASAP and prove that they are subanalytic. They hence enjoy the K L property. The global convergence of the iterates produced by ASAP is established using a result from [START_REF] Attouch | Convergence of descent methods for semi-algebraic and tame problems: proximal algorithms, forward-backward splitting, and regularized Gauss-Seidel methods[END_REF].

The ASAP algorithm addresses in a simple way a wide range of problems. Two particular applications of the algorithm to big-data air-born sequences of images are actually used by our industrial partner ONERA (the French aerospace lab).

Outline

Section 2 presents the optimization model. Section 3 is devoted to the prerequisites for understanding the ASAP algorithm and its convergence. The ASAP algorithm is outlined in section 4. Section 5 establishes the convergence of ASAP for biconvex coupling terms; a generic family of objective functions is proven to satisfy the K L property in subsection 5.3. Section 6 is on block multiconvex coupling terms. Applications for Hadamard based biquadratic coupling terms are presented in section 7.

The problem and examples

Notations We consider that x ∈ U and y ∈ V where U and V are real finite-dimensional spaces. The ith element of a vector or a matrix x (seen as a vector) reads as x i . A vector (a matrix) indexed for some purpose is denoted by x (i) . For an m × n real matrix w we systematically denote

w := w F = i,j w 2 i,j ,
noticing that if w is a vector (n = 1), the Frobenius norm • F boils down to the 2 norm. For subsets of real finite-dimensional spaces, capital italic letters are used. Given a nonempty set S ⊂ U , the distance of any point x + ∈ U to S is defined by

dist(x + , S) := inf{ x -x + | x ∈ S},
while the indicator function of S is given by

χ S (x) := 0 if x ∈ S, +∞ if x ∈ S.
The subdifferential of a convex function h is denoted by ∂h.

The optimization model for two block problems

We are interested in solving nonconvex minimization problems of the form J :

U × V → R J(x, y) := F (x) + G(y) + H(x, y), (8)
where U and V are real finite-dimensional spaces. 2 According to what was said in the introduction, we adopt the following blanket model assumption on the objective J: (i) H is continuously differentiable on its domain;

Assumption (M) (a) J : U × V → R ∪ {+∞} is lower bounded; (b) F : U → R and G : V → R
(ii) H is biconvex3 and differentiable on its domain.

The cases F = 0 or G = 0 are considered as well. From Assumption (M), the objective J is lowersemicontinuous. Further, H can be cast into the form

H(x, y) = H(x, y) + χ Dx (x) + χ Dy (y),
where H is (continuously) differentiable on U × V and D x ⊆ U and D y ⊆ V are nonempty closed sets. The optimization problem equivalently reads as

J(x, y) = H(x, y) + F (x) + χ Dx (x) + G(y) + χ Dy (y). (9)
In order to unify the presentation in cases (i) and (ii) we introduce

Q = 1 H satisfies (i) 2 H satisfies (ii). ( 10
)

Illustration: a general family of objective functions

A generic family of objective functions that can be minimized using the proposed algorithm is described next for the case of H biconvex:

J(x, y) := i f i ( A i x ) =:F (x) + j g j ( B j y ) =:G(y) + k h k (b k (x, y) -w k ) + χ Dx (x) + χ Dy (y) =:H(x,y)
with b k : U × V → W bilinear forms (D x , D y ) closed nonempty sets, A i and B i linear operators, and w ∈ W given data. Function h is differentiable; the most often, h(t) = t 2 . Some relevant choices for functions (f i , g j ) in accordance with Assumption (M)(b) are listed below.

Example 1. [Choices for (f i , g j ) in (2.2)] These functions are of the form ψ : R → R or (ψ) p : R → R where p is a rational number and include a parameter α > 0:

(i) ψ(t) := |t| 2 and (ψ(t)) p for p > 1 2 ;

(ii)

ψ(t) := |t| -α/2 if |t| > α t 2 /(2 α) if |t| ≤ α and (ψ(t)) p for p ∈ (0, 1]; (iii) ψ(t) := √ t 2 + α and (ψ(t)) p for p ∈ (0, 1]; (iv) ψ(t) := |t| -α log(1 + |t|/α); (v) ψ(t) := log(1 + t 2 /α); (vi) ψ(t) := t 2 /(α + t 2 ); (vii) ψ(t) = 1 -exp -t 2 /α .
When α is small, functions (ii)-(vi) are stiff near the origin and they provide smooth approximations of nonsmooth functions. In particular, (v), (vi) and (vii) can approximate the counting function 0 . Functions ψ in (iii)-(iv) are convex and for α small enough they are used to approximate the 1 norm [START_REF] Chen | Smoothing nonlinear conjugate gradient method for image restoration using nonsmooth nonconvex minimization[END_REF][START_REF] Beck | Smoothing and first order methods: a unified framework[END_REF] whereas ψ p for p ∈ (0, 1] are used to approximate the corresponding p "norm" for "sparse recovery" [START_REF] Lai | An unconstrained -q minimization with 0 < q ≤ 1 for sparse solution of underdetermined linear systems[END_REF][START_REF] Hintermüller | Nonconvex TV q -models in image restoration: Analysis and a trustregion regularization-based superlinearly convergent solver[END_REF][START_REF]Penalty methods for a class of non-Lipschitz optimization problems[END_REF].

3 Preliminary facts

Elements of subdifferential calculus

Here we recall some facts on subdifferential calculus in relation with the objective J. Given a function f : R m → R ∪ {+∞}, its domain is

dom f := {x ∈ R m | f (x) < +∞},
and f is proper if and only if dom f = ∅. Definition 1. [Subgradients of convex functions] Let h : R m → R ∪ {+∞} be a proper convex lower semicontinuous function and x + ∈ dom h. The subdifferential ∂h(x + ) of h at x + is the set of all p ∈ R m , called subgradients of h at x + , such that

∀ x ∈ R m h(x) ≥ h(x + ) + p, x -x + If x + / ∈ dom h, then ∂h(x + ) = ∅.
The subdifferential for nonconvex nonsmooth functions is defined below.

Definition 2. [33, Def. 8.3] Let h : R m → R ∪ {+∞} be a function. (a) The Fréchet subdifferential of h at x + ∈ dom h, denoted ∂h(x + ), is the set of vectors p ∈ R m such that ∀ x ∈ R m one has h(x) ≥ h(x + ) + p, x -x + + o( x -x + ). If x + / ∈ dom h, then ∂h(x + ) = ∅. (b) The (limiting-)subdifferential of h at x + ∈ dom h, written ∂h(x + ), is defined by ∂h(x + ) := p ∈ R m | ∃ x k → x + , h(x k ) → h(x + ), p k → p, p k ∈ ∂h(x k ) .
Remark 4. For any x + ∈ dom h, the sets ∂h(x + ) and ∂h(x + ) are closed, ∂h(x + ) is convex and ∂h(x

+ ) ⊂ ∂h(x + ) [33, Theorem 8.6]. If h is convex, then ∂h(x) = ∂h(x) as in Definition 1. If h is differentiable, ∂h(x + ) = {∇h(x + )}.
The next remark simplifies several parts in the analysis of our algorithm.

Remark 5. From Assumption (M), the function J is continuous at any point (x, y) ∈ int(dom J) and obeys J(x, y) = J(x, y) where the continuous function J is given by J(x, y) := F (x) + G(y) + H(x, y).

We can note that this holds true for the assumptions in [START_REF] Attouch | On the convergence of the proximal algorithm for nonsmooth functions involving analytic features[END_REF][START_REF] Xu | A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion[END_REF][START_REF]A globally convergent algorithm for nonconvex optimization based on block coordinate update[END_REF]. For y fixed, the partial subdifferential of J(•, y) at x is denoted by ∂ x J(x, y); for x fixed, ∂ y J(x, y) is defined in a similar way. Proposition 1. Let J obey Assumption (M). Then, for any (x, y) ∈ dom J,

∂J(x, y) = ∂ x J(x, y) × ∂ y J(x, y) = (∇F (x) + ∂ x H(x, y)) × (∇G(y) + ∂ y H(x, y))
where the symbol "×" stands for Cartesian product.

Proof . From (9), the function J(x, y) := F (x) + G(y) + H(x, y) is differentiable. This shows that for any (x, y) ∈ dom J one has ∂J(x, y) = ∇ J(x, y) + ∂ χ Dx (x) + χ Dy (y) .

Using subdifferential calculus for separable functions [33, Prop. 10.5, p. 426 ] yields

∂ χ Dx (x) + χ Dy (y) = ∂χ Dx (x) × ∂χ Dy (y).
Since ∇ J(x, y) = (∇F (x) + ∇ x H(x, y), ∇G(y) + ∇ y H(x, y)), and using the last equality together with Remark 5, leads to the stated result.

The Fermat's rule, extended to nonconvex/nonsmooth functions is given next. For nonconvex functions, the Fermat's rule is only a necessary conditions. First-order methods can find only points satisfying this rule. Definition 3. We say that (x * , y * ) is a critical point of J if (0, 0) ∈ ∂J(x * , y * ).

The set of the critical points of J will be denoted by crit(J).

Two facts on limit point sets

Any bounded sequence {z k } k∈N has a convergent subsequence; thus the set of all its limit points

L(z 0 ) := z * | ∃ {k j } j∈N strictly increasing such that z k j → z * as j → ∞ . ( 11
)
is nonempty and closed. The next two claims are proven in Appendix.

Lemma 1. Let {z k } k∈N be a sequence of R m such that z k+1 -z k → 0 as k goes to infinity. Suppose that the set L(z 0 ) of the limit points of {z k } k∈N is nonempty and bounded. Then lim k→∞ dist z k , L(z 0 ) = 0.

An immediate useful consequence is the following:

Corollary 1. Let {z k } k∈N be a sequence of R m such that z k+1 -z k → 0 as k goes to infinity. The following statements are equivalent:

(a) the set L(z 0 ) of the limit points of {z k } k∈N is nonempty and bounded;

(b) the sequence {z k } k∈N is bounded.

Elements of proximal gradient descent

Proximity operators were inaugurated by Moreau in 1962 [START_REF] Moreau | Proximité et dualité dans un espace hilbertien[END_REF] as a generalization of convex projection operators. They were studied in numerous works, see the monograph of Bauschke and Combettes [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF] for an overview in the convex setting. Extension to non-convex functions was introduced in [START_REF] Attouch | Proximal alternating minimization and projection methods for nonconvex problems: An approach based on the Kurdyka-Lojasiewicz inequality[END_REF] and used by many authors afterwards [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF][START_REF] Pock | Inertial proximal alternating linearized minimization (iPALM) for nonconvex and nonsmooth problems[END_REF][START_REF] Chouzenoux | A block coordinate variable metric forwardbackward algorithm[END_REF][START_REF]A globally convergent algorithm for nonconvex optimization based on block coordinate update[END_REF].

Definition 4. Let h : R m → R ∪ {+∞} be a proper, lower semicontinuous function. Given z ∈ R m and τ > 0, the proximity operator of h at z ∈ R m is defined as

prox τ h (z) = arg min x∈R m h(x) + 1 2τ z -x 2 . ( 12
)
Remark 6. Some important consequences of this definition are given below:

(a) for any τ > 0, the set prox τ h (z) is nonempty and compact [33, sec. 7]. (b) if h is convex, the minimizer prox τ h (z) is uniquely defined as being the minimizer of a strictly convex coercive function;

(c) if h is twice differentiable, it is uniquely defined if inf dom h ∂ 2 h < 1 τ . (d) for any z ∈ R m one has prox τ h (z) ∈ dom h, which follows from Fermat's rule.
For h convex, a fundamental inequality associated with the proximity operator is given next.

Lemma 2 (Proximal inequality). Let

h : R m → R ∪ {+∞} be a proper, lower semicontinuous convex function. Given z ∈ R m and τ > 0, if x + = prox τ h (z), then h(x) ≥ h(x + ) - 1 τ x -x + , x + -z ∀ x ∈ R m .
The proximal gradient descent, known also as proximal forward-backward, has a crucial role when minimising the sum of a convex and a smooth function.

Definition 5. Let h : R m → R ∪ {+∞} be a proper, lower semicontinuous and function and f : R m → R be a differentiable function with L ∇f -Lipschitz continuous gradient. The main iteration to minimize f + h using the proximal gradient method starting from any u ∈ R m is given by

x + ∈ prox τ h (u -τ ∇f (u)) . (13)
The point

x + satisfies x + ∈ dom h, see Remark 6. It is uniquely defined if inf dom h ∂ 2 h < 1 τ . If h is convex then x + is uniquely defined.
Remark 7. The proximal gradient descent can be seen as the minimization of h + f , where f is a quadratic approximation of f around the point u:

x + ∈ prox τ h (u -τ ∇f (u)) = arg min x∈R m h(x) + 1 2τ x -(u -τ ∇f (u)) 2 = arg min x∈R m h(x) + f (u) + x -u, ∇f (u) + 1 2τ x -u 2 =: f (z) . (14)
We recall a classical result whose proof can be found in the textbook [12, A. 24].

Lemma 3 (Descent lemma). Let f : R m → R be a differentiable function with L ∇f -Lipschitz continuous gradient. Then

f (u) ≥ f (x) -x -u, ∇f (u) - L ∇f 2 x -u 2 ∀ x, u ∈ R m . ( 15
)
The next lemma warrants a sufficient decrease of the objective after a proximal step.

Lemma 4. Let f : R m → R be a differentiable function with L ∇f -Lipschitz continuous gradient and h : R m → R ∪ {+∞} a lower semicontinuous and proper function. For τ > 0 fixed, let x + be defined by [START_REF] Bierstone | Semianalytic and subanalytic sets[END_REF].

(a) A general fact

f (u) + h(u) ≥ f (x + ) + h(x + ) + 1 2 1 τ -L ∇f x + -u 2 ∀ u ∈ R m (16) (b) If h is convex, then f (u) + h(u) ≥ f (x + ) + h(x + ) + 1 τ - L ∇f 2 x + -u 2 ∀ u ∈ R m (17)
Remark 8. Using the notation Q = 2 if H is biconvex and Q = 1 otherwise, Lemma 4 can be summarized as:

f (u) + h(u) ≥ f (x + ) + h(x + ) + 1 2 Q τ -L ∇f x + -u 2 ∀ u ∈ R m
Proof . The descent Lemma 3 applied to f with x := x + yields

f (u) ≥ f (x + ) -x + -u, ∇f (u) - L ∇f 2 x + -u 2 (18)
(a) The result can be found in [17, Lemma 2]. It is recalled for better readability. From the minimization over x in [START_REF] Bochnak | Real Algebraic Geometry[END_REF], one has

h(x) + x -u, ∇f (u) + 1 2τ x -u 2 ≥ h(x + ) + x + -u, ∇f (u) + 1 2τ x + -u 2 ∀ x ∈ R m
Choosing x = u and rearranging leads to

h(u) ≥ h(x + ) + x + -u, ∇f (u) + 1 2τ x + -u 2 (19)
Summing this with (18) entails ( 16). (b) Apply the proximal inequality to h with x := u and z := u -τ ∇f (u):

h(u) ≥ h(x + ) - 1 τ u -x + , x + -u + τ ∇f (u) = h(x + ) + 1 τ x + -u 2 + x + -u, ∇f (u) .
Then, adding equation [START_REF] Chen | Smoothing methods for nonsmooth, nonconvex minimization[END_REF] to the result obtained above leads to [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF].

4 A simple algorithm for block regularized objectives 4.1 Alternating Structure-Adapted Proximal gradient descent algorithm (ASAP)

According to Remark 7, the iterations of our algorithm sketched in subsection 1.3, see [START_REF] Auslender | Asymptotic properties of the Fenchel dual functionaland applications to decomposition problems[END_REF], are equivalent to minimizing the following quadratic approximations of F around the point x k-1 and of G around the point y k-1 , respectively:

x k = arg min x∈U F (x k-1 ) + x -x k-1 , ∇F (x k-1 ) + H(x, y k-1 ) + 1 2τ x -x k-1 2 , ( 20
)
y k = arg min y∈V G(y k-1 ) + y -y k-1 , ∇G(y k-1 ) + H(x k , y) + 1 2σ y -y k-1 2 . ( 21
)
Hence, ignoring the terms that are constant at iteration k, namely F (x k-1 ), x k-1 , ∇F (x k-1 ) , G(y k-1 ) and y k-1 , ∇G(y k-1 ) , and using the definition of proximal gradient descent (Definition 5), the proposed ASAP algorithm takes the compact form stated below. We recall that Q ∈ {1, 2} is defined in [START_REF] Beck | Smoothing and first order methods: a unified framework[END_REF].

Algorithm 1 ASAP: Alternating Structure-Adapted Proximal gradient descent Initialization: (x 0 , y 0 ) ∈ U × V and 0 < τ < Q/L ∇F , 0 < σ < Q/L ∇G General Step: for k = 1, 2, . . ., compute x k = prox τ H(•,y k-1 ) x k-1 -τ ∇F (x k-1 ) , (22)
y k = prox σH(x k ,•) y k-1 -σ∇G(y k-1 ) . ( 23
)
The step-sizes τ and σ are set according to Lemma 4 so that they ensure a sufficient decrease of the objective J at each step. They depend only on the Lipschitzians of ∇F and ∇G, and Q. One can always set Q = 1 but when H is biconvex Q = 2 gives twice larger step-sizes. According to the initialization, if F = 0 (resp., G = 0), then τ (resp., σ) is a positive number.

Remark 9. The ASAP algorithm confirms immediately two attractive points that are direct consequences of Remark 6:

(a) for any k ≥ 1, it holds that (x k , y k ) ∈ dom J; (b) when H is multiconvex, iterates (x k , y k ) are always uniquely defined, even if F or G are nonconvex and step-sizes are twice larger.

ASAP algorithm with N blocks

The ASAP algorithm can be applied to the case when the problem has N > 2 blocks. Then the variable z is split into N blocks z = z (1) , . . . , z (N ) and the optimization problems has the form

J z (1) , . . . , z (N ) := H z (1) , . . . , z (N ) + N i=1 F i z (i) , H z (1) , . . . , z (N ) = H z (1) , . . . , z (N ) + N i=1 χ D i (z (i) ). ( 24
)
Analogously to Assumption (M), J is lower-bounded, and for any i = 1, . . . , N , functions

F i obeys (M)(b) With Q = 2 if
H is multiconvex and Q = 1 otherwise, the corresponding algorithm is: If the update is sequential, the alternating system has the form

(z k-1 (1) , • • • , z k-1 (N ) ) → (z k (1) , z k-1 (2) , • • • , z k-1 (N ) ) → • • • → (z k (1) , • • • , z k (N ) ). ( 25
)
Other block updates can also be used; in [39, sec. 3] a random shuffling was shown to provide better numerical results for a scheme of the form (5).

Algorithm 2 ASAP for N block problems Initialization: Updating rule (UR), z 0 (i) i and 0 < τ i < Q/L ∇F i for i = 1, . . . , N General Step: for k = 1, 2, . . . pick i ∈ {1, 2, . . . , N } according to UR and compute

z k (i) = prox τ i H(z k (1) ,...,z k (i-1) ,•,z k-1 (i+1) ,...,z k-1 (N ) ) z k-1 (i) -τ i ∇F i (z k-1 (i) ) .
5 Convergence analysis of ASAP for two block problems

Essential convergence facts

Given an initial (x 0 , y 0 ), the alternating system we have to study is of the form (x k-1 , y k-1 ) → (x k , y k-1 ) → (x k , y k ). Our first result states the convergence of the sequence {J(x k , y k )} k∈N to a real number J * with a guaranteed decrease.

Proposition 3. Let Assumption (M) hold and let {(x k , y k )} k∈N be a sequence generated by ASAP. Based on the step-size initialization of the algorithm we set

ρ := 1 2 min Q τ -L ∇F , Q σ -L ∇G > 0, ( 26
)
where Q ∈ {1, 2} is defined in [START_REF] Beck | Smoothing and first order methods: a unified framework[END_REF].

(a) For every k ≥ 1 the following sufficient decrease property holds:

J(x k-1 , y k-1 ) ≥ J(x k , y k-1 ) + ρ x k -x k-1 2 ; J(x k , y k-1 ) ≥ J(x k , y k ) + ρ y k -y k-1 2 ; J(x k-1 , y k-1 ) ≥ J(x k , y k ) + ρ x k -x k-1 2 + y k -y k-1 2 ; ( 27
)
and hence J(x k-1 , y k-1 ) ≥ J(x k , y k-1 ) ≥ J(x k , y k ).

(b) The sequences {J(x k , y k )} k∈N and {J(x k , y k-1 )} k∈N are non-increasing and interlaced by [START_REF] Hintermüller | Nonconvex TV q -models in image restoration: Analysis and a trustregion regularization-based superlinearly convergent solver[END_REF], hence they converge to a value denoted by J * ;

(c) We have

+∞ k=1 y k -y k-1 2 + x k -x k-1 2 < +∞
and hence lim k→∞ x k+1 -x k = 0 and lim k→∞ y k+1 -y k = 0.

Proof . Applying Lemma 4 with f := F , h := H(•, y k-1 ), x + := x k and x := x k-1 , along with the definition of ρ in ( 26), shows that

J(x k-1 , y k-1 ) = F (x k-1 ) + G(y k-1 ) + H(x k-1 , y k-1 ) ≥ F (x k ) + G(y k-1 ) + H(x k , y k-1 ) + ρ x k -x k-1 2 = J(x k , y k-1 ) + ρ x k -x k-1 2 . ( 29
)
Lemma 4 with f := G, h := H(x k , •), x + := y k and x := y k-1 , and (26) yield

J(x k , y k-1 ) = F (x k ) + G(y k-1 ) + H(x k , y k-1 ) ≥ F (x k ) + G(y k ) + H(x k , y k ) + ρ y k -y k-1 2 = J(x k , y k ) + ρ y k -y k-1 2 . ( 30
)
Bringing ( 29) and ( 30) together leads to

J(x k-1 , y k-1 ) ≥ J(x k , y k-1 ) + ρ x k -x k-1 2 ≥ J(x k , y k ) + ρ x k -x k-1 2 + ρ y k -y k-1 2 (31)
which completes the proof of (a). It follows that the sequences {J(x k , y k )} k and {J(x k , y k-1 )} k are nonincreasing and interlaced by [START_REF] Hintermüller | Nonconvex TV q -models in image restoration: Analysis and a trustregion regularization-based superlinearly convergent solver[END_REF], and bounded from below because J is lower bounded. Therefore, they converge to the same finite number J * , which proves (b). Using the inequalities in [START_REF] Moreau | Proximité et dualité dans un espace hilbertien[END_REF] we also have

∀ k ≥ 1, J(x k-1 , y k-1 ) -J(x k , y k ) ≥ ρ x k -x k-1 2 + y k -y k-1 2 . ( 32
)
For K ≥ 1, summing (32) from k = 1 to K yields

∀ K ≥ 1, J(x 0 , y 0 ) -J(x K , y K ) ≥ ρ K k=1 y k -y k-1 2 + x k -x k-1 2 .
Statement (b) entails that

∀ K ≥ 1, 1 ρ J(x 0 , y 0 ) -J * ≥ K k=1 y k -y k-1 2 + x k -x k-1 2 .
Taking the limit as K → ∞ leads to

+∞ k=1 y k -y k-1 2 + x k -x k-1 2 ≤ 1 ρ J(x 0 , y 0 ) -J * ,
which establishes statement (c).

From the proof of Proposition 3(c) we can derive a global O(1/k) convergence rate for x k -x k-1 2 + y k -y k-1 2 k∈N . Corollary 2 (Convergence rate). Let Assumption (M) hold and {(x k , y k )} k∈N be a sequence generated by ASAP. Then for any K ∈ N it holds that

inf k≥K x k -x k-1 2 + y k -y k-1 2 ≤ 1 ρK J(x 0 , y 0 ) -J * .
The following result, together with Proposition 3(c), shows that the partial subgradients of J vanish when k goes to infinity. Proposition 4. Let Assumption (M) hold and let {(x k , y k )} k∈N be a sequence generated by ASAP. Then for any k ≥ 1, one has

(a) p k x ∈ ∂ x J(x k , y k-1
) and q k y ∈ ∂ y J(x k , y k ) given by

p k x := ∇F (x k ) -∇F (x k-1 ) + 1 τ (x k-1 -x k ) ∈ ∂ x J(x k , y k-1 ), (33)
q k y := ∇G(y k ) -∇G(y k-1 ) + 1 σ (y k-1 -y k ) ∈ ∂ y J(x k , y k ), (34)
such that

p k x ≤ L ∇F + 1 τ x k-1 -x k and q k y ≤ L ∇G + 1 σ y k-1 -y k . (35)
(b) Further, (q k x , q k y ) ∈ ∂J(x k , y k ) for q k y as in [START_REF] Shiota | Geometry of subanalytic and semialgebraic sets[END_REF] and

q k x = ∇ x H(x k , y k ) -∇ x H(x k , y k-1 ) + p k x . (36)
Proof . (a) The Fermat's rule for x k in (20) yields

1 τ (x k-1 -x k ) ∈ ∇F (x k-1 ) + ∂ x H(x k , y k-1 ). ( 37
)
From the original formula for J in (8) one has

∂ x J(x k , y k-1 ) = ∇F (x k ) + ∂ x H(x k , y k-1 ).
Subtracting [START_REF] Xu | Alternating proximal gradient method for sparse nonnegative Tucker decomposition[END_REF] from this expression yields the result in [START_REF] Rockafellar | Variational analysis[END_REF]. Similarly, using Fermat's rule for y k in ( 21) -the y-step of the algorithm -and (8) leads to [START_REF] Shiota | Geometry of subanalytic and semialgebraic sets[END_REF].

Using that ∇F has Lipschitzian L ∇F (see Assumption (M)(b)), it follows that

p k x ≤ L ∇F x k-1 -x k + 1 τ x k-1 -x k = L ∇F + 1 τ x k-1 -x k ,
while the Lipschitz-continuity of ∇G leads to the bound on q k y . (b) From the equivalent formula for J in [START_REF] Beck | An alternating semiproximal method for nonconvex regularized structured total least squares problems[END_REF],

p k x = ∇F (x k ) + ∇ x H(x k , y k-1 ) + g k x ∈ ∂ x J(x k , y k-1 )
, where by [START_REF] Rockafellar | Variational analysis[END_REF] in Proposition 4 one has

g k x = -∇F (x k-1 ) + 1 τ (x k-1 -x k ) -∇ x H(x k , y k-1 ) ∈ ∂χ Dx (x k ).
This result, together with ( 9) and shows that q k x ∈ ∂ x J(x k , y k ) reads as

q k x := ∇ x H(x k , y k ) + ∇F (x k ) + g k x = ∇ x H(x k , y k ) + ∇F (x k ) -∇F (x k-1 ) + 1 τ (x k-1 -x k ) -∇ x H(x k , y k-1 )
Inserting ( 33) into the last expression proves [START_REF] Tseng | Convergence of a block coordinate descent method for nondifferentiable minimization[END_REF].

The set of all limit points of a sequence {(x k , y k )} k∈N generated by ASAP starting from a point (x 0 , y 0 ) is denoted by L(x 0 , y 0 ); see [START_REF] Beck | On the convergence of block coordinate descent type methods[END_REF] in subsection 3.2.

Remark 10. According to Corollary 1, assuming that L(x 0 , y 0 ) is nonempty and bounded is equivalent to assume that the iterates {(x k , y k )} k∈N are bounded.

Boundedness of the iterates is a usual assumption, see e.g., [START_REF] Attouch | Proximal alternating minimization and projection methods for nonconvex problems: An approach based on the Kurdyka-Lojasiewicz inequality[END_REF][START_REF] Xu | A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion[END_REF][START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF][START_REF] Hesse | Proximal heterogeneous block implicit-explicit method and application to blind ptychographic diffraction imaging[END_REF][START_REF] Beck | An alternating semiproximal method for nonconvex regularized structured total least squares problems[END_REF]. This assumption holds for instance when the level sets of J are bounded or if J is coercive.

Below we summarize several facts on the limit point set of ASAP.

Proposition 5. Let Assumption (M) hold and let {(x k , y k )} k∈N be a sequence generated by ASAP which is assumed to be bounded. Let (x * , y * ) ∈ L(x 0 , y 0 ).

(a) there is a subsequence

x k j , y k j j∈N such that (x k j , y k j ) → (x * , y * ) as j → ∞; (b) lim k→∞ J(x k , y k ) = J(x * , y * ); (c) (0, 0) ∈ ∂J(x * , y * ) and thus (x * , y * ) is a critical point of J.
Proof . (a) follows from the definition of L(x 0 , y 0 ). (b) From Remark 9(a), (x k j , y k j ) belongs to dom J. From Remark 5, J(x k j , y k j ) = J(x k j , y k j ) where J is a continuous function. Therefore

lim j→∞ J(x k j , y k j ) = lim j→∞ J(x k j , y k j ) = J(x * , y * ).
Since dom J is closed one has (x * , y * ) ∈ dom J and J(x * , y * ) = J(x * , y * ). From Proposition 3(a)-(b) the sequence J(x k , y k ) k∈N is nonincreasing, interlaced by [START_REF] Hintermüller | Nonconvex TV q -models in image restoration: Analysis and a trustregion regularization-based superlinearly convergent solver[END_REF] and converges to J * . This, together with the above results completes statement (b).

(c) The cases when H satisfies (M)(c)(i) and (M)(c)(ii) need to be considered separately.

-Case (M)(c)(i). Using Proposition 4 we have (q

k j x , q k j y ) ∈ ∂J(x k j , y k j ) where q k j x = ∇ x H(x k j , y k j ) -∇ x H(x k j , y k j -1 ) + p k j x and q k j y = ∇G(y k j ) -∇G(y k j -1 ) + 1 σ (y k j -1 -y k j ) Since ∇ x H is continuous, ∇ x H(x k j , y k j ) -∇ x H(x k j , y k j -1
) → 0 as j → ∞. From [START_REF] Soncco | Fast and accurate multiplicative decomposition for fringe removal in interferometric images[END_REF] in Proposition 4 we have (p k j

x ) j → 0. These facts show that q k j x → 0. From Proposition 4, q k j y → 0. Therefore, (q Remark 4) this shows that (0, 0) ∈ ∂J(x * , y * ).

k j x , q k j y ) → 0 as k → ∞. Since ∂J is closed (see
-Case (M)(c)(ii) -H is biconvex. Since H is biconvex, The formula for J in ( 8), together with ( 33) and [START_REF] Shiota | Geometry of subanalytic and semialgebraic sets[END_REF] in Proposition 4 yield

p k j x -∇F (x k j ) ∈ ∂ x H(x k j , y k j -1 ) and q k j y -∇G(y k j ) ∈ ∂ y H(x k j , y k j ).
x → H(x, y k j -1 ) and y → H(x k j , y) are convex, lower semicontinuous and proper functions (see Assumption (M)(c)) the subgradient inequality (Definition 1) can be applied which leads to

∀ x ∈ U, H(x, y k j -1 ) ≥ H(x k j , y k j -1 ) + p k j x -∇F (x k j ), x -x k j , ∀ y ∈ V, H(x k j , y) ≥ H(x k j , y k j ) + q k j y -∇G(y k j ), y -y k j .
From [START_REF] Soncco | Fast and accurate multiplicative decomposition for fringe removal in interferometric images[END_REF] in Proposition 4 we have (p k j x ) j → 0 and (q k j y ) j → 0. By Proposition 3(c), (y k j -1 -y k j ) j∈N converges to 0. Since {y k j } j∈N converges to y * , it follows that {y k j -1 } j∈N converges to y * as well. Using the facts that ∇F , ∇G and H = H are continuous functions on int(dom J) and that ∂J is closed (see Assumption (M) and Remark 5) we can evaluate the limit in the inequalities above when j → +∞:

∀ x ∈ U, H(x, y * ) ≥ H(x * , y * ) + 0 -∇F (x * ), x -x * , ∀ y ∈ V, H(x * , y) ≥ H(x * , y * ) + 0 -∇G(y * ), y -y * .
Then, by the definition of the (partial) subgradients of H (see Definition 1), 0 ∈ ∇F (x * ) + ∂ x H(x * , y * ) and 0 ∈ ∇G(y * ) + ∂ y H(x * , y * ).

Invoking Proposition 1, it holds that (0, 0) ∈ ∂J(x * , y * ).

Proposition 6. Let Assumption (M) hold and let {(x k , y k )} k∈N be a sequence generated by ASAP which is assumed to be bounded. Then the following properties hold: (d) The sequences J(x k , y k-1 ) and J(x k , y k ) decrease and converge to a finite value J * by Proposition 3.

(a) L(x 0 , y 0 ) ⊂ crit(J); (b) lim k→∞ dist (x k , y k ), L(x 0 , y 0 ) = 0; (c) L(x 0 , y 0 ) is connected; (d) J is
According to (a), there may be critical points of J that cannot be reached when ASAP starts from a given initial (x 0 , y 0 ) which emphasizes the role of this initial guess. Statement (b) guarantees that for k large enough, the iterates (x k , y k ) are arbitrary close to a critical point of J. We conclude this subsection with the following remark: Remark 11. Proposition 5 ensures only subsequential convergence of ASAP to critical points of J. In practice this result might be enough thanks to Proposition 6. Yet, if J shares some special properties (e.g., all critical points are isolated), one may be able to conclude with convergence of the whole sequence (x k , y k ) towards a critical point of J. For an example, see [START_REF] Xu | A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion[END_REF]Cor. 2.4].

Subgradient convergence

Subgradient convergence can be proved if one considers an additional assumption on the regularity of the smooth part H of the coupling term H (see Assumption (M)(c) and ( 9)):

Assumption (H) ∇ H is one-sided locally Lipschitz continuous on bounded subsets of U × V in the sense that for each bounded subset

B U × B V ⊂ U × V there is a constant ξ > 0 such that ∀ (x, y), (x , y ) ∈ B U × B V , ∇ x H(x, y) -∇ x H(x, y ) ≤ ξ y -y .
A stronger assumption, namely that ∇ H is globally Lipschitz continuous on bounded subsets is often used in other papers, see e.g., [START_REF] Attouch | Proximal alternating minimization and projection methods for nonconvex problems: An approach based on the Kurdyka-Lojasiewicz inequality[END_REF][START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF][START_REF] Pock | Inertial proximal alternating linearized minimization (iPALM) for nonconvex and nonsmooth problems[END_REF][START_REF] Beck | An alternating semiproximal method for nonconvex regularized structured total least squares problems[END_REF][START_REF] Chouzenoux | A block coordinate variable metric forwardbackward algorithm[END_REF]. Assumption (H) can also be derived if H is twice continuously differentiable.

Proposition 7. Let Assumptions (M) and (H) hold. Let {(x k , y k )} k∈N be a sequence generated by ASAP which is assumed to be bounded. Then there exists ξ ∈ (0, ∞) such that for any k ≥ 1 one has

∃ q k x , q k y ∈ ∂J(x k , y k ) obeying (q k x , q k y ) ≤ β x k -x k-1 , y k -y k-1 ,
where

β := max    √ 2 L ∇F + 1 τ , L ∇G + 1 σ 2 + 2 ξ 2    .
Proof . Since the iterates are bounded, there exist

B U × B V ⊂ U × V bounded such that (x k , y k ) ∈ B U × B V for any k ∈ N.
Using assumption (H), the constant ξ below is finite:

∀ k ∈ N, ∇ x H(x k , y k ) -∇ x H(x k , y k-1 ) ≤ ξ y k -y k-1 . (38)
Using Proposition 4(a) and (b) yields

q k j x = ∇ x H(x k j , y k j ) -∇ x H(x k j , y k j -1 ) + p k j x
and q

k j y = ∇G(y k j ) -∇G(y k j -1 ) + 1 σ (y k j -1 -y k j )
Using [START_REF] Xu | A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion[END_REF] and the first inequality in [START_REF] Soncco | Fast and accurate multiplicative decomposition for fringe removal in interferometric images[END_REF] in Proposition 4 yields

q k x ≤ ∇ x H(x k , y k ) -∇ x H(x k , y k-1 ) + p k x ≤ ξ y k -y k-1 + L ∇F + 1 τ x k -x k-1 .
Thus (q k x , q k y ) ∈ ∂J(x k , y k ) where q k y is given in [START_REF] Shiota | Geometry of subanalytic and semialgebraic sets[END_REF]. The obtained result on q k x together with the second inequality in [START_REF] Soncco | Fast and accurate multiplicative decomposition for fringe removal in interferometric images[END_REF] in Proposition 4, and using that (a + b) 2 ≤ 2 a 2 + 2 b 2 , shows that

(q k x , q k y ) 2 = q k x 2 + q k y 2 ≤ ξ y k -y k-1 + L ∇F + 1 τ x k -x k-1 2 + L ∇G + 1 σ 2 y k -y k-1 2 ≤ 2 L ∇F + 1 τ 2 x k -x k-1 2 + L ∇G + 1 σ 2 + 2ξ 2 y k -y k-1 2 .
This gives the value of β in the statement which completes the proof.

5.3

The Kurdyka-Lojasiewicz (K L) property and the objective J

The Kurdyka-Lojasiewicz (K L) property was studied for the optimization of nonsmooth functions originally in [START_REF] Bolte | The lojasiewicz inequality for nonsmooth subanalytic functions with applications to subgradient dynamical systems[END_REF] and further on in [START_REF] Attouch | On the convergence of the proximal algorithm for nonsmooth functions involving analytic features[END_REF][START_REF] Attouch | Proximal alternating minimization and projection methods for nonconvex problems: An approach based on the Kurdyka-Lojasiewicz inequality[END_REF][START_REF] Attouch | Convergence of descent methods for semi-algebraic and tame problems: proximal algorithms, forward-backward splitting, and regularized Gauss-Seidel methods[END_REF].

Definition 6 (K L property). Let f : R m → R ∪ {+∞} be a proper lower semicontinuous function. The function f is said to have the Kurdyka-Lojasiewicz (K L) property at x * ∈ dom ∂f if there exist η ∈ (0, +∞], a neighborhood O(x * ) of x * and a constant κ > 0 such that

∀ x ∈ O(x * ), κ dist(0, ∂f (x)) ≥ |f (x) -f (x * )| θ , (39)
where θ ∈ [0, 1) and

O(x * ) := O(x * ) ∩ {x ∈ U | f (x * ) < f (x) < f (x * ) + η} . (40)
Observe that the K L property implies that all critical points belonging to O(x * ) share the same critical value f (x * ) since otherwise [START_REF]A globally convergent algorithm for nonconvex optimization based on block coordinate update[END_REF] would fail. The K L property requires assumptions only on the shape of the function around its critical points. The K L property does not require that the critical points are strict or connected. We recall that if f is a real analytic function on O(x * ), then the K L property holds with θ ∈ [1/2, 1) thanks to the Lojasiewicz gradient inequality [30, p. 92].

It is instructive to provide some examples and counter-examples.

Example 2.

(a) The function f : R → R given by f

(x) = (x -1) 2 /2 satisfies |f (x)|/ √ 2 = |f (x) -f (x * )| 1/2
for any x ∈ R \ {x * } where x * = 1 is its minimizer. Thus f satisfies the K L property at x * with θ = 1/2 and κ = √ θ. (b) The function h(x, y) = (xy -1) 2 is analytic, has a nonstrict minimizer at (x * , y * ) = (1, 1), and satisfies the K L property at (x * , y * ). (c) The functions below are infinitely differentiable but fail the K L property:

f (x) = exp -1/x 2 if x > 0, 0 if x ≤ 0, f (x) = x 2 sin(1/x) if x = 0, 0 if x = 0.
The first one is "infinitely flat" [5, p. 453] while the second presents "wild oscillations" [16, p. 569].

We recall the following fundamental fact:

Remark 12. Both real-analytic and semi-algebraic functions are semi-analytic functions and thus they are subanalytic functions, according to [START_REF] Bierstone | Semianalytic and subanalytic sets[END_REF].

It is quite obvious that, in general, semi-algebraic functions are not real-analytic (one could think of functions that are piecewisely built from real-analytic functions such as Example 1 (ii)). Conversely, some real-analytic functions are not semi-algebraic; this is the case of the exponential function. Functions that are both real-analytic and semi-algebraic are called Nash functions [START_REF] Bochnak | Real Algebraic Geometry[END_REF]; see, e.g., functions ψ in (i) and (iii) in Example 1.. We will use the result obtained by Bolte, Daniilidis, and Lewis in [START_REF] Bolte | The lojasiewicz inequality for nonsmooth subanalytic functions with applications to subgradient dynamical systems[END_REF] saying that all subanalytic functions enjoy the K L property at their critical points: Theorem 1. [15, Theorem 3.1] Let f : R m → R ∪ {+∞} be a function which is subanalytic with a closed domain and continuous on its domain. Let x * be a critical point of f . Then there is an exponent θ ∈ [0, 1), a neighborhood O(x * ) of x * and a constant κ > 0 such that

∀ x ∈ O(x * ), κ dist(0, ∂f (x)) ≥ |f (x) -f (x * )| θ . (41)
The objectives J we suggest in subsection 2.2 contain various functions and involves compositions and sums of them.

The following properties of subanalytic functions are taken from the monograph of Shiota [START_REF] Shiota | Geometry of subanalytic and semialgebraic sets[END_REF] and from the work of Denkowska and Denkowski [START_REF] Denkowska | A long and winding road to definable sets[END_REF], respectively. Proposition 8. Let f and g be two subanalytic functions. Then the following results hold:

(a) [START_REF] Shiota | Geometry of subanalytic and semialgebraic sets[END_REF]Chapter II.1] If f and g are lower-bounded, then f + g is subanalytic.

(b) [START_REF] Denkowska | A long and winding road to definable sets[END_REF]Proposition 2.46] If g maps bounded sets on bounded sets or if f -1 (X ) is bounded for any bounded subset X , then f • g is a subanalytic function.

Thanks to these results we can derive some useful properties of the objectives J in subsection 2.2. We begin with the terms composing J in (2.2). (b) According to Proposition 9, all functions x → f i ( A i x ), y → g j ( B i y ) and (x, y) → h( b(x, y) -w ) are subanalytic since f i , g j and h are subanalytic functions by (a). Thus J is a finite sum of lower-bounded subanalytic functions, hence J is subanalytic according to Proposition 8(a). Noticing also that J has a closed domain and is continuous on its domain, see Remark 5, Theorem 1 shows that J has the K L property at its critical points.

Convergence of ASAP to critical points under the K L property

We summarize that, under Assumptions (M) and (H), we have proved that any bounded sequence generated by ASAP satisfy the assumptions in [6, Theorem 2.9]:

1) there exists ρ ∈ (0, ∞) such that for any k ≥ 1 the sufficient decrease property holds (Proposition 3(a)):

J(x k , y k ) + ρ x k -x k-1 2 + y k -y k-1 2 ≤ J(x k-1 , y k-1 );
2) there exists a subsequence (x k j , y k j ) and a critical point (x * , y * ) of J such that the subsequential continuity towards a critical point holds (Proposition 5): lim j→∞ x k j , y k j = (x * , y * ) and lim j→∞ J x k j , y k j = J(x * , y * );

3) there exists β ∈ (0, ∞) such that ∀ k ≥ 1 the subgradient relative error condition holds (Proposition 7):

∃ q k x , q k y ∈ ∂J(x k , y k ) obeying (q k x , q k y ) ≤ β (x k -x k-1 , y k -y k-1 ) .
Theorem 3. Let Assumptions (M) and (H) hold and let {(x k , y k )} k∈N be a sequence generated by ASAP which is assumed to be bounded. Assume also that J has the K L property at a limit point (x * , y * ) ∈ L(x 0 , y 0 ). Then the following hold:

(a) {(x k , y k )} k∈N is a Cauchy sequence that converges to a critical point (x * , y * ) of J as k goes to infinity;

(b) moreover, +∞ k=0 x k+1 -x k + y k+1 -y k < +∞.
Proof . This theorem is a direct consequence of [6, Theorem 2.9]. Our Definition 6 endowed with the function ϕ : R → R given by fulfills [6, Def. 2.4]. All conditions needed in [6, Theorem 2.9] are the three statements given before the theorem, along with the K L property as stated.

Theorem 3 proves that for objectives J satisfying Assumptions (M) and (H), any bounded sequence generated by ASAP converges to a limit point (x * , y * ) ∈ L(x 0 , y 0 ) provided J has the K L property at (x * , y * ). This includes various nonconvex objective functions J that are continuous on their closed domain. A large family of nonconvex objective functions J was given in subsection 2.2. More generally, the objective functions J can be composed out of differentiable components that are real-analytic and semi-algebraic (nonconvex) functions, subanalytic (nonconvex) functions, strongly and uniformly convex functions.

Convergence of ASAP with H multiconvex

The variable z is split into N blocks, i.e., z = (z (1) , . . . , z (N ) ) ∈ U 1 × • • • × U N with U i finite-dimensional real spaces, and the objective J reads as

J z (1) , . . . , z (N ) := N i=1 F i z (i) + N i=1 χ D i (z (i) ) + H z (1) , . . . , z (N ) =:H(z (1) ,...,z (N ) ) (42)
The objective and the full algorithm were presented in subsection 4.2 where the main assumptions were briefly described. In this section we adopt the following hypothesis: Assumption (U) The update in Algorithm 2 is sequential, see [START_REF] Grippo | On the convergence of the block nonlinear Gauss-Seidel method under convex constraints[END_REF].

Then all statements in subsections 5.1 and 5.2 have a straightforward extension to the case when J has N > 2 blocks, using analogous arguments and proofs, but at the expense of heavier notation. Only the main statements are presented.

The next assumption is a direct extension of Assumption (M): for some constants -∞ ≤ a x < b x ≤ +∞ and -∞ ≤ a y < b y ≤ +∞. We define the function h : R 2 → R by h(x i , y i ) = 1 2 (x i y i -w i ) 2 so that H(x, y) = mn i=1 h(x i , y i ). Then the proximity operators are computed componentwisely prox τ H(•,y) (z) = prox τ h(•,y i ) (z i ) mn i=1 where prox τ h(•,y i ) (z i ) = arg min

Assumption (M') (a) J : U 1 × • • • × U N → R ∪ {+∞}
x i ∈[ax,bx] 1 2 (x i y i -w i ) 2 + 1 2τ (x i -z i ) 2 ∀ i.
The computation of prox σH(x,•) (z) is done in a similar way. The full algorithms is particularly simple:

Algorithm 3 ASAP with Hadamard product Initialization: (x 0 , y 0 ) ∈ U × V and 0 < τ < 2/L ∇F , 0 < σ < 2/L ∇G General Step: for k = 1, 2, . . ., compute

x k = min max τ w • y k-1 + z τ (y k-1 • y k-1 + 1)
, a x , b x with z = x k-1 -τ ∇F (x k-1 );

y k = min max σw • x k + z σ(x k • x k + 1)
, a y , b y with z = y k-1 -σ∇G(y k-1 ).

Fringe Separation in Interferometric Images

Sieleters [START_REF] Coudrain | SIELETERS, an airborne infrared dual-band spectro-imaging system for measurement of scene spectral signatures[END_REF] is a nonconventional infrared spectro-imaging device, which is based on interferometric imaging. Its purpose is to obtain hyperspectral images from a temporal image sequence, provided the images can be accurately registered. Due to the imaging system, the images provided by Sieleters have so-called interference fringes, which have a particular structure (see Fig. 7.2 bottom left). However, for many applications (such as stereo-matching, which is necessary to compute accurate image registration), these fringes have to be removed, so that one can recover the so-called panchromatic images. The latter can then be handled by conventional image processing techniques.

According to [START_REF] Soncco | Fast and accurate multiplicative decomposition for fringe removal in interferometric images[END_REF], a multiplicative model for the fringe formation was shown to be much more physically accurate than previous additive decomposition models. According to this model, an observed image w ∈ R m×n is given by w = x • y + perturbations, where x is the panchromatic image to recover, y -the fringe oscillations.and 1 stand for an image composed of ones. Two facts known from the physics of the observation device are that -the fringe oscillations have a range constraint, |y i,j | ≤ 1;

-the 1D Fourier transform of the columns of the fringe y, denoted by F(y), belongs to a known interval.

Applying these constraints enables the obtention of a good initialization (x 0 , y 0 ). Other observations are that the gradients of the panchromatic image and of the fringes are not sparse and that the perturbations under the multiplicative model are negligible. we focus on an optimization problem of the form J(x, y) = We were given 9 "ground-truth" 424×1000 images computed from an ONERA's Sieleters image sequence using a protocol based on the physical model. Fig. 7.2 shows 140 × 250 zooms. Algorithm 3 was successfully applied on Sieleters sequences, which typically consist in thousands of such airborne images.

Proposition 2 .

 2 [START_REF] Rockafellar | Variational analysis[END_REF] Theorem 10.1] Let f : R m → R ∪ {+∞} be a proper function. If f has a local minimum at x * , then 0 ∈ ∂f (x * ).

  finite and constant on L(x 0 , y 0 ). Proof . (a) follows from Proposition 5(c). (b) is a direct consequence of Lemma 1 and of Proposition 3(c). (c) is a consequence of Proposition 3(c) and the previous claim (b).

Proposition 9 .Theorem 2 .

 92 Let f : R → R ∪ {+∞} be a subanalytic function and U a finite dimensional real space. Then: (a) If A is a linear operator on U then x → f ( Ax ) is a subanalytic function. (b) If (V, W ) are finite dimensional real spaces and b := (b (1) , . . . , b (N ) ) : U × V → W a bilinear form, then (x, y) → f ( b(x, y) -w ) is a subanalytic function. Proof . (a) Since • is a distance function, the application x → Ax is semi-algebraic and thus subanalytic (Remark 12), and it also maps bounded sets on bounded sets. Then its composition with f is subanalytic according to Proposition 8(b). (b) For any k = 1, . . . , N , b (n) : U ×V → R is a bilinear form. Since (U, V ) are of finite dimension, b (k) can be represented using a real matrix A k so that b (k) (x, y) = A n x, y . Hence, b (k) is real polynomial. Consequently, (x, y) → b(x, y) -w 2 = n (b (k) (x, y) -w (k) ) 2 is semi-algebraic, and so is (x, y) → b(x, y) -w since x → x 1/2 is a semi-algebraic function. Using Proposition 8(b), it follows that (x, y) → f ( b(x, y) -w ) is subanalytic. For the family of J in subsection 2.2, the following hold: (a) all functions in Example 1 are subanalytic; (b) the objective function J in (2.2) is subanalytic on its domain. Hence, J satisfy the K L property at its critical points. Proof . (a) Functions ψ in (i), (ii) and (iii) in Example 1 are semi-algebraic and thus subanalytic. Further, since p is rational, then all (ψ) p in (i)-(iii) are semi-algebraic by composition [17, Ex. 4] and hence subanalytic. Functions ψ in (v), (vi) and (vii) are real-analytic, and thus subanalytic. Since t → |t| maps bounded sets on bounded sets, and since t → log(1 + t/α) is real-analytic, t → α log(1 + |t|/α) is subanalytic according to Proposition 8(b) Thus, function (iv) is the sum of the lower-bounded semi-algebraic function t → |t| and a subanalytic function, hence it is subanalytic according to Proposition 8(a).

7. 1

 1 The ASAP algorithm with Hadamard product Here x, y and w are seen as mn-length vectors. We focus on constraints (D x , D y ) of the formD x := {x ∈ U | a x ≤ x i ≤ b x } and D y := {y ∈ U | a y ≤ x i ≤ b y }

2 F

 2 -1 M F(y) 2 =:G(y) + η 2 x • y -w 2 + χ Dy (y)=:H(x,y)

Figure 1 :

 1 Figure 1: Result with Algorithm 3 (zooms). Clockwise: data, recovered panchromatic image, recovered fringe, reference

  are continuously differentiable and their gradients ∇F and ∇G are Lipschitz continuous with constants L ∇F and L ∇G , respectively; (c) H : U × V → R ∪ {+∞} has a two block closed domain , and one of the following holds:

  is lower bounded; (b) for any i = 1, . . . , N , F i : U i → R has a gradient which is Lipschitz continuous with constant L ∇F i ; (c) H has an N -block closed domain and on its domain (i) H is continuously differentiable;

(ii) z (i) → H(z) is convex and differentiable. Proposition 10. Let Assumptions (M') and (U) hold and let {z k } k∈N be a sequence generated by ASAP in Algorithm 2. Then the following holds:

A particular form of this algorithm designed for a specific application was proposed by Soncco, Barbanson, Nikolova et al.[START_REF] Soncco | Fast and accurate multiplicative decomposition for fringe removal in interferometric images[END_REF] without convergence proof.

For instance, in section 7 an application with U = V = R m×n is presented.

for fixed y ∈ V , x → H(x, y) is convex, and for fixed x ∈ U , y → H(x, y) is convex.

Acknowledgements

The authors would like to thank Jérôme Bolte for helpful comments on Remark 12.

This work has been partially funded by the French Research Agency (ANR) under grant No ANR-14-CE27-001 (MIRIAM).

(a) for every k ≥ 1 one has J(z k-1 ) ≥ J(z k ) + ρ z k -z k-1 2 where ρ := 1 2 min i=1,...,N

(b) assuming that {z k } k∈N is bounded and denoting z * ∈ L(z 0 ) a limit-point of {z k } k∈N , one has lim k→∞ J(z k ) = J(z * ) and z * is a critical point of J.

Proof . (a) As in the proof of Proposition 3(a), applying Lemma 4 with f := F i and h := H(z k (1) , . . . , z k (i-1) , •, z k-1 (i+1) , . . . , z k-1 (N ) ) yields for any i = 1, . . . , N

The result is obtained by

is defined for any i = 1, . . . , N and for all k ≥ 1 by:

By Fermat's rule, for any i = 1, . . . , N it holds that

Then the claim is proven using the pipeline established in the proofs of Propositions 4 and 5.

As in the two block case, strong subgradient convergence can be stated using an additional assumption on the smooth part H of the coupling term:

Assumption (H') ∇ H is partially locally Lipschitz continuous on bounded subsets in the sense that for each bounded subset

The papers dealing with N block coupling terms use a slightly stronger assumption to prove subgradient convergence [START_REF] Xu | A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion[END_REF][START_REF] Hesse | Proximal heterogeneous block implicit-explicit method and application to blind ptychographic diffraction imaging[END_REF][START_REF] Chouzenoux | A block coordinate variable metric forwardbackward algorithm[END_REF][START_REF]A globally convergent algorithm for nonconvex optimization based on block coordinate update[END_REF]. Proposition 11. Let Assumptions (M'), (U) and (H') hold. Let {z k } k∈N be a sequence generated by ASAP (Algorithm 2) which is assumed to be bounded. Then there exists ξ ∈ (0, ∞) such that for any k ≥ 1 one has

where

and γ i := L ∇F i + 1/τ i for any i = 1, . . . , N .

Proof . Using that the iterates are bounded and Assumption (M'), for any k ∈ N and for any

Following the proof of Proposition 7, the idea is to find

Inserting the expression for g k i in the above equation shows that

Using Assumption (H') and the Lipschitz continuity of ∇F i , the last equation yields

. Then, using yet again that (a+b) 2 ≤ 2 a 2 +2 b 2 , the sum of q k (i) 2 over i = 1, . . . , N is computed. Its square root is the desired q k and β follows from the obtained inequality.

We have proven that ASAP for N blocks satisfies all the three conditions in [6, Theorem 2.9]: sufficient decrease by Proposition 10(a), subgradient relative error by Proposition 11 and subsequential continuity towards a critical point by Proposition 10(b). Therefore we can state the following result: Theorem 4. Let Assumptions (M'), (U) and (H') hold. Let {z k } k∈N be a sequence generated by ASAP (Algorithm 2) which is assumed to be bounded. Assume also that J has the K L property at a limit point (z * ) ∈ L(z 0 ). Then the following assertions hold:

(a) {(z k )} k∈N is a Cauchy sequence that converges to a critical point (z * ) of J as k goes to infinity;

Let us consider general objectives of form (see subsection 2.2)

where b :

are linear mappings and f i,j are as in Example 1. The results in subsection 5.3 show that J in (44) is also a subanalytic function and thus fulfils the K L property.

An application for Hadamard based coupling terms

Here we focus on objective functions J of the form presented in subsection 2.2 for U = V = W = R m×n and a bilinear mapping b given by b(x, y) = x • y where "•" denotes the Hadamard (componentwise) product.

Following Assumption (M), the coupling term H reads as

Hadamard matrix products arise in lossy compression, tensor factorization, image processing, statistics, among others.

Appendix 7.3 Proof of Lemma 1

For simplicity, set L := L(z 0 ). Since L is closed and nonempty, the distance dist(z k , L) of z k to the set L is well defined. Suppose that dist(z k , L) does not go to zero as k → ∞, i.e., that there exist ε > 0 and a subsequence {z k j } j such that

Observe that { z k+1 -z k } k∈N is a Cauchy sequence of non-negative numbers converging to zero. This, together with the fact that z → dist(z, L) is a continuous function [33, p. 19.], entails that for ε > 0 there is K ∈ N such that for any k > K one has dist(z k+1 , L) -dist(z k , L) < ε. Therefore, there exists K ∈ N so that for any j > K one has k j > K and thus

This, together with (45), yields

From the definition of L, see [START_REF] Beck | On the convergence of block coordinate descent type methods[END_REF], there is an infinite number of points of z k satisfying dist(z k , L) ≤ ε.

Hence there exists k j + 1 ≤ k j+1 such that dist(z k j +1 , L) < ε; thus, a contradiction.

Proof of Corollary 1

If L(z 0 ) is nonempty and bounded, there exists M > 0 such that any z * ∈ L(z 0 ) satisfies z * ≤ M . From Lemma 1, the distance of z k to L(z 0 ) goes to zero as k goes to infinity, hence there exists k 0 ∈ N such that k ≥ k 0 implies that min z∈L(z 0 ) z k -z < M . Therefore, for any k ≥ k 0 , one has z k ≤ min z∈L(z 0 ) z k -z + z * ≤ 2M , hence (a) ⇒ (b). The converse claim is obvious.