Pauline Tan

Cmla

Fabien Pierre

Mila Nikolova

Inertial Alternating Generalized Forward-Backward Splitting for Image Colorization

In this paper, we propose a novel accelerated alternating optimization scheme to solve block-biconvex nonsmooth problems whose objectives can be split into smooth (separable) regularizers and simple coupling terms. The proposed method performs a Bregman distance based generalization of the well-known forward-backward splitting for each block, along with an inertial strategy which aims at getting empirical acceleration. We discuss the theoretical convergence of the proposed scheme and provide numerical experiments on image colorization.

Introduction 1.Image Colorization

In the previous decade, the number of papers dealing with the challenging problem of image colorization has literally exploded. This technique consists in transforming a gray-scale image (also known as black-and-white image) into a color one by addition of a chromatic information. This topic has many opportunities in the entertainment industry with the success of the recolorized historical documentaries. The challenge of this technique comes from the fact that turning a color into a gray level is not an invertible operation. Indeed the color space has three dimensions whereas the gray level space has only one. In this way, turning a gray-scale image into a color one requires additional information.

In the literature, there are two ways to add color information to a gray scale image. First, the manual approaches require some color strokes given by a user, which are then propagated over the whole image with diffusion techniques guided by the grayscale channel (see, e.g., [START_REF] Levin | Colorization using optimization[END_REF][START_REF] Yatziv | Fast image and video colorization using chrominance blending[END_REF][START_REF] Horiuchi | Colorization algorithm for grayscale image by propagating seed pixels[END_REF]). The major drawback of this approach is the huge amount of work needed from the user, especially when the scene represented on the image is complex or contains a lot of textured areas.

To tackle this issue, the exemplar-based methods use a color image as a reference and transfer the color to the gray-scale one based on texture criterion (see, e.g., [START_REF] Welsh | Transferring color to greyscale images[END_REF][START_REF] Ironi | Colorization by example[END_REF][START_REF] Ha | Variational models for image colorization via chromaticity and brightness decomposition[END_REF][START_REF] Horiuchi | Colorization algorithm using probabilistic relaxation[END_REF][START_REF] Kumar Gupta | Image colorization using similar images[END_REF][START_REF] Chen | Grayscale image matting and colorization[END_REF][START_REF] Charpiat | Automatic image colorization via multimodal predictions[END_REF]), morphological mappings [START_REF] Persch | Exemplar-based face colorization using image morphing[END_REF] or based on deep learning approach [START_REF] Zhang | Colorful image colorization[END_REF]. Obviously, the choice of the reference image is a first issue and influences the results. To take into account the various aspects of the textures for a reliable comparison, it is common in the literature to use many descriptors to transfer the colors between images [START_REF] Welsh | Transferring color to greyscale images[END_REF][START_REF] Pierre | Exemplar-based colorization in RGB color space[END_REF]. The main difficulty is the choice of the descriptors for textures, as well as the distance between the image patches in the case of patch-based approach. This choice can been done experimentally [START_REF] Welsh | Transferring color to greyscale images[END_REF] or with metric learning approach [START_REF] Pierre | Combinaison linéaire optimale de métriques pour la colorisation d'images[END_REF].

The alternative consists in the computation of multiple candidates followed by a method to select the final result. These candidates can be extracted from a patch-based approach [START_REF] Bugeau | Patch-based image colorization[END_REF] or based on a learning approach [START_REF] Zhang | Colorful image colorization[END_REF]. Some state-of-the art methods perform the choice among these candidates without taking into account the geometrical information of the result. For instance, the authors of [START_REF] Bugeau | Patch-based image colorization[END_REF] use a median to choose one candidate among the ones obtained by their patch based approach. As another example, the authors of [START_REF] Zhang | Colorful image colorization[END_REF] use an annealing mean to compute the final result from the different results given by a Convolutional Neural Network (CNN).

To choose a candidate and to enhance the regularity of the result, a first approach is the method of [START_REF] Charpiat | Automatic image colorization via multimodal predictions[END_REF] which is based on a graph-cut approach to select the final candidate. The authors of [START_REF] Bugeau | Variational exemplarbased image colorization[END_REF] have proposed a variational method to choose one candidate among given ones with a regularity hypothesis. This method has been improved then for image colorization by the authors of [START_REF] Pierre | Luminance-chrominance model for image colorization[END_REF]. The aim of these last two models is to choose, for each pixel of the image, a candidate such that the total variation of the resulting image is minimized.

This kind of approach based on a multiple criterion candidate extraction, followed by a total variation regularized choice, can be applied to other problems. For instance, it can be useful for optical flow estimation [START_REF] Fortun | Aggregation of local parametric candidates with exemplar-based occlusion handling for optical flow[END_REF][START_REF] Fortun | A variational aggregation framework for patch-based optical flow estimation[END_REF] or for video colorization [START_REF] Pierre | Interactive video colorization within a variational framework[END_REF].

These problems are usually nonconvex and often nonsmooth. To solve them, authors of [START_REF] Pierre | Luminance-chrominance model for image colorization[END_REF] have proposed a primal-dual algorithm to make their problem numerically tractable but without convergence proof. Indeed, the variational model used in [START_REF] Pierre | Luminance-chrominance model for image colorization[END_REF] is biconvex, that is, the functional depends on two block-variables, with convexity in each block. This partial convexity feature allows to apply schemes encountered in convex optimization which alternate partial (block) optimizations, since each subproblem remains convex. However, for nonconvex functionals, no convergence is guaranteed.

Biconvex Optimization

Despite some notable progress in the field of continuous nonconvex optimization, the minimization of a biconvex functional under separable biconvex constraints min (x,y)∈C X ×C Y J(x, y) with C X , C Y convex, J biconvex remains a difficult task. The partial convexity carried by the biconvex structure of the functional yet allows applying convex optimization algorithm in a block-optimization framework, but convergence may be hard to prove. Moreover, unlike in the convex case, such schemes can usually be guaranteed to converge to a critical point, which may not be a (even local) minimizer. In 2013, [START_REF] Xu | A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion[END_REF] and [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF] proposed an alternating proximal scheme (known as the PALM method in [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF]) which aims at solving nonsmooth and nonconvex optimization problems, by alternating the well-known forward-backward splitting on each partial problem. They provided a convergence analysis based on the Kurdyka-Lojasiewicz 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], and the proposed scheme presents in practice similar numerical performance as when used for convex optimization. However, the application hypotheses of this algorithm only allow to consider simple regularizers4 , in the sense that their Moreau proximity operator [START_REF] Moreau | Proximité et dualité dans un espace hilbertien[END_REF] prox

F (x 0 ) := arg min x∈X F (x) + 1 2 x -x 0 2
must have a closed form for any x 0 ∈ X , or be exactly computable in a reasonable time.

In practice, this limits the class of functionals which can be optimized by this scheme.

In the case where the regularizers are not simple but sufficiently smooth (or replaced by a smooth approximation), they can be incorporated in the coupling term, but may lead to smaller stepsizes (hence, slower convergence). Thus, in 2017, the authors in [START_REF] Nikolova | Alternating proximal gradient descent for nonconvex nonsmooth block-regularised problems[END_REF] proposed another alternating scheme, which is basically a mirror of PALM, in the sense that the forward and the backward steps are not applied on the same part of the functional. This led to the ASAP algorithm, whose convergence has been studied in details in [START_REF] Nikolova | Alternating proximal gradient descent for nonconvex nonsmooth block-regularised problems[END_REF] for a large class of nonsmooth and nonconvex functional (among others biconvex functionals). Like the PALM method, the ASAP algorithm yet presents two main drawbacks: the simplicity hypothesis for the functional (to ensure that the backward step is computable), and the convergence rate, which is empirically in O(1/K) (with K the number of iterations).

A way to relax this assumption is to consider more general proximity operators, typically by using the so-called Bregman distances [START_REF] Bregman | The relaxation method of finding the common point of convex sets and its application to the solution of problems in convex programming[END_REF][START_REF] Censor | Proximal minimization algorithm withdfunctions[END_REF]. Indeed, some functions may not be simple with respect to the usual Moreau proximity operator, but simple to a more general one defined thanks to an adapted Bregman distance. We will show in this paper that such a generalization can be carried out in the framework introduced in [START_REF] Nikolova | Alternating proximal gradient descent for nonconvex nonsmooth block-regularised problems[END_REF]. Besides, the convergence speed issue has already been successfully tackled in the convex optimization by using inertia [START_REF] Boris | Some methods of speeding up the convergence of iteration methods[END_REF][START_REF] Güler | On the convergence of the proximal point algorithm for convex minimization[END_REF][START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF][START_REF] Chambolle | Accelerated alternating descent methods for Dykstra-like problems[END_REF][START_REF] Goldstein | Fast alternating direction optimization methods[END_REF] and for some nonconvex proximal schemes [START_REF] Ochs | ipiano: Inertial proximal algorithm for nonconvex optimization[END_REF][START_REF] Pock | Inertial proximal alternating linearized minimization (iPALM) for nonconvex and nonsmooth problems[END_REF]. We will show in this paper that ASAP and its Bregman-distance based generalization can also incorporate such a speed-up strategy, which leads to our proposed method. Numerical experiments will prove that an acceleration can be observed.

Apart from the colorization problem studied in this paper, many other applications can be concerned by the proposed accelerated scheme. Indeed, most of the joint estimation problems (i.e. the estimation of two or more objects during the same process) increasingly encountered in image processing result from the combination of the variational model of each variable, which may usually be chosen convex, so that the resulting functional is nonconvex but still biconvex (or multiconvex), see e.g. [START_REF] Barbanson | Relief computation from images of a fourier transform spectrometer for interferogram correction[END_REF][START_REF] Tan | Correction par méthode variationnelle des non uniformités des détecteurs d'un interféromètre imageur[END_REF][START_REF] Burger | A variational model for joint motion estimation and image reconstruction[END_REF].

The paper is organized as follows: in the next section, we present the colorization method proposed in [START_REF] Pierre | Luminance-chrominance model for image colorization[END_REF] and the associated variational model. In particular, we analyze its regularity. Then, in Section 3, we recall some mathematical preliminaries, mainly about subdifferential calculus and Bregman generalization of the forward-backward splitting. In Section 4, we introduce the general form of our proposed algorithm and discuss its convergence properties for large classes of nonconvex nonsmooth functions. The parameter choice is discussed in Section 5. Eventually, in Section 6, we present the numerical application of our method for the image colorization problem, and compared its performances with three other methods.

Variational Approach for Image Colorization

Variational Model

The model of [START_REF] Pierre | Luminance-chrominance model for image colorization[END_REF], implemented in [START_REF] Pierre | Colociel. Dépôt Agence de Protection des Programmes No IDDN[END_REF], settles on a multiple candidates selection strategy. The model is written in the Y U V color-space, and the target image (the image to colorize) is considered to be the luminance (Y ) channel of the expected color image. Since they carry the color information, the aim is to compute the U and V channels of this image from the reference (or source) image.

The method proposed in [START_REF] Pierre | Luminance-chrominance model for image colorization[END_REF] estimates the U and V channels in a two-step procedure. In the first step, for each pixel of the target image, C color candidates are selected from the source image via a candidate extraction method (see Section 6). Then, a voting process is run, to select one candidate among the C candidates extracted at the previous step. Once the selection is done for each pixel, the resulting U and V channels are combined with the luminance channel to get the expected color image.

In [START_REF] Pierre | Luminance-chrominance model for image colorization[END_REF], the voting process is done in a variational framework. For any pixel x ∈ Ω, given C chrominance candidates, denoted by c i (x), the choice of a specific candidate is modelized thanks to a vector w(x) ∈ R C belonging to the simplex Σ C which represents the distribution of the votes (or the likelihood) for each candidate. Thus, at the end of the selection process, one retains the candidate c i (x) which obtained most votes, namely such as w i (x) is the largest coefficient of w(x). In [START_REF] Pierre | Luminance-chrominance model for image colorization[END_REF], the vote vector (and thus, its largest value) is chosen so that the resulting color image is feasible (i.e. within a given range R :

= {u | u(x) = (U (x), V (x)) ∈ [U min , U max ] × [V min , V max ]}
) and has a small total variation. Thus, it has been proposed to minimize the following variational model:

J(u, w) = TV C (u) + χ R (u) + χ Σ (w) + λ 2 Ω C i=1 w i (x) u(x) -c i (x) 2 dx. ( 1
)
In this model, u = (U, V ) represents the chrominance channels of the sought-after color image and x is the pixel position in the image domain Ω. λ > 0 is a scalar parameter which controls the regularization of the results. The coupled total variation TV C (u) is defined as follows:

TV C (u) = Ω γ ∇Y 2 + ∇u 2 (2)
and can be seen as the smoothed by γ ∇Y 2 total variation of the color image in the Y U V space, with weight γ > 0 for the luminance channel. As shown is [START_REF] Pierre | Luminance-chrominance model for image colorization[END_REF], this term enforces the correlation between the (given) luminance variations and the (estimated) chrominances variations. Thus, the optimization of the model in (1) yields the joint estimation of both the weight (votes) w and the chrominance of the sought-after color image.

As underlined in the introduction, the main difficulty of the selection process is the optimization of (1). Indeed, it is obvious that the constraints are bounded, convex and separable, the TV C regularizer is convex, but the coupling term is solely biconvex. Thus, the whole functional J is biconvex, and in particular nonconvex.

Regularity of the Problem

Biconvex optimization can be handled by various algorithms. Among them, two methods are noteworthy, namely the PALM method and the ASAP algorithm. The application of these two methods requires simplicity and smoothness conditions of the functional J. The critical conditions in the studied problem [START_REF] Attouch | Convergence of descent methods for semialgebraic and tame problems: proximal algorithms, forward-backward splitting, and regularized Gauss-Seidel methods[END_REF] are the smoothness of the coupled total variation and the simplicity and/or smoothness of the coupling term.

A simple way to ensure the smoothness of the coupled total variation regularizer for any grayscale image Y is to introduce a threshold α > 0 for the weighted gradient of the luminance Y , namely to replace γ ∇Y 2 by max{α, γ ∇Y 2 } in (2), so that the resulting regularized function TV α C has a Lipschitz continuous gradient given by

∇TV α C (u) = - div(∇u) max {γ ∇Y 2 , α} + ∇u 2 with modulus L ∇F ≤ |||∇||| 2 / √ α.
In practice, given that, except in completely flat areas in the image (which is unlikely for real noisy images), ∇Y is expected to be numerically nonnull, if α is chosen small enough, one has the identity TV α C (u) = TV C (u) for any u. Otherwise said, for proper choices of α, this shows that for real images, the coupled total variation is smooth with a Lipschitz gradient, with a modulus bounded thanks to α.

The remaining terms -the convex constraints and the coupling term -can be easily handled with standard implicit or explicit numerical schemes. Thus, both ASAP and PALM are directly relevant for this problem. However, it is noteworthy that the computation of the proximity operator which involves a simplex constraint may be timeconsuming, but can also be cleverly handled by generalized Bregman proximity operators.

Mathematical Preliminaries

Notations

Let X and Y be real finite-dimensional spaces. The ith element of a vector or a matrix x (seen as a vector) reads as x i . For an m × n real matrix w we denote

w := w F = i,j w 2 i,j , (3)
noticing that if w is a vector (n = 1), the Frobenius norm • F boils down to the 2 norm. Given a nonempty set S ⊂ X , the distance of any point x + ∈ X to S is defined by dist(x + , S)

:= inf{ x -x + | x ∈ S}, χ S (x) := 0 if x ∈ S, +∞ if x ∈ S.

Subdifferential and Partial Subdifferentials

Let us first recall the definition of the subdifferential of a convex function.

Définition 1 Let f : R m → R ∪ {+∞} be a proper 5 , convex, lower semicontinuous (l.s.c.) function and x + ∈ dom f . The subdifferential ∂f (x + ) of f at x + is the set of all p ∈ R m , called subgradients of f at x + , such that

∀ x ∈ R m f (x) ≥ f (x + ) + p, x -x + If x + / ∈ dom h, then ∂f (x + ) = ∅.
Its main interest is to characterize the minimizers x * (when they exist) of f , via the well-known Fermat's rule 0 ∈ ∂f (x * ) (which is the generalization of the first-order optimality condition for smooth functions 0 = ∇f (x * )). To extend this rule to nonconvex nonsmooth function (Proposition 1), we introduce the following extensions of the subdifferential:

Définition 2 [41, Def. 8.3] Let f : R m → R ∪ {+∞} be a proper function and let

x + ∈ dom f . 1. The Fréchet subdifferential ∂f (x + ) of f at x + is the set of all p ∈ R m such that f (x) ≥ f (x + ) + p, x -x + + o( x -x + ). 2. The (limiting-)subdifferential ∂f (x + ) of f at x + is the set of all p ∈ R m such that there exist a sequence {x k } k∈N ∈ (R m ) N converging to x + and a sequence {p k } k∈N ∈ (R m ) N converging to p satisfying ∀ k ∈ N, p k ∈ ∂f (x k ) and f (x k ) → f (x + ). If x + / ∈ dom f , then ∂f (x + ) = ∂f (x + ) = ∅.
Remark 1 If f is convex, then the definitions of the subdifferential in Definitions 1 and 2 coincide, and ∂f (x) = ∂f (x) for any x ∈ dom f . If f is continuously differentiable at x, one has ∂f (x) = {∇f (x)}.

As for smooth / convex functions, the subdifferential gives a necessary first-order optimality condition, known as Fermat's rule. Given a proper function h : R m × R n → R ∪ {+∞} and (x + , y + ) ∈ dom h, we define its partial subdifferentials ∂ x h(x + , y + ) and ∂ y h(x + , y + ) at (x + , y + ) respectively as the subdifferentials of the partial functions x → h(x, y + ) at x + and y → h(x + , y) at y + . Unlike the differentiable case, the links between the partial subdifferentials and the subdifferential may be quite complicated. In particular, in general, one has

∂h(x, y) = ∂ x h(x, y) × ∂ y h(x, y).
As explored in details in [START_REF] Nikolova | Alternating proximal gradient descent for nonconvex nonsmooth block-regularised problems[END_REF], this may make alternating schemes such as the proposed method fail to find critical points. Another issue which may arise when dealing with such schemes in the nonsmooth and nonconvex case is the following one. The notion of limiting subdifferential is introduced because it brings more robustness than that of the Fréchet subdifferential, in the sense that it satisfies the following closedness property: if

{x k } k∈N in dom f converges to x * ∈ dom f such that f (x k ) → f (x * ) and p k → p * with p k ∈ ∂f (x k ) for any k ∈ N, then one has p * ∈ ∂f (x *
). Such a closedness property is usually necessary to preserve, at convergence, useful property for sequences generated by first-order optimization schemes. However, as shown further in the convergence analysis, the proposed scheme may need a stronger closedness property on the partial subdifferentials.

Proposition 2 Let h : R m × R n → R ∪ {+∞} be a proper, biconvex function continuous on its domain. Then both its x-and y-partial subdifferential are parametrically closed at any point (x * , y * ) ∈ dom h, that is, for any sequence {(x k , y k )} k∈N which converges to (x * , y * ) such that h(x k , y k ) → h(x * , y * ) and for any p k

x → p x (resp. p k y → p y ) with p k x ∈ ∂ x h(x k , y k ) (resp. p k y ∈ ∂ y h(x k , y k )) for any k ∈ N, one has p x ∈ ∂ x h(x * , y * ) (resp. p y ∈ ∂ y h(x * , y * )). Proof . Let {(x k , y k )} k∈N be a sequence as in Proposition 2. Let k ∈ N. Since x → h(x, y k ) is a convex, l.s.c.

and proper functions, the subgradient inequality (Definition 1) leads to

∀ x ∈ R m , h(x, y k ) ≥ h(x k , y k ) + p k x , x -x k , If we let k → +∞ in the inequality above, we get by continuity ∀ x ∈ R m , h(x, y * ) ≥ h(x * , y * ) + p x , x -x *
which means that p x is a subgradient of the convex function x → h(x, y * ) at x * . Similar computations for the convex function y → h(x, y) completes the proof.

Note that, unlike the closedness of the subdifferential, the parametric closedness of the partial subdifferentials may not hold for general nonconvex functions, even for quite smooth functions (consider for instance h(x, y) = √ xy + χ R + (x) + χ R + (y) at (0, 0)).

Bregman Distances

Proximal methods have brought notable improvement in the field of continuous optimization, but the main issue is their applicability for general problems, due to the noncomputability of many proximity operators. To deal with this problem, we consider in this paper a generalization of the Moreau proximity operator, which aims at making this operator computable in many cases, as shown in Section 6.

Définition 3 (Bregman distance 6 [START_REF] Bregman | The relaxation method of finding the common point of convex sets and its application to the solution of problems in convex programming[END_REF]) Let b : R m → R ∪ {+∞} be a strictly convex function and differentiable on int(dom b). The Bregman "distance" associated to b is defined for any

(x, y) ∈ R m × R m by D b (x, y) := b(y) -b(x) -∇b(x), y -x if (x, y) ∈ int(dom b) × dom b +∞ otherwise.
Σ m := x = (x i ) ∈ R m | m i=1 x i = 1 and x i ≥ 0 .
We also consider the entropy function, given by b e (x) :=

     m i=1 x i log x i if x ∈ Σ m +∞ otherwise,
where 0 × log 0 = 0 by convention. Then, the Bregman distance associated to the entropy function is explicitly given for (x, y)

∈ (intΣ m ) × Σ m by D be (x, y) = m i=1 y i (log y i -log x i ) .
Proposition 3 Let b : R m → R ∪ {+∞} be a strongly convex function of modulus L, differentiable on int(dom b). Then for any (x, y) ∈ R m × R m , one has the following lower bound:

D b (x, y) ≥ L 2 x -y 2 . ( 4
)
Proof . The proof is a direct consequence of the strong convexity.

The following lemma will be further used in proofs.

Lemma 1 ([4]

) Let b : R m → R∪{+∞} be a strictly convex function and differentiable on int(dom b). Then, one has the following identity for any three points x, y ∈ int(dom b) and z ∈ dom b:

D b (x, z) + D b (y, x) -D b (y, z) = z -x, ∇b(y) -∇b(x) .

Bregman Proximal Gradient Descent

In this section, we will introduce a generalization of well-known Forward-Backward Splitting (FBS) based on a generalization of the Moreau proximity operator that uses the Bregman distances introduced above. This generalization was initially introduced in [START_REF] Censor | Proximal minimization algorithm withdfunctions[END_REF].

Définition 4 Let b : R m → R ∪ {+∞} be a strictly convex function and differentiable on int(dom b). Let h : R m → R ∪ {+∞} be a proper function and let f : R m → R be a differentiable function. Then we define for any

(u, v) ∈ int(dom b) × R m P b h (u; ∇f (v)) := arg min x∈R m h(x) + x -u, ∇f (v) + D b (u, x) . (5)
When b = • is the Frobenius norm, the optimization problem (5) when u = v is equivalent to the FBS or proximal gradient descent applied to h + f . Remark 4 If h is a proper convex function and if b is strongly convex, P b h (u; ∇f (v)) is nonempty. Indeed, the function involved in ( 5) is convex and lowerbounded by a strongly convex function. Hence, it is coercive and thus admits at least a minimizer on its domain.

The following proposition is a descent-like property for the Bregman FBS.

Proposition 4 Let b : R m → R∪{+∞} be a strictly convex function and differentiable on int(dom b). Let h : R m → R ∪ {+∞} be a proper convex function and let f : R m → R be a continuously differentiable function, with L ∇f -Lipschitz continuous gradient. Then, for any u ∈ int(dom b) and any τ > 0, if x + ∈ P b τ h (u; τ ∇f (v)), one has for any

x ∈ R m h(x + ) + f (x + ) + 1 τ D b (u, x + ) + D b (x + , x) ≤ h(x) + f (x) + 1 τ D b (u, x) + L ∇f + s 2 x -x + 2 + L 2 ∇f 2s x -v 2 .
Before giving the proof of Proposition 4, let us recall two useful lemmas.

Lemma 2 (Descent lemma [START_REF] Bertsekas | Nonlinear programming[END_REF]) Let f : R m → R be a differentiable function with L ∇f -Lipschitz continuous gradient. Then, for any

(x, u) ∈ R m × R m , one has f (x) ≥ f (u) -u -x, ∇f (x) - L ∇f 2 x -u 2 . ( 6
)
Lemma 3 (Young inequality) For any s > 0,

x, y ≤ s 2 x 2 + 1 2s y 2 .
Proof of Proposition 4. Fermat's rule for (5) yields

-∇f (v) - ∇b(x + ) -∇b(u) τ ∈ ∂h(x + ).
By definition of the subgradient of a convex function, one has for any

x ∈ R m h(x + ) -x -x + , ∇f (v) + ∇b(x + ) -∇b(u) τ ≤ h(x).
Applying Lemma 1 with (x, y, z) = (x + , u, x), we get

h(x + ) -x -x + , ∇f (v) + 1 τ D b (x + , x) + D b (u, x + ) -D b (u, x) ≤ h(x).
Adding ( 6) with u = x + , using Lemma 3 with some s > 0 to bound the scalar product and using the Lipschitz continuity of ∇f then yields the desired result.

4 Proposed Algorithm

Optimization Problem

We consider the following two-block unconstrained optimization problem min

(x,y)∈X ×Y J(x, y) := F (x) + G(y) + H(x, y). (7)
We assume that the objective function J satisfies some of the following assumptions:

Assumptions (H1) (a) J : X × Y → R ∪ {+∞}

Assumptions (H2)

(a) H is continuous on its closed domain;

(b) For any (x, y) ∈ dom H, the partial subdifferentials satisfy ∂ x H(x, y)×∂ y H(x, y) ⊂ ∂H(x, y);

(c) H is biconvex, i.e.
x → H(x, y) is convex for any y ∈ Y and y → H(x, y) is convex for any x ∈ X .

Remark 5 From Assumptions (H1) and (H2), J is continuous on its closed domain dom J = dom H, which means that for any convergent sequence {(x k , y k )} k∈N in dom J with limit (x * , y * ) ∈ dom J, one has J(x k , y k ) → J(x * , y * ).

Example 2 Assumption (H2)(b) is satisfied as soon as H is differentiable, or additively separable, or a sum of such functions. Otherwise said, Assumption (H2)(b) holds true if

H(x, y) = h(x, y) + f (x) + g(y),
with h differentiable, f and g proper l.s.c. functions.

Assumptions (H3) H : X × Y → R ∪ {+∞} can be split into H(x, y) = h(x, y) + f (x) + g(y)
, where (a) f : X → R ∪ {+∞} and g : Y → R ∪ {+∞} are continuous on their domain;

(b) h : X × Y → R is continuous;

(c) x → h(x, y) is continuously differentiable for any y ∈ dom g. Moreover, y → ∇ x h(x, y) locally Lipschitz continuous on dom f , in the sense that for each bounded subset B X ×B Y ⊂ dom f ×dom g, there is a constant ξ > 0 such that for any x ∈ B X and (y, y ) ∈ B 2 Y , one has

∇ x h(x, y) -∇ x h(x, y ) ≤ ξ y -y .
Obviously, (H3)(c) holds for Lipschitz continuously differentiable functions h. In this case, (H2) holds as soon as (H3) is satisfied (according to Example 2).

Assumption (H4)

The function J satisfies the Kurdyka-Lojasiewicz (K L) property on every point (x + , y + ) of its domain, namely there exist η ∈ (0, +∞], a neighborhood O(x + , y + ) of (x + , y + ) and κ > 0 such that for any (x, y) ∈ O(x + , y + ) ∩ N η (J(x + , y + ))

κ dist(0, ∂J(x, y)) ≥ J(x, y) -J(x + , y + ) θ , (8)
where θ ∈ [0, 1) and N η (t) := J -1 (t, t + η) .

Remark 6 Assumption (H4) is crucial to prove strong convergence for nonconvex and nonsmooth optimization schemes such as [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF]. The main apparent difficulty is to verify that the functional J fulfills [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF]. In practice, it is sufficient to show that they are subanalytic [9, Theorem As shown in the analysis below, adding hypotheses on the optimization problem (7) leads to stronger convergence results for the proposed method. For the sake of simplicity, we will mainly focus on the case where H is biconvex 7 . However, the hypotheses (H2) and (H3) given in this paper can be relaxed while yielding same convergence results. The interested reader may see [START_REF] Nikolova | Alternating proximal gradient descent for nonconvex nonsmooth block-regularised problems[END_REF] for more details.

Inertial Bregman Alternating Structure-Adapted Proximal Gradient Descent

We can now introduce the general form of the proposed algorithm in Algorithm 1.

It basically consists in alternating partial (block) minimization steps, each step being decomposed into two overrelaxation (inertial) steps (in the spirit of the strategy proposed by Nesterov for projected gradient descent), followed by a Bregman FBS. When b X and b Y are Frobenius norm,

α k 1 = α k 2 = β k 1 = β k 2 = 0 (no relaxation) and (τ k , σ k ) = (τ, σ) (constant
stepsizes), one recovers the ASAP (Alternating Structure-Adapted Proximal gradient descent) algorithm introduced in [START_REF] Nikolova | Alternating proximal gradient descent for nonconvex nonsmooth block-regularised problems[END_REF]. Thus, the proposed scheme can be seen as an inertial Bregman generalization of ASAP with varying stepsizes.

Remark 7 Let us make some preliminary remarks about Algorithm 1. for all k ≥ 0, do

xk = x k + α k 1 (x k -x k-1 ), (9) xk
= x k + α k 2 (x k -x k-1 ), (10)
x k+1 ∈ P b X τ k H(•,y k ) (x k ; τ k ∇F (x k )); ( 11
) ȳk = y k + β k 1 (y k -y k-1 ), ( 12
) ŷk = y k + β k 2 (y k -y k-1 ), ( 13
)
y k+1 ∈ P b Y σ k H(x k+1 ,•) (ŷ k ; σ k ∇G(ȳ k )). ( 14
)
end for (b) If b X and b Y are strongly convex and H lower bounded, then ( 11) and ( 14) admit at most one solution. The solution is uniquely defined iff the function minimized in the cited problems are not identically equal to +∞, that is, iff there exists (x, y) ∈ X × Y such that (x, y k ), (x k+1 , y) ∈ dom H.

In order to ensure that the iterations in Algorithm 1 are well-defined and to prove convergence results, we need to make some assumptions on the algorithm parameters, i.e. on the choice of the Bregman distances (b Note that all Assumptions (A1) hold when b X and b Y are defined thanks to a norm

X , b Y ) (Assumptions (A1)-(A2)), the inertial parameters {(α k 1 , α k 2 , β k 1 , β k 2 )} k∈N (Assumptions (A2)-( A3
• M . Assumptions (A2) (a) If dom b X = X (resp. dom b Y = Y), then α k 2 = 0 (resp. β k 2 = 0) for any k ∈ N. (b) If α k 2 > 0 (resp. β k 2 > 0), then ∇b X (resp. ∇b Y ) is Lipschitz continuous of modulus L ∇b X (resp. L ∇b Y ).
Assumption (A2)(a) guarantees that the overrelaxed sequences {x k } k∈N and {ŷ k } k∈N stay in int(dom b X ) and int(dom b Y ) respectively, so that the iterates {(x k , y k )} k∈N are well-defined (see Remark 4).

Assumptions (A3)

Let {s k } k∈N and {t k } k∈N be two sequences of positive numbers. Set for any k ≥ 1

B k X := L 2 ∇F (α k 1 ) 2 2 s k + L 2 ∇b X α k 2 2 τ k ≥ 0 and B k Y := L 2 ∇G (β k 1 ) 2 2 t k + L 2 ∇b Y β k 2 2 σ k ≥ 0
(if ∇b X is not Lipschitz continuous and α k 2 = 0, then we set L 2 ∇b X α k 2 = 0 by convention; same for L 2 ∇b Y β k 2 ) and

A k X := 2 L X -α k-1 2 2 τ k - L ∇F + s k-1 2 and A k Y := 2 L Y -β k-1 2 2 σ k - L ∇G + t k-1 2 . Assume that τ k , σ k , s k , t k > 0 and α k 1 , α k 2 , β k 1 , β k 2 ≥ 0 are chosen such that ρ X := inf k A k X -B k X > 0, ρ Y := inf k A k Y -B k Y > 0 and let ρ := min{ρ X , ρ Y }. Remark 8 For any α k 1 , α k 2 , β k 1 , β k 2 ≥ 0 and s k , t k > 0, the quantities B k X and B k Y are nonnegative. Moreover, they cancel iff α k 1 = α k 2 = 0 and β k 1 = β k 2 = 0 respectively.

Hence, when Assumptions (A3) hold, A k

X and A k Y are necessarily positive and upper bounded.

Eventually, we introduce the following assumption:

Assumption (A4) Assume that τ k , σ k , s k , t k > 0 and α k 1 , α k 2 , β k 1 , β k 2 ≥ 0 are chosen such that inf k A k X > sup k B k X and inf k A k Y > sup k B k Y with A k X , A k Y , B k X and B k Y defined in Assumption (A4).
It is obvious that Assumption (A4) is a particular case of Assumptions (A3). When all the sequences involved in Assumption (A4) are chosen constant, (A3) and (A4) are equivalent.

Convergence Analysis

Proposition 5 (Objective convergence) Let Assumptions (H1), (A1)-(A3) hold. Let {(x k , y k )} k∈N be a sequence generated by Algorithm 1. Then the following assertions hold:

(a) For every k ≥ 1, one has

J(x k+1 , y k+1 ) + A k X x k+1 -x k 2 + A k Y y k+1 -y k 2 ≤ J(x k , y k ) + B k X x k -x k-1 2 + B k Y y k -y k-1 2 (b) lim k→∞ x k-1 -x k = 0 and lim k→∞ y k+1 -y k = 0. (c) The sequence {J(x k , y k )} k∈N is convergent. Proof . (a) Let k ≥ 1. Applying Proposition 4 with b = b X , h = H(•, y k ), f = F , x + = x k+1 , x = x k , u = xk , v = xk , and s = s k > 0 yields J(x k+1 , y k ) + 1 τ k D b X (x k , x k+1 ) + D b X (x k+1 , x k ) ≤ J(x k , y k ) + 1 τ k D b X (x k , x k ) + L ∇F + s k 2 x k -x k+1 2 + L 2 ∇F 2s k x k -xk 2
By definition of xk and xk , one has

x k -xk = α k 1 (x k -x k-1
). Hence, the inequality above also reads

J(x k+1 , y k ) + 1 τ k D b X (x k , x k+1 ) + D b X (x k+1 , x k ) - L ∇F + s k 2 x k -x k+1 2 (15) ≤ J(x k , y k ) + 1 τ k D b X (x k , x k ) + L 2 ∇F (α k 1 ) 2 2s k x k -x k-1 2 .
Using Proposition 1 with x = x k , y = xk , and z = x k+1 , we then get

J(x k+1 , y k ) + 1 τ k D b X (x k , x k+1 ) + D b X (x k+1 , x k ) - L ∇F + s k 2 x k -x k+1 2 ≤ J(x k , y k ) + L 2 ∇F (α k 1 ) 2 2s k x k -x k-1 2 + 1 τ k ∇b(x k ) -∇b(x k ), x k+1 -x k .
First assume that α k 2 > 0. Then, Lemma 3 along with the Lipschitz continuity of ∇b X yield

J(x k+1 , y k ) + 1 τ k D b X (x k , x k+1 ) + D b X (x k+1 , x k ) - L ∇F + s k 2 x k -x k+1 2 ≤ J(x k , y k ) + L 2 ∇F (α k 1 ) 2 2s k x k -x k-1 2 + 1 τ k L 2 ∇b X 2α k 2 xk -x k 2 + α k 2 2 x k+1 -x k 2 .
Then, noticing that xk -

x k = α k 2 (x k-1 -x k
) and rearranging the terms in the inequality above lead to

J(x k+1 , y k ) + 1 τ k D b X (x k , x k+1 ) + D b X (x k+1 , x k ) - L ∇F + s k 2 + α k 2 2τ k x k+1 -x k 2 ≤ J(x k , y k ) + L 2 ∇F (α k 1 ) 2 2s k + L 2 ∇b X α k 2 2τ k x k -x k-1 2 . ( 16
)
If α k 2 = 0, ( 15) and ( 16) are the same, so ( 16) holds whenever α k 2 ≥ 0. Now, since b X is L X -convex, one can lowerbound the previous inequality using (4), which gives:

J(x k+1 , y k ) + 2 L X -α k 2 2τ k - L ∇F + s k 2 x k+1 -x k 2 ≤ J(x k , y k ) + L 2 ∇F (α k 1 ) 2 2s k + L 2 ∇b X α k 2 2τ k x k -x k-1 2 .
Using the notations introduced in Assumptions (A2), one finally gets

J(x k+1 , y k ) + A k X x k+1 -x k 2 ≤ J(x k , y k ) + B k X x k -x k-1 2 . ( 17
)
Similar computations yield

J(x k+1 , y k+1 ) + A k Y y k+1 -y k 2 ≤ J(x k+1 , y k ) + B k X y k -y k-1 2 . (18)
Summing ( 17) and ( 18) then yields (a).

(b) Since A k X > B K X ≥ 0 and A k Y > B K X ≥ 0, summing the inequality (a) for k = 0 to k = K -1 yields J(x K , y K ) + K-1 k=0 A k X x k+1 -x k 2 + A k Y y k+1 -y k 2 ≤ J(x 0 , y 0 ).
Using that J(x K , y K ) is lower bounded by inf J ∈ R and that A k X , A k Y ≥ ρ > 0 yields that the series x k+1 -x k 2 and y k+1 -y k 2 converge. Thus, the sequences

{ x k+1 -x k } k∈N and { y k+1 -y k } k∈N converge to zero. (c) Since A k X > B K X ≥ 0 and A k Y > B K Y ≥ 0, the inequality (a) implies that the sequence {J(x k , y k ) + B k X x k -x k-1 2 + B k Y y k -y k-1 2
} k∈N is nonincreasing. Moreover, it is lower bounded, as J is. Thus it converges. According to the previous point, { x k+1x k } k∈N and { y k+1 -y k } k∈N converge to zero. Thus, {J(x k , y k )} k∈N converges as well.

Remark 9 Unlike the ASAP method, Algorithm 1 may generate objective-convergent sequences which are not monotone.

The set of all limit points of a sequence {(x k , y k )} k∈N generated by Algorithm 1 starting from a point (x 0 , y 0 ) is denoted by L(x 0 , y 0 ). Proposition 6 (Subsequential convergence to critical points) Let Assumptions (H1), (H2) and (A1)-(A3) hold and let {(x k , y k )} k∈N be a sequence generated by Algorithm 1 which is assumed to be bounded. Let (x * , y * ) ∈ L(x 0 , y 0 ).

(a) there is a subsequence (x kj , y kj ) j∈N converging to (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. (d) lim k→∞ dist((x k , y k ), critJ) = 0.
Proof . (a) is a consequence of the boundedness assumption. (b) Since J is continuous on its closed domain, and (x k , y k ) ∈ dom H = dom J (since there are defined as minimizers), one has that (x kj , y kj ) → (x * , y * ) ∈ dom J and lim j→∞ J(x kj , y kj ) = J(x * , y * ).

(c) Using Fermat's rule on the definition of x k+1 and y k+1 yields after rearrangement

q k+1 x := ∇F (x k+1 ) -∇F (x k ) - ∇b X (x k+1 ) -∇b X (x k ) τ k ∈ ∇F (x k+1 ) + ∂ x H(x k+1 , y k ) = ∂ x J(x k+1 , y k ) p k+1 y := ∇G(y k+1 ) -∇G(ȳ k ) - ∇b Y (y k+1 ) -∇b Y (ŷ k ) σ k ∈ ∇G(y k+1 ) + ∂ y H(x k+1 , y k+1 ) = ∂ y J(x k+1 , y k+1 ).
Note that, according to Proposition 5(b), one has

x k+1 -xk = x k+1 -x k -α k 1 (x k -x k-1 ) → 0 x k+1 -xk = x k+1 -x k -α k 2 (x k -x k-1 ) → 0 y k+1 -ȳk = y k+1 -y k -β k 1 (y k -y k-1 ) → 0 y k+1 -ŷk = y k+1 -y k -β k 2 (y k -y k-1 ) → 0
Then, the Lipschitz continuity of ∇F , ∇G, ∇b X , and ∇b Y implies that q k x → 0 and p k y → 0. One eventually gets the desired result by applying Proposition 2 with the subsequences {q kj

x } j and {p kj y } j , and using that, by hypothesis, ∂ x J(x * , y * )×∂ y J(x * , y * ) ⊂ ∂J(x * , y * ). (d) Suppose that {dist((x k , y k ), critJ)} k∈N does not go to zero as k → ∞. Then there exist M > 0 and (k j ) j∈N such that for any j ∈ N, dist (x kj , y kj ), critJ > M . However, since {(x kj , y kj )} j∈N is a subsequence of the bounded sequence {(x k , y k )} k∈N , it has convergent subsequence {(x kj n , y kj n )} n of limit (x * , y * ) ∈ critJ (according to (c)). Thus,

M < dist (x kj n , y kj n ), critJ ≤ (x kj n , y kj n ) -(x * , y * ) -→ n→∞ 0
which leads to a contradiction.

Remark 10 If dom J is bounded, then all the sequences generated by Algorithm 1 are bounded.

To prove the strong convergence of the iterates generated by Algorithm 1, we need to make additional assumptions and to prove some preliminaries lemmas. First introduce the following notations for ε x , ε y ≥ 0 and

X = (x 1 , x 2 ) ∈ U 2 and Y = (y 1 , y 2 ) ∈ V 2 : Ψ (εx,εy) (X, Y ) := J(x 1 , y 1 ) + ε x x 1 -x 2 2 + ε y y 1 -y 2 2 .
Assume that Assumption (A4) holds. Then, one can assume that

inf k A k X > ε x > sup k B k X and inf k A k Y > ε y > inf k B k Y ( 19
)
Corollaire 1 Let Assumptions (H1) and (A1)-(A4) hold. Let {(x k , y k } k∈N be a sequence generated by Algorithm 1. For any k ≥ 0, set X k = (x k , x k-1 ) and Y k = (y k , y k-1 ). Then, there exists a > 0 such that

Ψ (εx,εy) (X k+1 , Y k+1 ) + a (X k+1 -X k , Y k+1 -Y k ) 2 ≤ Ψ (εx,εy) (X k , Y k )
where (ε x , ε y ) satisfies [START_REF] Denkowska | A long and winding road to definable sets[END_REF].

Proof . According to Proposition 5(a), one has

J(x k+1 , y k+1 ) + ε x x k+1 -x k 2 + ε y y k+1 -y k 2 + (A k X -ε x ) x k+1 -x k 2 +(ε x -B k X ) x k -x k-1 2 + (A k Y -ε y ) y k+1 -y k 2 + (ε y -B k X ) y k -y k-1 2 ≤ J(x k , y k ) + ε x x k -x k-1 2 + ε y y k -y k-1 2 Thus, one can choose a := inf k {A k X -ε x , ε x -B k X , A k Y -ε y , ε y -B k Y }
, which is positive according to [START_REF] Denkowska | A long and winding road to definable sets[END_REF].

Corollaire 2 Let Assumptions (H1)-(H3) and (A1)-(A3) hold. Let {(x k , y k )} k∈N be a sequence generated by Algorithm 1 which is assumed to be bounded. Then for any k ≥ 0 one has (p k+1 x , p k+1 y ) ∈ ∂J(x k+1 , y k+1 ) and ξ > 0 such that

p k+1 x ≤ L ∇F + L ∇b X τ k x k+1 -x k + L ∇F α k 1 + L ∇b X τ k α k 2 x k -x k-1 + ξ y k -y k+1 p k+1 y ≤ L ∇G + L ∇b Y σ k y k+1 -y k + L ∇G β k 1 + L ∇b Y σ k β k 2 y k -y k-1 .
Proof . Since J(x, y k ) = F (x)+h(x, y k )+f (x)+constant, with x → h(x, y k ) continuously differentiable, simple subdifferential calculus shows that

∂ x J(x k+1 , y k ) = ∇F (x k+1 ) + ∇ x h(x k+1 , y k ) + ∂f (x k+1 ) = ∇F (x k+1 ) + ∇ x h(x k+1 , y k+1 ) + ∂f (x k+1 ) +∇ x h(x k+1 , y k ) -∇ x h(x k+1 , y k+1 ).
Hence, using the notations introduced in the proof of Proposition 6(c), one has (p k+1

x , p k+1 y ) ∈ ∂J(x k+1 , y k+1 ) with

p k+1 x := q k+1 x -∇ x h(x k+1 , y k ) + ∇ x h(x k+1 , y k+1 )
which satisfies the following bounds

p k+1 x ≤ ∇F (x k+1 ) -∇F (x k ) + 1 τ k ∇b X (x k+1 ) -∇b X (x k ) + ∇ x h(x k+1 , y k ) -∇ x h(x k+1 , y k+1 ) p k+1 y ≤ ∇G(y k+1 ) -∇G(ȳ k ) + 1 σ k ∇b Y (y k+1 ) -∇b Y (ŷ k ) .
Using the Lipschitz continuity of ∇F , ∇G, ∇b X , ∇b Y , and the local Lipschitz continuity of ∇ x h ({(x k , y k )} k∈N being bounded by assumption), one shows that there exists ξ > 0 such that

p k+1 x ≤ L ∇F x k+1 -xk + L ∇b X τ k x k+1 -xk +ξ y k -y k+1 p k+1 y ≤ L ∇G y k+1 -ȳk + L ∇b Y σ k y k+1 -ŷk
and one completes the proof using the definition of the relaxed sequences.

The following corollary shows that Algorithm 1 generates a sequence of subgradient which can be bounded.

Corollaire 3 Let Assumptions (H1)-(H3) and (A1)-(A3) hold. Let {(x k , y k } k∈N be a sequence generated by Algorithm 1 which is assumed to be bounded. For any k ≥ 0, set X k = (x k , x k-1 ) and Y k = (y k , y k-1 ). Then there exists ξ > 0 such that for any k ≥ 0 one has P k+1 ∈ ∂Ψ εx,εy (X k+1 , Y k+1 ) such that

P k+1 ≤ ξ (X k+1 -X k , Y k+1 -Y k )
where (ε x , ε y ) satisfies [START_REF] Denkowska | A long and winding road to definable sets[END_REF].

Proof . First note that, for any

X = (x 1 , x 2 ) ∈ X 2 and Y = (y 1 , y 2 ) ∈ Y 2 , ∂ X Ψ εx,εy (X, Y ) is the set of all P X = (p X,1 , p X,2 ) ∈ X 2 such that p X,1 ∈ ∂ x J(x 1 , y 1 ) + 2 ε x (x 1 -x 2 ) p X,2 = -2 ε x (x 1 -x 2 ) while ∂ Y Ψ εx,εy (X, Y ) is the set of all P Y = (p Y,1 , p Y,2 ) ∈ Y 2 such that p Y,1 ∈ ∂ y J(x 1 , y 1 ) + 2 ε y (y 1 -y 2 ) p Y,2 = -2 ε y (y 1 -y 2 ) Moreover, Ψ εx,εy (X, Y ) can be split into Ψ εx,εy (X, Y ) = S(X, Y ) + N (X, Y ) with S(X, Y ) = F (x 1 ) + G(y 1 ) + h(x 1 , y 1 ) + ε x x 1 -x 2 2 + ε y y 1 -y 2 2
a continuously differentiable function and N (X, Y ) = f (x 1 ) + g(y 1 ) which is separable. Thus, one can check that

∂ X Ψ εx,εy (X, Y ) × ∂ Y Ψ εx,εy (X, Y ) = ∂Ψ εx,εy (X, Y )
Hence, one has P k+1 ∈ ∂Ψ εx,εy (X k+1 , Y k+1 ) when setting

P k+1 := p k+1 x , -2 ε x (x k+1 -x k ), p k+1 y , -2 ε y (y k+1 -y k ) .
Let us bound P k+1 . Pythagora's theorem yields

P k+1 2 = p k+1 x 2 + 4 ε 2 x x k+1 -x k 2 + p k+1 y 2 + 4 ε 2 y y k+1 -y k 2 .
Using that (a + b) 2 ≤ a 2 /2 + b 2 /2 (Lemma 3 with s = 2), together with the bounds in Corollary 2 yields

P k+1 2 ≤ 1 4 L ∇F + L ∇b X τ k 2 + 4 ε 2 x x k+1 -x k 2 + 1 4 L ∇F α k 1 + L ∇b X τ k α k 2 2 x k -x k-1 2 + 1 2 L ∇G + L ∇b Y σ k 2 + 1 2 ξ 2 + 4 ε 2 y y k+1 -y k 2 + 1 2 L ∇G β k 1 + L ∇b Y σ k β k 2 2 y k -y k-1 2 .
Choosing

ξ 2 = max 1 4 L ∇F + L ∇b X τ k 2 + 4 ε 2 x , 1 4 L ∇F α k 1 + L ∇b X τ k α k 2 2 , 1 2 L ∇G + L ∇b Y σ k 2 + 1 2 ξ 2 + 4 ε 2 y , 1 2 L ∇G β k 1 + L ∇b Y σ k β k 2 2
completes the proof.

Proposition 7 (Strong convergence) Let Assumptions (H1)-(H4) and (A1)-(A4) hold. Let {(x k , y k } k∈N be a sequence generated by Algorithm 1 which is assumed to be bounded. Then {(x k , y k } k∈N is a convergent Cauchy sequence.

Proof . For any k ≥ 0, set X k = (x k , x k-1 ) and Y k = (y k , y k-1 ). Then the sequence {Ψ εx,εy (X k , Y k )} k∈N satisfies the three conditions needed in [1, Theorem 2.9], namely (a) sufficient decrease condition (Corollary 1);

(b) relative error condition (Corollary 3);

(c) continuity condition ({(X k , Y k )} k∈N is bounded).
Eventually, one may check that Ψ εx,εy is K L whenever the objective function J is.

Inertial Parameter Choice

In this section, we detail some feasible choices for the algorithm parameters so that Assumptions (A1)-(A4) are satisfied. We focus here on the constant case and also give some trick to allow large inertial parameters along with large stepsizes.

Constant Case

Consider the case when all the parameters are constant. Then, set

B 0 X := L 2 ∇F (α 1 ) 2 2 s + L 2 ∇b X α 2 2 τ ≥ 0 and B 0 Y := L 2 ∇G (β 1 ) 2 2 t + L 2 ∇b Y β 2 2 σ ≥ 0
(we recall that if ∇b X is not Lipschitz continuous and α 2 = 0, then L 2 ∇b X α 2 = 0 by convention; same for L 2 ∇b Y β 2 ) and

A 0 X := 2 L X -α 2 2 τ - L ∇F + s 2 and A 0 Y := 2 L Y -β 2 2 σ - L ∇G + t 2
In this case, Assumption (A3) reduces to find τ, σ, α 1 , α 2 , β 1 , β 2 > 0 and s, t > 0 such that [START_REF] Pock | Inertial proximal alternating linearized minimization (iPALM) for nonconvex and nonsmooth problems[END_REF], we may first choose (s, t) and the best choice is that which maximizes the upperbound for J(x k+1 , y k+1 ) -J(x k , y k ), namely that which maximizes both

ρ X := A 0 X -B 0 X > 0 and ρ Y := A 0 Y -B 0 Y > 0 According to
L ∇F + s 2 x k+1 -x k 2 + L 2 ∇F (α 1 ) 2 2 s x k -x k-1 2 and L ∇G + t 2 y k+1 -y k 2 + L 2 ∇G (β 1 ) 2 2 t y k -y k-1 2
Thus, the best choice for (s, t) is given by

L ∇F α 1 x k -x k-1 x k+1 -x k , L ∇G β 1 y k -y k-1 y k+1 -y k
According to [START_REF] Pock | Inertial proximal alternating linearized minimization (iPALM) for nonconvex and nonsmooth problems[END_REF], we may choose

s = L ∇F α 1 and t = L ∇G β 1
Thus, one has to find

α 1 , α 2 , β 1 , β 2 > 0 such that 2 -α 2 τ -L ∇F (1 + α 1 ) > α 2 τ + L ∇F α 1 and 2 -β 2 σ -L ∇G (1 + β 1 ) > β 2 σ + L ∇G β 1 i.e. 2 τ > L ∇F + 2 α 2 τ + L ∇F α 1 and 2 σ > L ∇G + 2 β 2 σ + L ∇G β 1
Note that, since the left-hand side terms are nonnegative, this implies that τ and σ should be chosen so that τ < 2/L ∇F and σ < 2/L ∇G . One has

ρ k x = 2 τ -L ∇F -2 α 2 τ + L ∇F α 1 and ρ k y = 2 σ -L ∇G -2 β 2 σ + L ∇G β 1
Let us consider two simple cases, namely the case where α 1 = α 2 =: α (resp. β 1 = β 2 =: β) and the the case where α 1 = 0 and α 2 =: α (resp. β 1 = 0 and β 2 =: β). In the first case, the feasibility conditions are

ρ X = 2 τ -L ∇F -2 α 1 τ + L ∇F > 0 and ρ Y = 2 σ -L ∇G -2 β 1 σ + L ∇G > 0
which means that, for given stepsizes,

1 ≥ 2/τ -L ∇F 2/τ + 2 L ∇F > α ≥ 0 and 1 ≥ 2/σ -L ∇G 2/σ + 2 L ∇G > β ≥ 0
Given the idea that the inertial parameters should be as large as possible (within a reasonable range) to have a visible effect on the algorithm convergence, one can choose for instance

α * (τ ) = 0.999 × 2/τ -L ∇F 2/τ + 2 L ∇F and β * (σ) = 0.999 × 2/σ -L ∇G 2/σ + 2 L ∇G
Then, we see that, the larger the stepsizes, the smaller the inertial parameters, since

α * (τ ) ≈ 0.999 × 2/τ -L ∇F 3 L ∇F and β * (σ) ≈ 0.999 × 2/σ -L ∇G 3 L ∇G
In the second case, the feasibility conditions become

1 - τ L ∇F 2 > α and 1 - σL ∇G 2 > β
One can see that, once again, these upperbounds decrease to zero as the stepsizes increase to their maximum values.

Adaptive Choice for Larger Inertial Parameters

For some problems, the feasible sets for inertial parameters leads to either small stepsizes (τ, σ) or small inertial parameters (α, β) (in the constant parameter cases). In both cases, the convergence is ensured by the results in the previous section, but the overrelaxation does not enhance the observed convergence rate. However, it may happen that larger inertial parameters than those allowed by the feasibility set lead in practice to empirical convergence. In particular, some none-null large inertial parameters (e.g. α = 1) can be chosen along with large stepsizes (i.e. close to their upperbound), which is theoretically not permitted by the analysis in the previous subsection. In order to guarantee convergence in such cases, the following trick may be applied.

In what follows, we consider the case when (α k 1 , β k 1 ) = (0, 0) and (α k 2 , β k 2 ) = (α k , β k ). Suppose that empirical convergence is observed for (α, β) = (1, 0). The idea is to adopt the following update rule for the inertial parameters: for any k, one set β k = 0 and

α k =    1 while (x k+1 , y k+1 ) -(x k , y k ) ≤ M (1 -ε) k 0.999 × 1 - τ L ∇F 2 otherwise ( 20
)
for given M > 0 and ε ∈ (0, 1). Then, two cases can occur:

(a) for any k, α k = 1. This means that (x k+1 , y k+1 ) -

(x k , y k ) ≤ M (1 -ε) k , that is, (x k , y k
) k∈N is a Cauchy sequence, thus it converges.

(b) there exists a k 0 such that

(x k0+1 , y k0+1 ) -(x k0 , y k0 ) ≤ M (1 -ε) k0
In this case, the choice for α k is feasible for k ≥ k 0 . Since the convergence results are asymptotic, this implies that all the convergence results stated in the previous section hold. In particular, under the proper hypotheses, (x k , y k ) k∈N converges.

Note that if the quantities M and ε are chosen so that M is large enough and 1 -ε is close enough to 1, the condition (x k+1 , y k+1 ) -(x k , y k ) ≤ M (1 -ε) k may hold for a sufficient number of iterations. This means that, in practice, one can always choose α k = 1 at least for the first iterations, without compromising the asymptotic convergence properties of the algorithm.

Numerical experiments

We recall that the studied colorization method consists in two steps: the candidate extraction step, which aims at extracting from the source image C color candidates for each pixel, and the candidate selection step via a vote process, which selects the best candidate under regularity assumptions, in the variational approach described in Section 2.

Candidate Extraction

As described in [START_REF] Bugeau | Patch-based image colorization[END_REF], for each pixel of the target image, a patch-based method is used to extract color candidates. Basically, the algorithm searches in the Y channel of the source image the closest patch with respect to a given criterion. The U and V chrominance channels corresponding to the pixel center of this retained patch is considered as a possible candidate. To speed-up the search of candidates, PatchMatch algorithm [START_REF] Barnes | Patchmatch: A randomized correspondence algorithm for structural image editing[END_REF] has been implemented in [START_REF] Pierre | Colociel. Dépôt Agence de Protection des Programmes No IDDN[END_REF]. The criterion used for the comparisons of the patches have been chosen experimentally, inspired from [START_REF] Bugeau | Patch-based image colorization[END_REF]. The number and the variety of the criterion have been reduced from the original version, because it has been experimentally observed in the work of [START_REF] Pierre | Colociel. Dépôt Agence de Protection des Programmes No IDDN[END_REF] that some criterion were especially low to compute and do not increase significantly the quality of the results. In the experiments, we have focused on four criterion:

• SSD (Euclidean norm) between patches with size 5×5;

• absolute difference between standard-deviation of patches with size 3×3 and 5×5;

• L1-norm between cumulative histograms (15 bins) of patches with size 11×11.

Applying these criterion provide four (C = 4) candidates. Let us remark that this candidate number is not restricted by the model of [START_REF] Pierre | Luminance-chrominance model for image colorization[END_REF] but by computational time considerations.

Candidate Selection

As discussed in Section 2, Model (1) fits with Model [START_REF] Bolte | Clarke subgradients of stratifiable functions[END_REF] with the following identifications:

F (u) = TV α C (u), G(w) = 0 and H(u, w) = λ 2 Ω C i=1 w i (x) u(x) -c i (x) 2 2 dx + χ R (u) + χ Σ (w)
In our implementation, the discrete version of ∇u = (∂ x U, ∂ y U, ∂ x V, ∂ y V ) used in the coupled total variation is defined thanks to horizontal (resp. vertical) forward finite differences, with symmetric boundary conditions. If Y is a noisy image, there exists α > 0 so that TV α C (u) = TV C (u) for any u. In practical case, when choosing α = 1, max γ ∇Y 2 , α = γ ∇Y 2 in a large part of pixels (for instance about 92% for Figure 1(b)). Thus this regularization of the coupled total variation has few influence on the result.

Hence, J satisfies Assumptions (H1). Moreover, one can check that, according to Example 2, Assumptions (H2) are fulfilled, while Assumptions (H3) hold when setting f = χ R and g = χ Σ . Eventually, Assumption (H4) is satisfied according to Remark 6. As a consequence, the constraints Σ and R being bounded, one can apply Algorithm 1 with a guaranteed strong convergence, as soon as the algorithm parameters satisfy Assumptions (A1)-(A4).

In order to compute the update of w, we use the Bregman distance associated to the entropy function to avoid an optimization on the simplex. The explicit solution is then obtained thanks to the KKT conditions. For the u-update, we consider the usual Moreau proximity operator (meaning that the associated Bregman distance is the Euclidean squared distance). With this setting, the inertial parameters associated to w are set to be null, and we choose a constant stepsize, which can be set arbitrarily large; we chose σ = 100. For the u-sequence, the inertial parameter sequence is chosen according to Subsection 5.2, namely equal to (α k 1 , α k 2 ) = (1, 1). The stepsize is then chosen to be large, e.g. τ = 1.999/L ∇F . Hence, the updates for u and w are explicitly given in Algorithm 2.

Algorithm 2 Inertial Bregman-based proximal gradient descent for image colorization

Input: Choose w 0 i (x) = 1/C for any x ∈ Ω and i = 1, • • • , C, u 0 (x) = C i=1 w 0 i (x)c i (x). Let u -1 = u 0 . for all k ≥ 0, do ūk = 2 u k -u k-1 , ûk = 2 u k -u k-1 , u k+1 = proj R       ûk -τ ∇TV α C (ū k ) + τ λ C i=1 w k i c i 1 + τ λ       ; w k+1 i (x) = w k i (x) exp   -σ λ C j=1 u k+1 (x) -c j (x) 2   C i=1 w k i (x) exp   -σ λ C j=1 u k+1 (x) -c j (x) 2   end for

Experimental Results and Discussion

In this section, some experimental results are proposed and an alternative optimization strategy is considered. The implementation is done in Matlab. In practical experiments, λ = 10 -3 and γ = 25 are reliable parameters for all images. The values of the limits of R are chosen as follows: U max = -U min = 111, and V max = -V min = 157.

From a visually qualitative point of view, let us remark that the results performed with Algorithm 2 are the same as with the primal-dual algorithm developed in [START_REF] Pierre | Luminance-chrominance model for image colorization[END_REF]. For instance, in Figure 1, we can see a source image, a target image and the results provided by Algorithm 2 and by [START_REF] Pierre | Luminance-chrominance model for image colorization[END_REF]. Visually, there is no difference and numerically, in average, the difference is pixel-wise less than 1/256 of the image range. Thus we can conclude that the local minimum computed by our proposed method is as reliable as the result computed by the method of [START_REF] Pierre | Luminance-chrominance model for image colorization[END_REF].

In Figure 2 the curves show the value of the functional with respect to the iteration of the algorithm. To see the benefit of using inertial scheme and Bregman distance, we also run the ASAP algorithm [START_REF] Nikolova | Alternating proximal gradient descent for nonconvex nonsmooth block-regularised problems[END_REF], which is a particular instance of the proposed method without inertia (i.e. null inertial parameters) and with usual proximity operator (i.e. the Bregman distances are the Euclidean distances). We also compared our method with the PALM method [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF][START_REF] Xu | A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion[END_REF], which reduces for this problem to an alternating projected gradient scheme. In our experiments, the projection onto the simplex is done using the KKT conditions. As expected, the proposed method (in purple) performs better that the ASAP algorithm (in yellow) in terms of energy convergence. The observed gain is entirely due to the use of inertia. Using Bregman distance to compute the update of the weights w has a notable incidence on the computational time: despite the addition of the overrelaxation steps for u, 300 iterations take 71.07s to compute for the proposed algorithm, while they take 81.32s for ASAP for the image tested in Figure 2. The PALM algorithm (in red) shows the slowest convergence of the bench, along with larger computational time (130.35s for the same image). The result provided by the primal-dual approach of [START_REF] Pierre | Luminance-chrominance model for image colorization[END_REF] is as reliable as the one computed by our Algorithm 2. In addition, with our proposed method, the convergence is guaranteed.

When comparing the energy evolution for the primal-dual approach of [START_REF] Pierre | Luminance-chrominance model for image colorization[END_REF] (in blue) and the proposed method (in yellow), one can see that Algorithm 2 makes the energy to decrease fast whereas the primal-dual algorithm may increase and oscillates in the first 100 iterations. Nevertheless, it loses this advantage in comparison of the primaldual approach after 330 iterations (see Figure 2(b)). The comparison of the histograms after 500 iterations (Figure 2(c) and (d)) shows that the local minimum is not reached by Algorithm 2 but more likely approached by the primal-dual algorithm, since it has been shown in [START_REF] Pierre | Luminance-chrominance model for image colorization[END_REF] that, at convergence, the minimizer w * contains only 0 and 1. Nevertheless, without any convergence guaranties, the primal-dual algorithm may have not converged. A way to benefit from the empirical convergence speed of the primal-dual algorithm while keeping the convergence guarantee provided by our proposed method is to consider the following strategy: first we use the primal-dual algorithm for 500 iterations and then we start Algorithm 2 from the values computed by the 500 first iterations of the primaldual algorithm. The final result is visually the same as the one of the two previous approaches (Algorithm 2, and primal-dual) and numerically, the difference between the two algorithms and our strategy is pixel-wise less than 1/256 of the image range in average. Figure 2: Value of the functional with respect to the number of iterations of the different algorithms (above) and the histograms of the weights after 500 iterations (bottom). In yellow, the energy of ASAP, in red, the one of [START_REF] Pierre | Luminance-chrominance model for image colorization[END_REF] and in blue, the one of our strategy.

In Figure 3, some additional results are available. One can see the source images (left column) the target ones (middle column) and the results computed with our strategy (right column). Visually, there is no difference with the results of [START_REF] Pierre | Luminance-chrominance model for image colorization[END_REF] but the convergence is guaranteed.

Some additional results can be found on http://www.fabienpierre.fr/colorisation/.

Conclusion

In this paper, we proposed a variant of an accelerated alternating proximal descent scheme introduced in [START_REF] Nikolova | Alternating proximal gradient descent for nonconvex nonsmooth block-regularised problems[END_REF] based on the Bregman distances. We provided a convergence analysis in the case when the coupling term is biconvex, along with a guide for the inertial parameter choice. This algorithm was applied to the image colorization problem. The numerical experiments confirm the empirical acceleration of the proposed scheme compared with the original one, as well as the benefit of the use of the Bregman distance in terms of computation times. Comparisons with two other optimization schemes, namely a Chambolle-Pock inspired primal-dual algorithm and the PALM algorithm prove that the proposed method gives comparable results as those obtained in the primal-dual framework with a guaranteed convergence (unlike the algorithm proposed in [START_REF] Pierre | Luminance-chrominance model for image colorization[END_REF]), while outperforming the PALM scheme, which is known as the state-of-the-art optimization scheme for nonconvex and nonsmooth problems. In future works, one could extend our algorithm to others biconvex problems such as joint image restoration and motion estimation [START_REF] Conti | A bv functional and its relaxation for joint motion estimation and image sequence recovery[END_REF].

Acknowledgement

This work has been partially funded by the French Research Agency (ANR) under grant No ANR-14-CE27-001 (MIRIAM) and supported by grants from Région Île-de-France.

Source image.

Target image Result

Figure 3: Results with our strategy. The result is first computed with the method of [START_REF] Pierre | Luminance-chrominance model for image colorization[END_REF] and refined with Algorithm 2. With this strategy, the speed of convergence is roughly the same as in [START_REF] Pierre | Luminance-chrominance model for image colorization[END_REF] and the convergence is theoretically guaranteed.

Proposition 1 [ 41 ,

 141 Theorem 10.1] Let f : R m → R ∪ {+∞} be a proper function. If f has a local minimum at x * , then 0 ∈ ∂f (x * ). Such points are called critical points, and the set of the critical points of f is denoted by critf .

Remark 2 Remark 3 Example 1

 231 Although we refer to D b as a distance, it does not satisfy the definition of a distance. Indeed, it does not fulfill the symmetry condition, since D b (x, y) = D b (y, x) in general. However, the other properties hold. Since b is convex, D b (x, y) is nonnegative for any (x, y) ∈ int(dom b) × dom b. The strict convexity of C also guarantees that D b (x, y) > 0 for x = y. For any x ∈ int(dom b), the function y → D b (x, y) is convex and continuous on dom b. Let us give some examples. (a) Let M be a definite positive matrix, and let x 2 M = x, M x . Then, D • M (x, y) = x -y 2 M /2 is the Bregman distance associated to • M . In particular, if M is the identity, one recovers the Euclidean squared distance. (b) For any m ∈ N * , let Σ m denotes the simplex of dimension m, defined by

  is lower bounded; (b) F : X → R and G : Y → R are continuously differentiable, and their gradients are Lipschitz continuous, of modulus L ∇F and L ∇G , respectively; (c) H : X × Y → R ∪ {+∞} is proper, l.s.c. and lower bounded.

Algorithm 1

 1 (a) If xk ∈ int(dom b X ), then x k+1 ∈ dom b X , and if ŷk ∈ int(dom b Y ), then y k+1 ∈ dom b Y . Inertial Bregman ASAP Input: b X : X → R ∪ {+∞} and b Y : Y → R ∪ {+∞} strictly convex, differentiable on int(dom b X ) and int(dom b Y ) respectively. Choose (x 0 , y 0 ) ∈ (int(dom b X ) × Y) ∩ dom J. Let x -1 = x 0 and y -1 = y 0 .

  )), and the stepsizes{(τ k , σ k )} k∈N (Assumptions (A3)).Assumptions (A1) (a) b X : X → R ∪ {+∞} and b Y : Y → R ∪ {+∞} are continuously differentiable on int(dom b X ) and int(dom b Y ), respectively; (b) b X and b Y are strongly convex of respective modulus L X and L Y .

  (a) Source image. (b) Target image. (c) Result of [34]. (d) Result of Algorithm 2.

Figure 1 :

 1 Figure1: The result provided by the primal-dual approach of[START_REF] Pierre | Luminance-chrominance model for image colorization[END_REF] is as reliable as the one computed by our Algorithm 2. In addition, with our proposed method, the convergence is guaranteed.

  (a) Evolution of the functionnal for the different alorithms. (b) Zoom on (a).

  Weights histogram of[START_REF] Pierre | Luminance-chrominance model for image colorization[END_REF].

E-mail: pauline.tan@cmla.ens-cachan.fr

E-mail: fabien.pierre@inria.fr

E-mail: nikolova@cmla.ens-cachan.fr

In this paper, we will refer to functions which only depend on one variable x or y as regularizers, while those which depends on both (x, y) are called coupling terms. Obviously, the additive decomposition of a functional into regularizers and coupling terms is not unique, since one always may add regularizers to a coupling term.

Note that the biconvexity of the whole objective J is not required in our analysis, and that for some weak convergence results, it is even not needed for the coupling term H.