Mathieu Chalvidal

Emilie Chouzenoux

Jean-Baptiste Fest
email: jean-baptiste.fest@centralesupelec.fr

Claire Lefort

Block Delayed Majorize-Minimize Subspace Algorithm for Large Scale Image Restoration

Keywords: Distributed optimization, asynchronous implementation, block alternating method, majorization-minimization, non-convex optimization, image restoration, depth-variant blur, multiphoton microscopy

Modern image acquisition devices, from microscopes to medical imaging machines, require to deal with increasingly large amount of data. To limit the dependence of an optimization algorithm on the dimension of the problem, distributed algorithms have been developed. In these schemes, at each iteration only a subset of the variables are updated simultaneously allowing to distribute computations on different nodes (or machines). The implementation of distributed algorithms requires to pay careful attention to the cost of communication, which can be reduced and controlled by resorting to an asynchronous implementation. However, asynchronous implementation raises challenging questions, in terms of convergence analysis, as the communication delays may introduce instabilities in the algorithm behavior. In this work, we propose an asynchronous majoration-minimization (MM) algorithm for solving large scale differentiable non-convex optimization problems. The proposed algorithm runs efficient MM memory gradient updates on block of coordinates, in a parallel and possibly asynchronous manner. We establish the convergence of the resulting sequence of iterates under mild assumptions. The performance of the algorithm is illustrated on the restoration of 3D images degraded by depth-variant 3D blur, arising in multiphoton microscopy. Significant computational time reduction, scalability and robustness are observed on synthetic data, when compared to state-of-the-art methods. Experiments on the restoration of real acquisitions of a muscle structure illustrate the qualitative performance of our approach and its practical applicability.

Introduction

Large-scale optimization algorithms, benefiting from fast convergence, capable of utilizing modern computing infrastructures, and dealing with distributed datasets are becoming compulsory for solving inverse problems in modern imaging [START_REF] Chung | Large-Scale Inverse Problems in Imaging[END_REF]. The evergrowing need for fast processing solutions that can operate on high-dimensional problems (i.e implying a huge number of variables) calls for the development of parallel methods harnessing the power of distributed computational architectures. In addition, the expansion of IoT systems and remote highly parallel computing induce new network issues with specific constraints. For instance, instabilities may occur whenever the volume of data dwarfs the memory capacity of a single agent or when the processing power is shared (potentially unevenly) between devices [START_REF] Konečný | Federated optimization: Distributed machine learning for on-device intelligence[END_REF]. Several classes of so-called distributed optimization methods, have been investigated under various assumptions on the computing scenario and on the optimization problem itself, that we review hereafter (see also [START_REF] Zheng | A review of distributed optimization: Problems, models and algorithms[END_REF][START_REF] Yang | A survey of distributed optimization[END_REF]).

Distributed optimization approaches inherit from block alternating methods. In the latter, at each iteration, only a subset of the variables are updated, by minimizing the objective function with respect to only those variables, the others being fixed. The blocks are selected sequentially following a cyclic (or quasi-cyclic) order or a random rule. Exact minimization with respect to a given block of variables is rarely possible in a closed form. It is not even desirable as it may lead to convergence issues [START_REF] Tseng | Convergence of a block coordinate descent method for nondifferentiable minimization[END_REF]. More efficient and stable block alternating schemes rely on a so-called majoration-minimization (MM) strategy [START_REF] Jacobson | An expanded theoretical treatment of iteration-dependent majorize-minimize algorithms[END_REF]. It consists in building, at each iteration, a majorizing approximation for the objective function within the active block of variables, whose minimizer has a more tractable form. Many powerful algorithms fall within this framework, such as BSUM [START_REF] Hong | A unified algorithmic framework for blockstructured optimization involving big data: With applications in machine learning and signal processing[END_REF], PALM [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF], NMF [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF], to name a few. By exploiting the structure of the objective function, block alternating MM methods can reach fast convergence rate [START_REF] Duval | Euclid in a taxicab: Sparse blind deconvolution with smoothed l1/l2 regularization[END_REF][START_REF] Fessler | Grouped coordinate descent algorithms for robust edge-preserving image restoration[END_REF][START_REF] Nesterov | Efficiency of coordinate descent methods on huge-scale optimization problems[END_REF][START_REF] Nutini | Coordinate descent converges faster with the Gauss-Southwell rule than random selection[END_REF] while offering theoretical guarantees in non-convex cases [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF][START_REF] Chouzenoux | A block coordinate variable metric forwardbackward algorithm[END_REF][START_REF] Bonettini | A block coordinate variable metric linesearch based proximal gradient method[END_REF].

When the problem size increases, as in 3D microscopy imaging [START_REF] Lefort | FAMOUS: a fast instrumental and computational pipeline for multiphoton microscopy applied to 3d imaging of muscle ultrastructure[END_REF] and astronomy [START_REF] Onose | Scalable splitting algorithms for big-data interferometric imaging in the SKA era[END_REF][START_REF] Prato | A blind deconvolution method for ground based telescopes and Fizeau interferometers[END_REF], running block alternating methods gets inefficient. Parallel implementations have been devised, where the block updates are performed simultaneously, allowing to distribute computations on different nodes (or machines) [START_REF] Cadoni | A block parallel majorize-minimize memory gradient algorithm[END_REF][START_REF] Sotthivirat | Image recovery using partitioned-separable paraboloidal surrogate coordinate ascent algorithms[END_REF][START_REF] Richtárik | Parallel coordinate descent methods for big data optimization[END_REF]. Implementation on parallel architecture requires to pay attention to communication cost. The latter can be reduced by resorting to an asynchronous parallel implementation, yielding the so-called distributed optimization approach. Each computation node has now its own iteration loop, so local variables are updated without the need to wait for distant variables update. This however raises challenging questions, in terms of convergence analysis, as the communication delays may introduce instabilities. A plethora of recent works have focused on proposing distributed optimization algorithms with assessed convergence, based on stochastic proximal primal [START_REF] Grishchenko | Asynchronous distributed learning with sparse communications and identification[END_REF][START_REF] Lian | Asynchronous parallel stochastic gradient for nonconvex optimization[END_REF][START_REF] Mishchenko | A distributed flexible delay-tolerant proximal gradient algorithm[END_REF] or primal-dual [START_REF] Pesquet | A class of randomized primal-dual algorithms for distributed optimization[END_REF][START_REF] Hannah | On unbounded delays in asynchronous parallel fixed-point algorithms[END_REF][START_REF] Zhang | Asynchronous distributed ADMM for consensus optimization[END_REF][START_REF] Chorobura | Random coordinate descent methods for nonseparable composite optimization[END_REF][START_REF] Onose | Scalable splitting algorithms for big-data interferometric imaging in the SKA era[END_REF][START_REF] Abboud | Distributed algorithms for proximity operator computation with applications to video processing[END_REF] techniques. Recent contributions in the field of federated learning are also highly related [START_REF] Horváth | Fedshuffle: Recipes for better use of local work in federated learning[END_REF][START_REF] Loizou | Revisiting randomized gossip algorithms: General framework, convergence rates and novel block and accelerated protocols[END_REF][START_REF] Wang | A field guide to federated optimization[END_REF]. However, as the aforementioned works rely on specific fixed-point analysis tools involving Fenchel-Rockafellar duality [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF], the proposed algorithms are limited to convex (sometimes even strongly convex) optimization and often require specific probabilistic assumptions on the block update rule difficult to meet in practice. In the context of MM algorithms, although the need for distributed implementation strategies is crucial (see the discussion in [START_REF] Hong | A unified algorithmic framework for blockstructured optimization involving big data: With applications in machine learning and signal processing[END_REF] and the specific examples in [START_REF] Tuck | Distributed majorization-minimization for Laplacian regularized problems[END_REF][START_REF] Fan | Majorization minimization methods for distributed pose graph optimization with convergence guarantees[END_REF]), theoretical results regarding convergence guarantees of MM technique in a distributed context are rather scarce. Let us first mention the work of [START_REF] Davis | The asynchronous PALM algorithm for nonsmooth nonconvex problems[END_REF][START_REF] Davis | The sound of APALM clapping: Faster nonsmooth nonconvex optimization with stochastic asynchronous palm[END_REF], that proposes an asynchronous version of PALM, with proven convergence of the iterates in non-convex case, and good practical behaviour [START_REF] Thouvenin | Partially asynchronous distributed unmixing of hyperspectral images[END_REF]. The convergence of distributed MM methods was also explored in the recent works [START_REF] Cannelli | Asynchronous optimization over graphs: Linear convergence under error bound conditions[END_REF][START_REF] Li | Surrogate-based distributed optimisation for expensive black-box functions[END_REF]. However, the analysis of [START_REF] Li | Surrogate-based distributed optimisation for expensive black-box functions[END_REF] is limited to the convex case. In [START_REF] Cannelli | Asynchronous optimization over graphs: Linear convergence under error bound conditions[END_REF], the analysis covers non-convex terms in the objective function, but it only shows the convergence of the sequence of objective function values, and not the convergence of the iterates themselves (thus, the results is weaker than the one of [START_REF] Davis | The asynchronous PALM algorithm for nonsmooth nonconvex problems[END_REF]).

In this work, we aim at solving a smooth optimization problem of the form minimize

xPR N f pxq, (1)
where f : R N Þ Ñ R is (Fréchet) differentiable but non necessarily convex. In the context of inverse problems in imaging, f typically reads as the sum of a data fidelity term (e.g., a least-squares term) measuring the discrepancy between an acquired, degraded (e.g., blurry, noisy) image, and its estimate (usually, through a linear observation operator), and a regularization term incorporating prior information on the sought solution [START_REF] Chung | Large-Scale Inverse Problems in Imaging[END_REF][START_REF] Bertero | Introduction to Inverse Problems in Imaging[END_REF] (see also our Section 5). We introduce the block delayed MM memory gradient (BD3MG) algorithm for the resolution of Problem [START_REF] Abboud | Distributed algorithms for proximity operator computation with applications to video processing[END_REF]. BD3MG is a distributed MM algorithm designed for an efficient implementation on a multi-CPU computing architecture, such as a high performance calculation unit. Our contributions are::

' Introduction of the BD3MG algorithm, that implements an advanced distributed asynchronous update rule within the block alternating MM method we recently proposed in [START_REF] Fest | Convergence analysis of block majorize-minimize subspace approaches[END_REF].

' Proof for the convergence of BD3MG iterates to a stationary point of f under mild assumptions (in particular, no convexity is assumed), using recent tools of Lyapunov analysis [START_REF] Wilson | Lyapunov Arguments in Optimization[END_REF].

' Illustration of the performance of BD3MG by means of various experiments on a real inverse problem of 3D image restoration arising in the context of multiphotonic miscroscopy.

The paper is organized as follows. Section 2 introduces our notations, recalls the principle of MM schemes and finally presents our proposed algorithm. Section 3 states our mathematical assumptions for the convergence analysis and presents preliminary : A preliminary version of this work has been presented in the conference proceedings [START_REF] Chalvidal | Block distributed 3MG algorithm and its application to 3D image restoration[END_REF]. The convergence result was weaker, and stated without proof. The experimental validation was limited to a single, simpler, numerical scenario.

technical propositions and lemmas. Section 4 presents our main theoretical contribution, dedicated to the convergence analysis of the proposed BD3MG scheme. Section 5 illustrates the qualitative and computational performance of BD3MG in the applicative context of 3D image deblurring in the presence of a depth-variant 3D blur. Section 6 concludes the paper.

Proposed algorithm

Notations

Throughout the paper, we consider the euclidean space R N endowed with the usual scalar product x¨|¨y (or, equivalently, ¨J¨) and the norm } ¨}. 0 N is the vector with null entries of R N . I N is the identity matrix of R N . We use the short notation rr1, N ss, to denote t1, 2, . . . , N u, i.e. the set of integers from 1 to N . S N denotes the set of symmetric matrices of R N ˆN , and S N `(resp. S N ``) the set of positive (resp definite positive) symmetric matrices. Given some M P S N ``, we denote by } ¨}M the induced weighted euclidean norm, such that, for all v P R N , }v} 2 M " v J M v. We use the Loewner orders symbols ă and ĺ, to compare real symmetric matrices pA, Bq P pS N q 2 i.e., A ĺ B (resp. A ă Bq is verified when difference B ´A belongs to S N `(resp. S N

``).

Let us introduce extra notations, that will be useful to present block coordinate optimization strategy. Most notations hereafter are reminiscent from [START_REF] Fest | Convergence analysis of block majorize-minimize subspace approaches[END_REF]. Let S Ă rr1, N ss. Ź We denote by S its complementary set rr1, N sszS, |S| its cardinal and ´R|S| , x¨, ¨yt he resulting euclidean space (with a slight abuse of notation). Moreover, we also denote by S |S| , S |S| `, S |S| ``respectively the set of symmetric, symmetric positive, and symmetric definite positive matrices of R |S|ˆ|S| . Ź Let x " px i q iPrr1,N ss P R N . We denote x pSq " px i q iPS P R |S| the vector gathering the entries of x with indexes within the set S of coordinates.

Ź Let x P R N . ∇f pxq is the gradient of f evaluated at x. Moreover, ∇ pSq f pxq "

´"∇f pxq ‰ i ¯iPS P R |S| denotes the partial gradient of f with respect to the coordinates with indexes in S, evaluated at x. Ź Let M P S N . We denote the (symmetric) sub-matrix M pSq " `Mi,i ˘iPS P S |S| . If M pSq P S |S| `, we define the induced weighted euclidean norm as } ¨}M pSq . Ź For any x P R N , we introduce the restriction of f to the block S and vector x as the function f pSq p. , xq : v P R |S| Þ Ñ f puq where u is related to pv, xq through the relations u pSq " v and u pSq " x pSq .

Block MM principle

MM approach for the resolution of Problem ( 1) is a generic iterative procedure where each iteration amounts to minimizing (exactly or not) a surrogate for f satisfying a majorizing property [START_REF] Sun | Majorization-minimization algorithms in signal processing, communications, and machine learning[END_REF][START_REF] Zhang | Surrogate maximization/minimization algorithms and extensions[END_REF][START_REF] Jacobson | An expanded theoretical treatment of iteration-dependent majorize-minimize algorithms[END_REF][START_REF] Hunter | A tutorial on MM algorithms[END_REF]. The theoretical and practical properties of an MM algorithm greatly depend on (i) the family of considered surrogates, (ii) the procedure to minimize it. In this work, we focus on quadratic MM techniques, where f is such that it can be upper bounded by quadratic functions (typically, f is Lipschitz differentiable). In such context, the inner step of an MM algorithm amounts to minimizing a quadratic function on R N or, otherwise stating, to invert an N ˆN system. In the large scale context, this is not desirable and various approaches have been proposed to cope with the curse of dimensionality in MM quadratic methods [START_REF] Chouzenoux | A majorize-minimize strategy for subspace optimization applied to image restoration[END_REF][START_REF] Chouzenoux | A majorize-minimize subspace approach for ℓ 2 ´ℓ0 image regularization[END_REF][START_REF] Chouzenoux | A stochastic majorize-minimize subspace algorithm for online penalized least squares estimation[END_REF][START_REF] Hong | A unified algorithmic framework for blockstructured optimization involving big data: With applications in machine learning and signal processing[END_REF][START_REF] Cadoni | A block parallel majorize-minimize memory gradient algorithm[END_REF][START_REF] Sotthivirat | Image recovery using partitioned-separable paraboloidal surrogate coordinate ascent algorithms[END_REF][START_REF] Fest | Convergence analysis of block majorize-minimize subspace approaches[END_REF]. In particular, to limit the dependence of the MM algorithm on the dimension of the problem, block alternating approaches have been developed. In these schemes, at each iteration only a subset of the variables are updated [START_REF] Hong | A unified algorithmic framework for blockstructured optimization involving big data: With applications in machine learning and signal processing[END_REF], giving rise to so-called block MM algorithms, that we describe hereafter.

Define a partition T of rr1, N ss. Let x P R N , and some block index S P T. The block MM approach requires to build a majorizing surrogate for the restriction f pSq p¨, xq. Let us assume the existence of a mapping A :

x P R N Þ Ñ Apxq P S N ``such that p@x P R |S| q Q S pv, xq " f pxq `x∇ pSq f pxq, v ´xpSq y `1 2 › › v ´xpSq › › 2 A pSq pxq , (2)
fulfills the majorizing condition p@S P Tqp@v P R |S| qp@x P R N q Q S pv, xq ě f pSq pv, xq

Note that, by [START_REF] Allain | On global and local convergence of half-quadratic algorithms[END_REF],

Q S px pSq , xq " f pSq px pSq , xq. (4)
The existence of such mapping A can be ensured under mild assumptions. For instance, it is satisfied as soon as f is Lipschitz differentiable. Morever, [START_REF] Chouzenoux | A block coordinate variable metric forwardbackward algorithm[END_REF]Remark 2.4] shows that, as soon as the above mapping holds for S " R N , it stays valid for any block subset S Ă rr1, N ss. Examples of constructions of majorant mappings have been extensively discussed in [START_REF] Sun | Majorization-minimization algorithms in signal processing, communications, and machine learning[END_REF][START_REF] Chouzenoux | A majorize-minimize subspace approach for ℓ 2 ´ℓ0 image regularization[END_REF][START_REF] Sotthivirat | Image recovery using partitioned-separable paraboloidal surrogate coordinate ascent algorithms[END_REF] for optimization problems arising in the fields of inverse problems, image processing and telecommunication.

Once the block majorant approximations (2) satisfying (3) are built, the block MM (B2M) algorithm reads [START_REF] Jacobson | An expanded theoretical treatment of iteration-dependent majorize-minimize algorithms[END_REF] (also called BSUM in [START_REF] Hong | A unified algorithmic framework for blockstructured optimization involving big data: With applications in machine learning and signal processing[END_REF]):

p@k P Nq $ ' ' ' ' & ' ' ' ' % Choose S k P T, x k`1 pS k q P arg min vPR |S k | Q S k pv, x k q, x k`1 pS k q " x k pS k q . (B2M)
Hereabove, `Sk ˘kPN is a sequence of subsets (i.e., blocks) of T whose construction is carried out upstream. The most current strategy is to adopt a cyclic rule, where each element of T is selected sequentially until the end of the partition list, and then the loop is repeated until convergence of the algorithm. A more flexible option is to adopt a so-called quasi-cyclic (or acyclic) rule where each S P T must be updated at least once per K iterations period. The interest of scheme (B2M) and more generally block coordinate methods notably lies in the large scale context involving a very huge N , for which dealing all coordinates of the current iterate may be too high time consuming and even infeasible for memory reasons. However, block MM methods require a sequential update of the blocks and thus, by construction, might require many iterations to reach convergence. To limit this issue, (block) diagonal mappings have been considered for instance in [START_REF] Sotthivirat | Image recovery using partitioned-separable paraboloidal surrogate coordinate ascent algorithms[END_REF][START_REF] Cadoni | A block parallel majorize-minimize memory gradient algorithm[END_REF]. The underlying idea is to choose the mapping so that the inner minimization problem in (B2M) is separable, and thus can be performed in parallel over the entries of the selected block. This yields the so-called block parallel MM schemes that take advantage of recent technological advances in parallel computing on multicore architectures. In particular, these methods can tailor the number of available processors to the computational load. However, such block diagonal structure may be detrimental to the approximation quality of the surrogates, and thus reduce again the practical convergence rate. In the present work, we opt for not making any extra structural assumption on the majorant mappings, thanks to the introduction of two catalizers into (B2M), namely (i) a subspace acceleration approach, (ii) a distributed asynchronous update strategy, that we describe hereafter.

Subspace acceleration

Our first catalyst is to introduce a subspace acceleration [START_REF] Sun | Majorization-minimization algorithms in signal processing, communications, and machine learning[END_REF], in (B2M). This strategy has been initially introduced for full-batch MM algorithms (i.e., without any block coordinate strategy) in [START_REF] Chouzenoux | A majorize-minimize strategy for subspace optimization applied to image restoration[END_REF]. Convergence analysis can be found in [START_REF] Chouzenoux | A majorize-minimize subspace approach for ℓ 2 ´ℓ0 image regularization[END_REF][START_REF] Chouzenoux | SABRINA: a stochastic subspace majorization-minimization algorithm[END_REF][START_REF] Chouzenoux | A stochastic majorize-minimize subspace algorithm for online penalized least squares estimation[END_REF][START_REF] Chouzenoux | Convergence rate analysis of the majorize-minimize subspace algorithm[END_REF][START_REF] Chouzenoux | A local MM subspace method for solving constrained variational problems in image recovery[END_REF] under various situations. We recently extended this strategy to cope with block coordinate updates with the form of (B2M) [START_REF] Fest | Convergence analysis of block majorize-minimize subspace approaches[END_REF], leading to the B2MS (Block MM Subspace) scheme that we present hereafter.

Starting with the (B2M) iteration, the subspace acceleration resort to performing the minimization of the majorant function within the current block S k in a constrained vectorial subspace spanned by a small number M k ě 1 of search directions. This reads:

p@k P Nq $ ' ' ' ' ' ' ' ' ' & ' ' ' ' ' ' ' ' ' % Choose S k P T, Set D k P R M k ˆ|S k | , v k P arg min vPR M k Q S k px k pS k q `Dk v, x k q, x k`1 pS k q " w k pS k q `Dk v k , x k`1 pS k q " x k pS k q , . (B2MS)
Hereabove, for every k P N, D k P R M k ˆ|S k | is the so-called subspace matrix. It stacks, row-wise, M k ě 1 vectors of dimension |S k |, spanning a vectorial subspace within which we seek for a minimizer of the majorant function Q S k p. , x k q (i.e., our next iterate). The advantage is to reduce again the dimensionality of the inner MM problems, without jeopardizing the convergence rate [START_REF] Chouzenoux | Convergence rate analysis of the majorize-minimize subspace algorithm[END_REF]. Several choices for the subspace matrix are discussed in [START_REF] Chouzenoux | A majorize-minimize strategy for subspace optimization applied to image restoration[END_REF][START_REF] Chouzenoux | A stochastic majorize-minimize subspace algorithm for online penalized least squares estimation[END_REF][START_REF] Chouzenoux | A local MM subspace method for solving constrained variational problems in image recovery[END_REF]. Intensive comparisons in the fields of inverse problems, image processing and machine learning (e.g., [START_REF] Florescu | A majorize-minimize memory gradient method for complex-valued inverse problems[END_REF][START_REF] Chouzenoux | SABRINA: a stochastic subspace majorization-minimization algorithm[END_REF]), have shown the superiority of the socalled memory gradient subspace which seems to reach the best compromise between simplicity and efficiency. In the context of (B2MS), this amounts to defining, for every k P N, the memory gradient matrix D k " " ∇ pS k q f px k q, x k pS k q ´xk´1 pS k q ı (with the convention x ´1 " 0 N ), so that M k " 2. When combined with a block diagonal majorant mapping, (B2MS) becomes equivalent to the BP3MG method considered in [START_REF] Chalvidal | Block distributed 3MG algorithm and its application to 3D image restoration[END_REF] for 3D image deblurring. The convergence properties of (B2MS) have recently been studied in [START_REF] Fest | Convergence analysis of block majorize-minimize subspace approaches[END_REF].

Block Delayed Majorize-Minimize Memory Gradient (BD3MG)

The second catalyst we introduce is the main contribution of this paper, namely the introduction of a distributed asynchronous update rule within (B2MS). Our motivation is to make the algorithm well suited to an implementation on a multi-core / multiprocessor architecture, while not being endangered by potential communication delays within the computing units. Let us consider a computing architecture with C units (or cores), each of them being able to communicate (i.e., send or receive) information to a master node. The architecture thus considered is forming a star graph as presented in Figure (1c). The two other graph topologies are discarded from this present study (see, for example, [START_REF] Abboud | Distributed algorithms for proximity operator computation with applications to video processing[END_REF] for an efficient distributed method running on a generic hypergraph topology). The proposed method BD3MG is presented in Algorithms 1-2, describing the iterations of the master (i.e., node 0) and a given worker/node c P rr1, Css, respectively. Let us describe these two algorithms. Each computation node c P rr1, Css has its own iteration loop, so that it can keep updating its local variables without having to wait for the update of the other, distant variables. It updates (independently from the other nodes) a subset of coordinates S c P T (which can change over the process) by applying an MM iteration including a memory gradient acceleration. Each node "books" its running block S c so that no other worker overwrites the associated coordinates. Conversely, any other S P TztS c u remains free to be updated by other workers. Communication steps are performed in order to maintain convergence to a minimizer of the globally shared objective function f and to control the propagation of errors. Basically, even if the other workers are still busy on their tasks, every time a worker c P rr1, Css ends one MM iteration on its running block S c , it sends a feedback to the master. As a response, the latter updates it with most recent available information, and assigns it a new task.

We denote px k q kPN the sequence of iterates gathered by the master loop. For any given node index c P rr1, Css and k P N, S k c denotes the block of coordinates processed by worker c during step k. We impose, by construction, that two nodes do not update the same block of coordinates at the same time, so that we ensure the no-overlap condition p@k P Nq p@pc, c

1 q P rr1, Css 2 q S k c X S k c 1 " H. (5)
At iteration k P N, worker c k P rr1, Css, updates the block of coordinates S k c k and sends to the master a vector d k of size ˇˇS k c k ˇˇ. The corresponding indexes of variable x k within block S k c k are then incremented with d k while the others remain unchanged, thus defining x k`1 . The master then defines a new set of coordinates S k`1 c k to be treated by worker c k , so as to satisfy the non-overlap rule. The master informs worker c k of this new running set of coordinates, and sends him the most recent information x k`1 and the difference px k ´xk´1 q pS k`1 c k q . Meanwhile, the other workers keep processing their allocated indexes. The master then waits until a new worker (possibly the same one) sends a new increment.

Let us now make a focus on the worker loop described in Algorithm 2. Remark that, even if worker c has access to some properties of function f (i.e., the expression for its gradient and for its majorizing approximation pQ S q SPS ), it is not informed of all the path of the master iterates. It can only rely on the data the master sends to it to perform its local task. From the viewpoint of the worker, a triplet set px, S, dq P R N ˆT ˆR|S| is received from the master and must be used to perform its MM update with memory gradient acceleration. The worker is in charge of first building the new memory gradient matrix Dpx, S, dq " r´∇ pSq f pxq | ds P R |S|ˆ2 .

and then compute the MM increment

d 1 P R |S| defined as d 1 " Dpx, S, dqu (7)
with u P arg min

vPR 2 Q S px `Dpx, S, dqv, xq. (8)
Note that the uniqueness of the solution for problem ( 7)-( 8) is not guaranteed in general.

To overcome such an obstacle, we follow the strategy in [START_REF] Chouzenoux | A majorize-minimize strategy for subspace optimization applied to image restoration[END_REF], and retain the pseudoinverse solution given by u " ´´Dpx, S, dq J A pSq pxq Dpx, S, dq ¯: Dpx, S, dq ∇ pSq f pxq,

Algorithm 1. BD3MG algorithm -Master loop Initialization.

(a) Set k " 0 and x 0 P R N . (b) Set S 0 1 , ..., S 0 C P T such that @pc, c

1 q P rr1, Css 2 , S 0 c X S 0 c 1 " H. (c) 0-th transmission: For every c P rr1, Css, send `x0 , S 0 c , 0 |S 0 c | ˘to worker c While a stopping criterion is not met: (Wait for a feedback from a worker) (a) pk `1q-th reception: Receive d k from a worker c k . (b) Update x k`1 ´Sk c k ¯" x k ´Sk c k ¯`d k and x k`1 ´Sk c k ¯" x k ´Sk c k ¯. (c) Set S k`1 1 , . . . , S k`1 C P T such that S k`1 c k P TztS k c u c‰c k and, p@c P rr1, Csszc k q, S k`1 c " S k c . (d) pk `1q-th transmission: Send ˆxk`1 , S k`1 c k , px k`1 ´xk q ´Sk`1 c k ¯ṫo worker c k . (e) k " k `1 End While Output. Vector x k .
where : referes to the Moore-Penrose pseudo-inverse. Such solution notably verifies the normal equation (e) Send d 1 " Dpx, S, dqu to the Master.

@ ∇ pSq f pxq, Dpx, S, dqu D " ´}Dpx, S, dqu} 2 A pSq pxq . (10

Equivalent form for BD3MG

The way we introduced our scheme BD3MG in the previous subsection was "implementation-oriented". In order to study its convergence behaviour, we must exhibit an equivalent form of it, mimicking the one of its non distributed counterpart, (B2MS). To do so, it is necessary to formalize the information gap between the master and the workers during the iterative process.

As we have already mentioned, all the information available to a worker (except those on f and pQ S q SPS ) is sent to it by the master only after it produces a feedback. For a given k P N, worker c k does not receive any information between the pk `1q-reception and the previous one it made. During this time, its counterparts c P 1, C ztc k u may have performed additional updates to the master without c k being informed. This results in an information mismatch, that we propose to formalize through a vector x ι k where

p@k P Nq ι k " $ & % 0 if k " 0, max ´␣ℓ P rr1, kss | c ℓ´1 " c k ( Y t0u ¯, otherwise. (11)
This vector corresponds to the iteration index of the working variable of worker c k , which does not necessarily matches with the vector x k manipulated by the master.

Let us list herebelow some situations of interest given the value of ι k at some iteration k P N:

' If ι k " 0, and k ą 0, it means that

␣ ℓ P rr1, kss | c ℓ´1 " c k (
is an empty set. Hence, the worker c k never returned any feedback to the master before the iteration k. Note that ι 0 " 0 by construction.

' If ι k " k, we thus have c k´1 " c k . Hence, worker c k was in charge of the two most recent updates, namely the pk `1q-th and the k-th ones. As a consequence, to prepare the pk `1q-th update, worker c k received vector x k from the master.

' More generally, if ι k ą 0, it follows that worker c k at least returned one feedback to the master before iteration k. And we have the relation

c ι k ´1 " c k .
Moreover, the non-overlap rule translates into

p@k P Nq x ι k pS k c k q " x k pS k c k q . (12)
In particular, when k ą 1 and ι k " k ´1, we have c k´2 " c k . Worker c k has proceeded to the pk ´1q-th and pk `1q-th master's reception. The k-th reception was made some another worker ck who received vector x k from the master. However, since worker c k was still proceeding block S k´1 c k , the master was not able to update the associated coordinate for computing x k from x k´1 for worker c k , i.e x k´1

pS k c k q " x k pS k c k q .
More generally, when it comes with dealing with asynchronous algorithms, the use of a specific indexes with similar roles than our ι k (k P N) is often necessary to build a theoretical delay model and thus to formulate an equivalent scheme being more compact and easier to analyse [START_REF] Davis | The asynchronous PALM algorithm for nonsmooth nonconvex problems[END_REF].

With this aim in mind, let us introduce the shorter notations

p@k P Nq $ & % B k " S k c k , D k " D ´xι k , B k , `xι k ´xι k ´1˘p B k q ¯, (13)
and

D k " D ´xι k , B k , `xι k ´xι k ´1˘p
B k q ¯with convention x ´1 " 0 N . Then, the master/worker BD3MG loops from Algorithms 1-2 can be rewritten equivalently in a single compact scheme as:

p@k P Nq $ ' ' ' ' ' ' & ' ' ' ' ' ' % Let c k P rr1, Css, u k " ´´pD k q J A pB k q px ι k q D k ¯: pD k q J ∇ pB k q f px ι k q, x k`1 pB k q " x k pB k q `Dk u k , x k`1 pB k q " x k pB k q , (14)
where we noticed that Equality (12) now reads (using ( 13))

p@k P Nq x ι k pB k q " x k pB k q . ( 15
)
For every k P N, according to ( 14), we have the optimality equation:

´∇pB k q f px ι k q ¯J ´xk`1 pB k q ´xk pB k q ¯" ´› › ›x k`1 pB k q ´xk pB k q › › › 2 A pB k q px ι k q . ( 16
)
The two next Sections are dedicated to establish the convergence of the iterates produced by ( 14).

Link with existing works

Let us discuss the links between our proposed scheme BD3MG and existing methods from the literature. When ι k " k for any k P N in BD3MG, the algorithm identifies with our block alternating scheme B2MS [START_REF] Fest | Convergence analysis of block majorize-minimize subspace approaches[END_REF] where the blocks of variables were updated sequentially in a non parallel (thus, not asynchronous) manner. This present paper can thus be viewed as an extension of the framework and of the convergence analysis of [START_REF] Fest | Convergence analysis of block majorize-minimize subspace approaches[END_REF] to the distributed setting. Other related methods are [START_REF] Davis | The asynchronous PALM algorithm for nonsmooth nonconvex problems[END_REF][START_REF] Cannelli | Asynchronous optimization over graphs: Linear convergence under error bound conditions[END_REF][START_REF] Li | Surrogate-based distributed optimisation for expensive black-box functions[END_REF], and our convergence analysis relies on similar tools than the one from [START_REF] Davis | The asynchronous PALM algorithm for nonsmooth nonconvex problems[END_REF]. Assuming zero-valued nonsmooth terms in [START_REF] Davis | The asynchronous PALM algorithm for nonsmooth nonconvex problems[END_REF][START_REF] Cannelli | Asynchronous optimization over graphs: Linear convergence under error bound conditions[END_REF][START_REF] Li | Surrogate-based distributed optimisation for expensive black-box functions[END_REF] (i.e., the objective function is differentiable), these methods identify with particular instances of BD3MG that (i) would not implement any subspace acceleration (i.e., D k " I N in ( 14)), (ii) would rely on the simple Lipschitz-based majorant metric (i.e., A pB k q px ι k q " LI |B k | in ( 14)) in the case of [START_REF] Davis | The asynchronous PALM algorithm for nonsmooth nonconvex problems[END_REF]. As a consequence, assuming differentiability of all terms, our convergence analysis presented in the next section thus also covers the schemes of [START_REF] Davis | The asynchronous PALM algorithm for nonsmooth nonconvex problems[END_REF][START_REF] Cannelli | Asynchronous optimization over graphs: Linear convergence under error bound conditions[END_REF][START_REF] Li | Surrogate-based distributed optimisation for expensive black-box functions[END_REF]. Up to our knowledge, our work is the first to show convergence of the iterates of a distributed MM algorithm involving generic quadratic surrogates and subspace acceleration, in the non-convex setting. Finally, we would like to point out that the 3MG update performed in Alg. 2 identifies with a nonlinear conjugate gradient (NLCG) update, for a specific (and closed form) pair of stepsize and conjugacy parameters (see discussion in [START_REF] Chouzenoux | A local MM subspace method for solving constrained variational problems in image recovery[END_REF]Sec. 1]). Therefore, our work can also be understood as the first convergence analysis of a distributed NLCG method in the non-convex setting. A comparative study will be conducted in our experimental section to illustrate the superiority of BD3MG with respect to the aforementioned existing methods in terms of convergence speed.

Assumptions and preliminary results

Assumptions

In order to analyse the asymptotic behaviour of the sequence px k q kPN generated by scheme [START_REF] Chorobura | Random coordinate descent methods for nonseparable composite optimization[END_REF], we introduce technical assumptions both on function f and on the parameters of BD3MG method. Assumption 1. Function f is coercive, continuously differentiable on R N , and has a L-Lipschitzian gradient with L ą 0, i.e. p@px, yq P pR N q 2 q }∇f pxq ´∇f pyq} ď L}x ´y}.

Assumption 1 ensures the existence of solutions for Problem (1) (by coercivity). Moreover, [START_REF] Chouzenoux | A majorize-minimize subspace approach for ℓ 2 ´ℓ0 image regularization[END_REF] in Assumption 1 guarantees the existence of a quadratic function ( 2) satisfying (3), setting A : x Þ Ñ LI N . Another direct consequence is p@S P Tqp@px, yq P pR N q 2 q }∇ pSq f pxq ´∇pSq f pyq} ď L}x ´y},

since }∇ pSq f pxq ´∇pSq f pyq} ď }∇f pxq ´∇f pyq} for all S P T and px, yq P pR N q 2 . Assumption 2. (Boundedness of delay) For every k P N, and every worker c k P rr1, N ss, the set S k c is not empty and there exists τ P N ˚such that p@k ě τ q rr1, N ss "

k´1 ď j"k´τ B j , (19)
using the notation in [START_REF] Chalvidal | Block distributed 3MG algorithm and its application to 3D image restoration[END_REF].

Assumption 2 gives an upper bound on the delay τ . Each of the block of variables should be updated within a time frame of at most τ iterations and thus the workers must follow a certain regularity. Such a condition follows a similar goal than quasicyclic rule frequently assumed in block coordinate methods [START_REF] Fest | Convergence analysis of block majorize-minimize subspace approaches[END_REF][START_REF] Hong | A unified algorithmic framework for blockstructured optimization involving big data: With applications in machine learning and signal processing[END_REF]. From Assumption 2, we deduce the following Proposition, which appears fundamental for the rest of our convergence study. It guarantees that, for a given k P N, the vector treated by worker c k before its feedback (i.e the pk `1q-th master's reception) is not "too old" compared to the iteration index. Proposition 3.1. Under Assumption 2, for every k ě τ , the index ι k given in [START_REF] Cannelli | Asynchronous optimization over graphs: Linear convergence under error bound conditions[END_REF] belongs to rrk ´τ `1, kss.

Proof. Let k ě τ , where τ ą 0 is defined in Assumption 2. Inequality ι k ď k directly comes from Definition [START_REF] Cannelli | Asynchronous optimization over graphs: Linear convergence under error bound conditions[END_REF]. We prove the lower bound on ι k by contradiction. Let us suppose that ι k ď k ´τ . Two situations may arise.

p@j P rr0, k ´1ssq S j c k X S j c j " S k c k X B j (20)
" H.

Since S k c k is non empty by Assumption 2, condition (20) ensures the existence of some

i k P rr1, N ss verifying i k R k´1 Ť j"0 B j contradicting k´1 Ť j"k´τ B j " rr1, N ss, as k ě τ .
Case 2: ι k ą 0. We have c ι k ´1 " c k and a finite induction leads to p@j P rrι k , kssq

S ι k c k " S ι k c ι k ´1 " S j c ι k ´1 " S j c k . (22)
Majoration ι k ď k ´τ implies that p@j P rrk ´τ, kssq

S ι k c k " S j c k . (23)
Non-overlap rule (5) with c k´τ , . . . , c k´1 ‰ c k then gives

p@j P rrk ´τ, k ´1ssq S j c k X S j c j " S ι k c k X B j (24)
" H.

Since S ι k c k is non empty, Condition (24) thus ensures the existence of

i k P rr1, N ss verifying i k R k´1 Ť j"k´τ B j which contradicts k´1 Ť j"k´τ B j " rr1, N ss. Assumption 3. (Curvature of majorizing matrix) (i) The mapping A : x P R N Þ Ñ Apxq P S N
``is such that (3) holds. Moreover, there exists ν ą 0 such that, for all S P T and v P R |S| ,

0 ă A pSq pvq ĺ νI |S| . (26)
(ii) There exists ν ą 0 such that, for all k P N,

Γ k c " A pB k q px ι k q ´1 2 A pB k q px k q ľ ˜L? τ p1 `τ q 2 `ν¸I |B k | . (27)
Assumption 3(i) is standard in optimization literature dealing with MM methods involving quadratic surrogates [START_REF] Chouzenoux | A majorize-minimize strategy for subspace optimization applied to image restoration[END_REF]. Assumption 3(ii) assumes that the spectrum of the difference between delayed and exact majorizing matrices of the partial quadratic majoring functions is strictly greater than a certain constant. This hypothesis controls the length of the MM increments performed by each worker. It aims at ensuring consistency between the asynchronous updates, by directly relating the worst-case curvature of the function f (parameterized by the Lipschitz constant L) and the worstcase communication delay (parameterized by the constant τ ). Condition [START_REF] Denis | Fast approximations of shiftvariant blur[END_REF] is key to ensure a condition descent for the general process generated by BD3MG scheme (see subsection 4.1). Assumption 3(ii) becomes redundant with Assumption 3(i) in the case when no delay occurs (i.e., τ " 0). A detailed constructive example on how to meet Assumption 3(ii) will be provided in our experimental Section 5.

Technical lemmas

We conclude this section by presenting some preliminary results that will be useful for our convergence analysis. Lemma 3.1. Under Assumption 2, for every k ě τ ,

› › ›x k ´xι k › › › 2 ď τ k ÿ j"k´τ `1 › › ›x j ´xj´1 › › › 2 . ( 28
)
Proof. Let k P N. If ι k " k, inequality ( 28) is trivial. For the rest of the proof we thus suppose ι k ď k ´1. According to the definition of the euclidean norm we have

› › ›x k ´xι k › › › 2 " N ÿ i"1 ´xk i ´xι k i ¯2 . (29)
Then, for all i P 1, N , the Jensen's inequality leads to

´xk i ´xι k i ¯2 " ¨k ÿ j"ι k `1 ´xj i ´xj´1 i ¯' 2 ď pk ´ιk q k ÿ j"ι k `1 ´xj i ´xj´1 i ¯2 . (30)
Moreover, Proposition 3.1 ensures that ι k belongs to k ´τ `1, k . As a consequence

`@i P 1, N ˘´x k i ´xι k i ¯2 ď τ k ÿ j"k´τ `1 ´xj i ´xj´1 i ¯2 . (31)
We then replace (31) in [START_REF] Escande | Sparse wavelet representations of spatially varying blurring operators[END_REF], which yields

› › ›x k ´xι k › › › 2 ď τ k ÿ j"k´τ `1 N ÿ i"1 ´xj i ´xj´1 i ¯2 . (32)
Relation [START_REF] Duval | Euclid in a taxicab: Sparse blind deconvolution with smoothed l1/l2 regularization[END_REF] directly comes from the identification of the inner sum of (32) as › › x j ´xj´1 › › 2 for all j P k ´τ `1, k .

Lemma 3.1 provides a bound on the residual between x k and the delayed vector x ι k updated by worker c k at iteration k P N. The right term in [START_REF] Duval | Euclid in a taxicab: Sparse blind deconvolution with smoothed l1/l2 regularization[END_REF] can be understood as the extra information available to the master, when compared to the one available to worker c k . This Lemma will allow to establish a descent condition on the BD3MG process in the next Section. Lemma 3.2. Under Assumptions 1 and 3(i), for every k P N,

}∇ pB k q f px ι k q} 2 ď ν 2 }x k`1 ´xk } 2 . ( 33
)
Proof. Let k P N. Let us analyse the quantity f px ι k q ´QpB k q px k`1 pB k q , x ι k q.

On the one hand, function

Ψ k : α P R Þ Ñ Q B k px k pB k q
´Dk αe, x ι k q with e " p1, 0q J is a second degree convex polynomial with a unique minimizer that reads

p α k " }∇ pB k q f px ι k q} 2 }∇ pB k q f px ι k q } 2 A pB k q px ι k q . ( 34
) Since u k is a minimizer of Q B k px ι k pB k q `Dk ., , x ι k q " Q B k px k pB k q
`Dk ., x ι k q, with x ι k pB k q " x k pB k q by Equation ( 11), we deduce that

Q B k px k`1 pB k q , x ι k q ď Ψ k pp α k q " f px ι k q ´1 2 p α k }∇ pB k q f px ι k q} 2 . ( 35
)
From Assumption 3(i), p α k verifies p α k ě ν ´1. Equation ( 35) can thus be rewritten as

f px ι k q ´QpB k q px k`1 pB k q , x ι k q ě 1 2ν }∇ pB k q f px ι k q} 2 . ( 36
)
On the other hand, using [START_REF] Chouzenoux | SABRINA: a stochastic subspace majorization-minimization algorithm[END_REF] from Definition (2), and Equation ( 16) yield

f px ι k q ´QpB k q px k`1 pB k q , x ι k q " A ∇ pB k q f px ι k q , x k`1 pB k q ´xι k pB k q E `1 2 › › ›x k`1 pB k q ´xι k pB k q › › › 2 A pB k q px ι k q " A ∇ pB k q f px ι k q , x k`1 pB k q ´xk pB k q E `1 2 › › ›x k`1 pB k q ´xk pB k q › › › 2 A pB k q px k q " 1 2 › › ›x k`1 pB k q ´xk pB k q › › › 2 A pB k q px ι k q . ( 37
)
The combination of ( 36) and (37) leads to

}∇ pB k q f px ι k q} 2 ď ν › › ›x k`1 pB k q ´xk pB k q › › › 2 A pB k q px ι k q . ( 38
)
Finally, Equation (33) comes using Assumption 3(i), and in particular,

› › ›x k`1 pB k q ´xk pB k q › › › 2 A pB k q px ι k q ď ν › › ›x k`1 pB k q ´xk pB k q › › › 2 (39)
" ν}x k`1 ´xk } 2 . ( 40
)
Lemma 3.2 generalizes the decreasing behavior of standard MM schemes [START_REF] Chouzenoux | A majorize-minimize subspace approach for ℓ 2 ´ℓ0 image regularization[END_REF][START_REF] Fest | Convergence analysis of block majorize-minimize subspace approaches[END_REF] to the asynchronous context. It is not directly invoked in our main convergence proof but serves as an intermediary to show the following technical result. Lemma 3.3. Under Assumptions 1 and 3(i), for all k ě 2τ ,

}∇f px k q} ď Lτ k ÿ j"k´2τ `1 }x j ´xj´1 } `ν k ÿ j"k´τ }x j`1 ´xj }. ( 41
) that is p1 ´rP q ř n k"k ˚uk ď rP ř k ˚´1 k"0 u k `ř`8 k"0 v k .
With 0 ă 1 ´rP ă 1, we deduce the summability of `uk ˘kPN .

Lemma 3.4 is a technical result to ensure the convergence of some real series. Several variants of inequality [START_REF] Kim | Blind depth-variant deconvolution of 3d data in wide-field fluorescence microscopy[END_REF] have been used to prove the finite length of iterative processes and then their convergence [START_REF] Davis | The asynchronous PALM algorithm for nonsmooth nonconvex problems[END_REF][START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF].

Convergence results

Let us now state our main theoretical results, that relate to the convergence properties of BD3MG iterates. Our proof line is reminiscent from [START_REF] Davis | The asynchronous PALM algorithm for nonsmooth nonconvex problems[END_REF][START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF] and follows similar steps that we summarize hereafter. First, starting from the majoration property (3) and using Lemma 3.1, we will establish a descent inequality. The latter is the key point of the rest of the proof. In particular, it will allow to show convergence of `f px k q ˘kPN . Then, Lemma 3.3 will ensure that `∇F px k q ˘kPN converges to 0 N , and usual topological properties will serve to show that the set of cluster points Cpx k q kPN q of px k q kPN lies in the set of stationary point of f . Finally, we will exhibit a Lyapunov function from our descent equality and will resort to the Kurdyka-Lojasewicz (KL) inequality [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF] to prove our main theorem, showing the convergence of the BDM3G iterates and providing a rate of convergence.

Descent inequality

Proposition 4.1. Under Assumptions 1-2-3, there exists a positive sequence pν k q kěτ such that p@k ě τ q f px k`1 q `νk`1 ď f px k q `νk ´ν}x k`1 ´xk } 2 .

Proof. By definition of the majorant function (3), for every k P N,

f px k`1 q ď f px k q `A∇ pB k q f px k q, x k`1 pB k q ´xk pB k q E `1 2 }x k`1 pB k q ´xk pB k q } 2 A pB k q px k q . ( 50
)
Decomposing the scalar product term then yields, for every k P N,

f px k`1 q ď f px k q `Rk `A∇ pB k q f px ι k q, x k`1 pB k q ´xk pB k q E `1 2 }x k`1 pB k q ´xk pB k q } 2 A pB k q px k q , (51)
with R k "

A ∇ pB k q f px k q ´∇pB k q f px ι k q, x k`1 pB k q ´xk pB k q E .
Let τ defined as in Assumption 2. A majoration of R k for every k ě τ comes by using successively Cauchy-Schwartz inequality, L gradient-Lipschitz inequality from Assumption 1, and Lemma 3.1:

p@k ě τ q R k ď L}x k ´xι k } }x k`1 pB k q ´xk pB k q } ď L 2 ? τ }x k ´xι k } 2 `L? τ 2 }x k`1 pB k q ´xk pB k q } 2 , ď L ? τ 2 k ÿ j"k´τ `1 }x j ´xj´1 } 2 `L? τ 2 }x k`1 pB k q ´xk pB k q } 2 . ( 52
)
We then set, for all k ě τ , ν k " L ?

τ 2 k ř j"k´τ `1pj ´k `τ q}x j ´xj´1 } 2 . Since x k`1 pB k q ´xk pB k q " x k`1 ´xk , (52) also reads p@k ě τ q R k ď ν k ´νk`1 `Lτ ? τ 2 }x k`1 pB k q ´xk pB k q } 2 `L? τ 2 }x k`1 pB k q ´xk pB k q } 2 , " ν k ´νk`1 `L? τ p1 `τ q 2 }x k`1 pB k q ´xk pB k q } 2 . ( 53
)
Moreover, Equation ( 16) ensures that p@k ě τ q

A ∇ pB k q f px ι k q, x k`1 pB k q ´xk pB k q E " ´}x k`1 pB k q ´xk pB k q } 2 A pB k q px ι k q . (54)
Replacing both ( 53) and ( 54) in [START_REF] Loizou | Revisiting randomized gossip algorithms: General framework, convergence rates and novel block and accelerated protocols[END_REF] gives, for all k ě τ , f px k`1 q `νk`1 ď f px k q `νk `L? τ p1 `τ q 2 }x k`1

pB k q ´xk pB k q } 2 ´}x k`1 pB k q ´xk pB k q } 2 Γ k c " f px k q `νk ´}x k`1 pB k q ´xk pB k q } 2 Γ k c ´L? τ p1`τ q 2 I |B k | , (55)
with Γ k c defined in Assumption 3(ii). Inequality ( 49) is a direct consequence of Assumption 3(ii) remarking that }x k`1 pB k q ´xk pB k q } " }x k`1 ´xk }.

General behaviour

We now state our first convergence Theorem for BD3MG algorithm.

Theorem 4.1. Under Assumptions 1-2-3, sequence `f px k q ˘kPN generated by BD3MG converges to a finite limit f 8 . Moreover, `∇f px k q ˘kPN converges to 0 N .

Proof. Coercivity of f (Assumption 1) and the descent inequality (49) guarantee that `f px k q `νk ˘kPN is a decreasing and lower-bounded sequence. It is thus converging to a real value f 8 . Equation ( 49) then directly leads to ř `8 k"0 }x k`1 ´xk } 2 ă `8. On the first hand, using the same notation pν k q kPN introduced in our proof of Proposition 4.1, we have

p@k ě τ q ν k ď L ? τ 2 `8 ÿ j"k`τ ´1 }x j ´xj´1 } 2 . (56)
Thus, the sequence pν k q kPN converges to 0 and so, by Proposition 4.1, `f px k q ˘kPN goes to f 8 . On the other hand, Lemma 3.3 gives p@k ě 2τ q }∇f px k q} ď Lτ `8 ÿ j"k´2τ `1 }x j ´xj´1 } ``8 ÿ j"k´τ }x j`1 ´xj },

which enables to conclude that `∇f px k q ˘kPN converges to 0 N . Proposition 4.2. Under Assumptions 1-2-3, C `px k q kPN ˘, defined as the set of accumulation points of `xk ˘kPN , is non empty , compact and is a subset of the set of stationary points of f . Moreover, f takes value f 8 on C `px k q kPN ˘.

Proof. Coercivity of f (by Assumption 1) and convergence of `f px k q ˘kPN to f 8 (by Theorem 4.1) show that `xk ˘kPN is a bounded sequence and C `px k q kPN ˘is non empty and compact. Convergence of `∇f px k q ˘kPN to 0 N (by Theorem 4.1) guarantees that every point x ˚P C `px k q kPN ˘is a stationary point of f . Moreover, using again `f px k q ˘kPN converging to f 8 yields f 8 " f px ˚q for every x ˚P C `px k q kPN ˘which concludes the proof.

Lyapunov-based asymptotical analysis

In order to establish the convergence of the iterates of BD3MG, we will follow an analysis relying on the use of a Lyapunov function. Such proof technique has also been used in [START_REF] Davis | The asynchronous PALM algorithm for nonsmooth nonconvex problems[END_REF][START_REF] Zavriev | Heavy-ball method in nonconvex optimization problems[END_REF][START_REF] Wilson | Lyapunov Arguments in Optimization[END_REF]. The idea is to exhibit a function, related to the loss function f but non necessarily equals to it, that decreases monotonically along the iterative process. Given Inequality [START_REF] Li | Surrogate-based distributed optimisation for expensive black-box functions[END_REF], a natural choice is

L : Z " ¨Z0 . . . Z τ ‹ ‹ ' P R pτ `1qN Þ Ñ f pZ 0 q `L? τ 2 τ ÿ ℓ"1 pτ ´ℓ `1q}Z ℓ ´Zℓ´1 } 2 . ( 58
)
Let us denote, for every k ě τ , Z k " ¨xk . . .

x k´τ ‹ ‹ ' P R pτ `1qN , with x k the k-th BD3MG
iterate. Then, the descent condition from Proposition 4.1 can be rewritten as

p@k ě τ q LpZ k`1 q ď LpZ k q ´ν}x k`1 ´xk } 2 . ( 59
)
The structure of L allows to build an upper bound of its gradient norm along the iterates, where the bound depends only on the differences of the past iterates: Lemma 4.1. There exists ρ ą 0 such that

p@k ě τ q }∇LpZ k q} 2 ď ρ k ÿ j"k´τ `1 }x j ´xj´1 }. ( 60
)
' φp0q " 0.

' φ P C 1 pp0, ζqq and is continuous in 0.

' @s P p0, ζq, φ 1 psq ą 0.

We are then ready to introduce the so-called KL property. [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] Definition 4.1. [KL property] A differentiable function g : R d Ñ R, with d ě 1, satisfies the Kurdyka-Lojasiewicz (KL) property on E Ă R d if, for every r

x P E and every bounded neighborhood V of r

x, there exist ζ ą 0 and φ P Φ ζ such that every x P E X tx s.t. |gpxq ´gp r xq| ă ζu,

}∇gpxq}φ 1 p|gpxq ´gpr xq|q ě 1. ( 67
)
We also recall the following Lemma:

Proposition 4.3. Under Assumptions 1-2-3, if L defined in (58) fulfills the KL property on R pτ `1qN then, considering g " L, C " C `pZ k q kPN ˘with LpCq " tf 8 u, there exists ϵ L , ζ L and ϕ L P Φ ζ L such that L satisfies [START_REF] Tseng | Convergence of a block coordinate descent method for nondifferentiable minimization[END_REF].

Proof. This is a direct consequence of Lemma 4.2. Continuity of L is clear. We still have to verify the compactness of C `pZ k q kPN ˘and that L is constant valued on that set. C `pZ k q kPN ˘is closed. Moreover, it is straightforward to show that this set is included in the Cartesian product

" C `px k q kPN ˘ıτ`1
, where C `px k q kPN ˘is compact. C `pZ k q kPN ȋs thus bounded and, finally, it is compact.

Let Z P C `pZ k q kPN ˘. We have LpZ k q " f px k q `νk for all k P N. From our proof of Theorem 4.1, it follows that sequence `LpZ k q ˘kPN converges to f 8 . Continuity of L finally ensures that f 8 " LpZq. This proves that f is constant valued on C `pZ k q kPN (and equals to f 8 ).

Convergence of the iterates

We are now ready to state our second convergence Theorem for BD3MG algorithm, characterizing the convergence of px k q kPN . Theorem 4.2. Let assume that Assumptions 1-2-3 hold. Assume furthermore that the Lyapunov function L in (58) satisfies the KL property on R pτ `1qN . Then, sequence px k q kPN is of finite length, i.e :

`8 ÿ k"0 }x k`1 ´xk } ă `8, (69)
and converges to a stationary point of f .

Proof. Let us start considering the case when there exists some k 0 P N where LpZ k 0 q " f 8 . Since `LpZ k q ˘kPN is decreasing sequence converging to f 8 (see proof of Proposition 4.3), it follows that LpZ k q " f 8 for all k ě k 0 . Descent inequality (59) then gives p@k ě k 0 q }x k`1 ´xk } 2 ď ν ´1 ´LpZ k q ´LpZ k`1 q ¯" 0, [START_REF] Wang | A field guide to federated optimization[END_REF] ensuring that px k q kPN has a finite length and x k , k ě 0, is a stationary point of f .

We now suppose that, for all k P N, LpZ k 0 q ‰ f 8 . We aim at exhibiting a uniform KL inequality on sequence `LpZ k q ˘kPN . To do so, let us peruse the quantities ϵ L , η L , φ L arising from Proposition 4.3. On the one hand, the decrease of `LpZ k q ˘kPN implies that, for all k P N, LpZ k q ą f 8 . The set C `px k q kPN ˘is non empty (see proof of Proposition 4.2), so is the set C `pZ k q kPN ˘.

Let Z P C `pZ k q kPN ˘an element of such set i.e., a cluster point of pZ k q kPN . From Proposition 4.3, LpZq " f 8 . Hence, LpZ k q ´LpZq ą 0 for all k P N.

On the other hand, `LpZ k q ˘kPN converges to f 8 " LpZq. The boundedness of px k q kPN also ensures this of pZ k q kPN . We deduce the existence of some k 1 ě 2τ such that

p@k ě k 1 q 0 ă LpZ k q ´LpZq ă η L , d ˆZk , C ´pZ k q kPN ¯˙ă ϵ L . (71)
From Proposition 4.3, the uniform KL property on L holds i.e., p@k ě k 1 q }∇LpZ k q} ´φL ¯1 ´LpZ k q ´LpZq ¯ě 1.

Moreover, setting ∆ k " φ L `LpZ k q ´LpZq ˘´φ L `LpZ k`1 q ´LpZq ˘for all k P N, concavity of φ L and the descent inequality (59) ensure that

p@k ě k 1 q ∆ k ě ´φL ¯1 ´LpZ k q ´LpZq ¯´LpZ k q ´LpZ k`1 q ě ν}x k`1 ´xk } 2 ´φL ¯1 ´LpZ k q ´LpZq ¯. ( 73
)
The combination of the latter with (72) leads to

p@k ě k 1 q }x k`1 ´xk } 2 ď ν ´1∆ k }∇LpZ k q}. ( 74
)
By Lemma 4.1, we can upper bound the gradient term in [START_REF] Zavriev | Heavy-ball method in nonconvex optimization problems[END_REF]. This gives

p@k ě k 1 q }x k`1 ´xk } 2 ď ρν ´1∆ k k ÿ j"k´2τ `1 }x j ´xj´1 }. ( 75
)
Passing to the root and using the classical identity ? ab ď a{c `bc{4, with a " k ř j"k´2τ `1 }x j ´xj´1 } for all k ě k 1 , b " ∆ k , both positive for all k ě k 1 and some c ą 0 is generic, leads to

p@k ě k 1 q }x k`1 ´xk } ď ? ρν ´1{2 c k ÿ j"k´2τ `1 }x j ´xj´1 } `c? ρν ´1{2 4 ∆ k . ( 76
)
Since p∆ k q kPN is summable (as a telescopic sequence), we can apply Lemma 3.4 with some c ą 2τ ? ρν ´1{2 so that 2τ ? ρν ´1{2 c ă 1. This shows that sequence px k q kPN has a finite length. This finite length property entails that px k q kPN is a Cauchy sequence and thus a converging one. The final conclusion directly comes from Proposition 3.1, ensuring that every accumulation point of px k q kPN is a stationary point of f .

Discussion

Under the KL condition for the Lyapunov function L defined in [START_REF] Onose | Scalable splitting algorithms for big-data interferometric imaging in the SKA era[END_REF], we were able to demonstrate the convergence of sequence px k q kPN to a stationary point of f . Let us notice that f satisfying KL property does not necessary imply that L does. Still, our assumption on L can be verified in practice for a wide class of functions f . For instance, following the discussion in [24, section 6], if f is semi-algebraic [START_REF] Bochnak | Real Algebraic Geometry[END_REF][START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF], then the required condition on L in Theorem 4.2 is satisfied, with function φ L " κp.q 1´θ for a some pκ, θq P R ˚ˆp0, 1q. Such situation will be met in our experimental settings in Section 5. Extending Theorem 4.2 to any KL function f would be an interesting avenue for future work but up to our knowledge, it does not seem straightforward.

Application to 3D image restoration

5.1. Problem statement 5.1.1. Observation model. We focus on the inverse problem of restoring a vectorized 3D volume x of size N " N X ˆNY ˆNZ given blurry and noisy observation y P R N . We consider a depth-variant blur operator H P R N ˆN and additive i.i.d. Gaussian noise with standard deviation σ ą 0, so that the observed volume is related to x through,

y " Hx `b, (77)
with vector b P R N accounting for the noise. The goal is to solve the inverse problem of estimating x given y and H. Depth-variant blurs are commonly encountered in 3D microscopy [START_REF] Preza | Depth-variant maximum-likelihood restoration for threedimensional fluorescence microscopy[END_REF][START_REF] Hadj | Modeling and removing depth variant blur in 3d fluorescence microscopy[END_REF][START_REF] Kim | Blind depth-variant deconvolution of 3d data in wide-field fluorescence microscopy[END_REF][START_REF] Jezierska | Spatially variant psf modeling in confocal macroscopy[END_REF], due to optical aberrations. They are particular cases of spatially-variant blurs [START_REF] Chakrabarti | Analyzing spatially-varying blur[END_REF][START_REF] Nagy | Restoring images degraded by spatially variant blur[END_REF]. The degradation operator H raises specific challenges due to its high computational cost. Several strategies have been investigated in the case of 2D spatially variant blur maps encountered for instance in astronomical imaging [START_REF] Denis | Fast model of space-variant blurring and its application to deconvolution in astronomy[END_REF][START_REF] Denis | Fast approximations of shiftvariant blur[END_REF][START_REF] Escande | Sparse wavelet representations of spatially varying blurring operators[END_REF]. The extension to 3D maps of these methods is however not covered up to our knowledge. This motivates the use of a distributed optimization approach for solving the inverse problem (77).

Objective function

We adopt a variational strategy, which consists in seeking for an estimate of x that minimizes a penalized least squares criterion f . A hybrid regularization term is employed incorporating prior knowledge on the smoothness and the range of the sought solution. The objective function reads:

p@x P R N q f pxq " S ÿ s"1 f s pL s xq, (78)
where, for every s P t1, . . . , Su, L s P R PsˆN , P s P N ˚, and f s is a function from R Ps to R. f 1 ˝L1 represents the data fidelity term while the other terms are regularization terms. Here, we set S " 4 and

' P 1 " N , L 1 " H, f 1 " 1 2 } ¨´y} 2 , ' P 2 " N , L 2 " I N , f 2 " η d 2 rx min ,xmaxs N , ' P 3 " 2N , L 3 " rpV X q J pV Y q J s J , f 3 " λ N ř n"1 b r.s 2 n `r.s 2 N `n `δ2 , ' P 4 " N , L 4 " V Z , f 4 " κ} ¨}2 .
Hereabove, pη, λ, δ, κq P p0, `8q 4 are hyper-parameters. The linear operators V X , V Y , V Z P R N ˆN are discrete gradient operators along X (horizontal), Y (vertical), and Z (longitudinal) directions of the 3D volume. Function d 2 rx min ,xmaxs N states for the squared distance to set rx min , x max s N Ă R N , with px min , x max q P R 2 minimal and maximal bounds on the sought intensity values. The later term can be viewed as an exterior penalty function [START_REF] Chouzenoux | A local MM subspace method for solving constrained variational problems in image recovery[END_REF]. ;.

Majorant mapping.

In order to implement BD3MG, we must build a majorant mapping ensuring the majorization condition (3). Function f fits within the class of half-quadratic majorizing constructions initially introduced in [START_REF] Geman | Constrained restoration and the recovery of discontinuities[END_REF][START_REF] Geman | Nonlinear image recovery with half-quadratic regularization[END_REF] and later analysed in [START_REF] Nikolova | Analysis of half-quadratic minimization methods for signal and image recovery[END_REF][START_REF] Allain | On global and local convergence of half-quadratic algorithms[END_REF][START_REF] Chouzenoux | A majorize-minimize strategy for subspace optimization applied to image restoration[END_REF]. A general structure for the majorant mapping of (78) is

p@x P R N q Apxq " S ÿ s"1 pL s q J Diag 1ďpďPs ␣ rω s pL s xqs p ( L s , (79)
where, for every s P t1, . . . , Su, ω s : R Ps Ñs0, `8r Ps is a majorizing potential that depends on the properties of pf s q 1ďsďS [20, Tab. I]. In our case, for every s P t1, . . . , 4u, each of these terms is

f s is β s -Lipschitz differentiable with $ ' ' ' ' ' & ' ' ' ' ' % β 1 " 1, β 2 " 4η, β 3 " λδ ´1, β 4 " 2κ. (80)
;

Function x P R N Þ Ñ d 2 E
pxq is 2-Lipschitz differentiable a soon as E is non-empty closed and convex set [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF]. Denoting by p E the orthogonal projection operator, its gradient then corresponds to

x P R N Þ Ñ 2px ´PE pxqq.
Then, from descent lemma [START_REF] Bertsekas | Nonlinear programming[END_REF], a valid choice is ω s p¨q " αβ s 1 Ps with some α ě 1 [START_REF] Chouzenoux | A majorize-minimize strategy for subspace optimization applied to image restoration[END_REF]. We adopt this simple strategy for functions f 1 , f 2 and f 4 , which yields

$ ' ' & ' ' % ω 1 p¨q " α1 N ω 2 p¨q " 4αη1 N ω 4 p¨q " 2ακ1 N . (81)
Regarding function f 3 , a more sophisticated majorization is adopted, inherited from half-quadratic strategies [START_REF] Allain | On global and local convergence of half-quadratic algorithms[END_REF][START_REF] Nikolova | Analysis of half-quadratic minimization methods for signal and image recovery[END_REF]:

p@v P R 2N q ω 3 pvq " λ » - - - - ˆ1{ b v 2 n `v2 N `n `δ2 ˙1ďnďN ˆ1{ b v 2 n `v2 N `n `δ2 ˙1ďnďN fi ffi ffi ffi fl . (82)
5.1.4. Distributed implementation We implement BD3MG algorithm as presented in Section 2.4. We split the 3D volume into 2D slices along the depth axis z P t1, . . . , N Z u, and consider each 2D slice as an individual block upon which workers can compute an update. Assuming a lexicographic ordering of the voxels, this means that the following partition is adopted:

T " ␣ pi ´1qN X N Y `1, iN X N Y | 1 ď i ď N Z ( . (83)
BD3MG is implemented on a star graph of workers with a specific master aggregating the current solution. For a given number of active cores C tot " C `1 of the computer (or of the cluster), one is used as the master process to manage the computation split between the workers on a first-come, first-served basis (i.e., Algorithm 1), while all the Cp" C tot ´1q others, are computing updates asynchronously on planar blocks (i.e., Algorithm 2).

Validity of Assumptions.

Let us discuss the validity of Assumptions 1, 2 and 3 for the considered problem and implementation.

Assumption 1. Function f in (78) is differentiable. Moreover, it has a L-Lipschitzian gradient with L " ř S s"1 β s |||L s ||| 2 , where ||| ¨||| denotes the spectral norm over matrices and pβ s q 1ďsďS were given in the previous subsection. According to [64, Prop. 2.5], a sufficient condition for f to be coercive is KerpHq " t0 N u. This latter is verified in our experiments, since H is a full-rank operator. Thus, Assumption 1 holds.

Assumption 2. This assumption relates to the practical implementation of BD3MG. It requires that every feedback from a worker happens within a finite number of iterations. This can be upper bounded by N o , the number of necessary operations to compute one iteration of BP3MG scheme (i.e the pessimistic situation for BD3MG where one single worker has to focus on all the coordinates of 1, N ). As a result, at least one feedback is produced on every period T o , the necessary machine time to perform the N o operations (supposing every workers as the same performances). Assumption 2 is thus valid as soon as every blocks of the partition (83) is taken into consideration by one of the workers following a regularly time interval. In our practical implementation, we spread the N Z 2D slices regularly among the C workers, so that each worker c P t1, . . . , Cu is on average in charge of N X N Y N Z {C coordinates at every stage of the process. The master loop is then defined in such a way that each worker updates the 2D slices it has in charge, in a cyclic manner. This allows to control a clock time T c that guarantees every block S P T to be updated at least once every T c iteration. Then, Assumption 2 holds with τ " tT c {T o u `1.

Assumption 3. This assumption relates to the majorant mapping. To prove this assumption, we proceed in three steps. On the one hand, we have,

p@x P R N q Apxq ľ pL 2 q J Diag 1ďpďP 2 ␣ ω 2 pL 2 xq ( L 2 ľ αηI N . (84)
On the other hand, according to definition (79) and those of ω 1 , . . . , ω 4 p@x P R N q Apxq ĺ ¨S ÿ

s"1

|||L s ||| 2 max 1ďpďPs rω s pL s xqs p 'I N ĺ νI N , (85) with
ν " α

¨S ÿ

s"1

β s |||L s ||| 2 '. (86)
Considering ( 84), ( 85)-( 86) and the fact that any sub-matrix M pSq pS Ă 1, N ) of a (symmetric) positive matrix M remains positive, the chosen mapping A thus respects conditions imposed by Assumption 3(i). Moreover, for all px, yq P pR N q 2 , Apxq ´1 2 Apyq

" S ÿ s"1 pL s q J Diag 1ďpďPs # ˆrω s pL s xqs p ´1 2 rω s pL s yqs p ˙+ L s (87) " α 2 ÿ sPt1,2,4u pL s q J L s `pL 3 q J Diag 1ďpďP 3 # ˆrω 3 pL 3 xqs p ´1 2 rω 3 pL 3 yqs p ˙+ L 3 ľ α 2 L J 2 L 2 `pL 3 q J Diag 1ďpďP 3 # ˆrω 3 pL 3 xqs p ´1 2 rω 3 pL 3 yqs p ˙+ L 3 ľ ηpαqI N with ηpαq " α 2 ´λ 2δ |||L 3 ||| 2 . ( 88
)
and, under the same previous remark on the block positivity preservation, Assumption 3(ii) is verified considering α large enough (i.e so as for ηpαq to strictly exceed bound

L ? τ p1`τ q 2
Convergence result In a nutshell, Assumptions 1-2-3 are fulfilled in our experiments, so that Theorem 4.1 holds. Moreover, function f is semi-algebraic, hence so is the Lyapunov function L (see discussion in Sec. 4.5). Thus, Theorem 4.2 holds.

Comparative analysis on a controlled scenario

We first set x as the 3D microscopic image FlyBrain § with size N " 256ˆ256ˆ57. The linear operator H models a 3D depth-varying Gaussian blur, with kernel size 5 ˆ5 ˆ11. For each depth z P t1, . . . , N Z u, the blur kernel is characterized by different variance and rotation parameters pσ X pzq, σ Y pzq, σ Z pzq, φ Y pzq, φ Z pzqq, following the model from [START_REF] Wirjadi | Approximate separable 3D anisotropic Gauss filter[END_REF]. In practice, the values of these five parameters are chosen randomly through a uniform distribution over r0, 3s ˆr0, 3s ˆr0, 4s ˆr0, 2πs ˆr0, 2πs, sampled independently for every z. A zero-mean white Gaussian noise with standard deviation σ " 4 ˆ10 ´2 is then added to the blurred volume. The regularization parameters pλ, δ, κ, ηq " p1, 1, 10 ´1, 10 ´3q are chosen empirically so as to maximize the Signal-to-Noise Ratio (SNR) of the restored volume. Moreover, we set px min , x max q " p0, 1q, equal to the range of the ground truth image. In order to illustrate the acceleration induced by the proposed BD3MG, we run a comparative analysis between different versions of the optimization scheme, in the spirit of an ablation study. Namely, we propose to compare BD3MG with three methods listed hereafter.

' The 3MG algorithm [START_REF] Chouzenoux | A majorize-minimize strategy for subspace optimization applied to image restoration[END_REF][START_REF] Chouzenoux | A majorize-minimize subspace approach for ℓ 2 ´ℓ0 image regularization[END_REF] is considered as the baseline. At each iteration, this algorithm builds the majorant mapping as in Sec. 5.1.3 and computes memory gradient updates on the full volume, without any parallelization.

' The Asynchronous Block Gradient Descent (ABGD) algorithm implements the method from [START_REF] Niu | Hogwild: A lockfree approach to parallelizing stochastic gradient descent[END_REF]. It performs parallel asynchronous gradient descent updates over the slices of the volume. We adopt here the same parallelization settings as for our BD3MG. Updates correspond to the standard gradient descent on the selected planar blocks, using a fixed step-size µ ensuring convergence of the iterative scheme, namely µ " 0.99{p1 `κ `2λ{δ `2κq.

' The BP3MG algorithm from [START_REF] Cadoni | A block parallel majorize-minimize memory gradient algorithm[END_REF][START_REF] Fest | Convergence analysis of block majorize-minimize subspace approaches[END_REF] runs a synchronous version of BD3MG algorithm. The master process carries out the main loop of [START_REF] Cadoni | A block parallel majorize-minimize memory gradient algorithm[END_REF]Alg 4.3]. At each iteration k P N, it selects C block indices (following a cyclic rule) and sends to each worker c P 1, C the required data allowing it to update S k c , the associated block. Workers wait for each other to finish their tasks, combine their respective updates into a unique vector px j q jPS k 1 Y...YS k C and finally send the latter to the master. This approach could be interpreted as a special case of BD3MG with a single worker (potentially composed of several subworkers) working on sending its update to the central process (potentially composed of several sub-updates) S k " tS k c u cPC . Since there is no mismatch in information between central process and workers in this synchronous version, the delay vector i k always equals k. § https://imagej.nih.gov/ij/images/ All methods are implemented in Python using the built-in Multiprocessing library as well as Numpy and Scipy for both data manipulation and scientific computing. The experiments of this section are conducted on an Intel® Xeon(R) W-2135 CPU with 12 cores clocked at 3.70GHz. All the versions were initialized with x 0 " 0 N leading to an initial value f px 0 q " 91292.92. During computation and at every iteration k P N ˚, we monitor the cost function f px k q, the normalized increment }x k`1 ´xk }{}x k }, and the signal to noise ratio (SNR, in dB) defined as SNR " 20 log 10 ˆ}x} }x ´xk } ˙.

(89)

We set a stopping criterion to }x k`1 ´xk } ď ε }x k } with ε a small threshold. The obtained solution is denoted as x f . In this example, when using BD3MG, a threshold of ε " 10 ´3 yields a final SNR of 18.33 dB and f px f q " 1261.59. This is on par with solutions found for the other tested algorithms. BD3MG exhibits a faster practical convergence (See Table 1 and Figure 2 for comparison). Slices of the reconstructed volume are displayed in Figure 3, revealing fine details of the image recovered by the restoration procedure.

Effect of an imbalanced computing power

In order to further demonstrate the advantages of BD3MG over its synchronous counterpart BP3MG, we tested the methods under different computing environments by synthetically modeling stochastic delays in the computing loop of workers. More specifically, the same restoration task and computer characteristics than in the previous section is considered, and we introduce artificial perturbation in the computing environment by randomly "freezing" some worker processes for a certain amount of time (i.e., delay) following the three scenarios below: ' Type I: One of the workers is consistently affected by a delay that follows a uniform distribution Upr0, 1sq (in sec.). The other cores are not affected by any delay.

' Type II: Two worker cores are not affected by any delay while the others 9 agents are delayed in the following fashion:

3 cores hold a delay following a uniform distribution Upr0, 1sq. The results are summarized in Table 2 and Figure 4. In all three scenarios, BD3MG outperforms its synchronous version BP3MG, in terms of computation time while reaching similar final criterion value and SNR. The criteria decrease is faster for BD3MG due to accelerated iterations of the workers. A more efficient handling of the workload is performed, as shown in Figure 4 where CPU idle time is consistently lower for BD3MG than for BP3MG. We note that in asymmetric settings such as (Type II) and (Type III), BD3MG proved to be particularly efficient in reducing the synchronicity constraint of BP3MG for "fast" workers. The comparable results for BD3MG on all three scenarios further suggest that the proposed algorithm is more robust to an imbalance in the computing power of workers.

Scalability assessment.

In order to assess the scalability properties of BD3MG, we further analyse the speed-up generated by the number of cores available. We consider the restoration problem of the 3D image Aneurysm} of size N " 256 ˆ256 ˆ154, under the same } http ://www.mathworks.com/matlabcentral/fileexchange /25987-showvol-isosurface-render/ degradation operator and noise level than in the previous example. Figure 5 presents the acceleration ratio between the required computation time for one single worker versus the computation time of up to 30 active workers in reaching the stopping criterion ε " 10 ´3. The regularization parameters are set empirically to pλ, δ, κ, ηq " p1, 1, 10 ´1, 10 ´3q to maximize the final SNR. The computations were performed using HPC resources from the Oscar -Ocean State Center for Advanced Resources of the Center for Computation and Visualization, Brown University. The hardware is an Intel Corei9 CPU with up to 48 physical cores at 3.3 GHz GHz and 300G of RAM. Results found in Figure 5 illustrate the great potential of scalability of the proposed algorithm. As the number of core increases, a mild saturation effect is observed (in agreement with Amdahl's law [START_REF] Patterson | Computer architecture: a quantitative approach[END_REF]).

Application to real data from multiphoton microscopy

We finally illustrate the performance of BD3MG on a restoration task of real multiphoton microscopy data specifically acquired for this experiment. Multiphoton microscopy is an interesting solution for the 3D and submicrometric characterization of biomedical structures, it is label-free and contactless [START_REF] Göbel | Imaging cellular network dynamics in three dimensions using fast 3d laser scanning[END_REF]. Such a solution takes advantage of optical sectioning, an optical property resulting from the nonlinear optical processes involved. 3D images are produced with sub-micrometer resolution without slicing the sample. We use an instrumental solution for acquisition rested on a commercial system from Olympus (BX61WI) coupled with a multiphoton water immersion objective (Olympus XLPLN25XWMP, 25×, NA 1.05). A laser system, emitting femtosecond pulses centred at 810 nm with 10 nm of spectral bandwidth, is used for production of the nonlinear phenomena of second harmonic generation (SHG) and two-photon fluorescence (TPF). The biomedical sample is made of a whole mouse muscle, the Extensor digitorum longus (EDL), isolated from tendon to tendon. Sub- micrometric fluorescent microspheres emitting in the green range are included into the EDL and spread homogeneously all along the whole muscle structure. Under such an experimental protocol, the production of two 3D images is obtained. The first channel contains the SHG from the myosin of the muscle and the second channel displayed the TPF of microspheres used for calibrating the instrumental PSF. A hundred of 2D image slices of SHG and TPF are produced, with 0.1 µm resolution along depth axis Z and 0.049µmˆ0.049µm resolution over X ´Y horizontal-vertical axis. The acquisition recording starts 140 µm under the sample surface for a total sample thickness of 180 µm. For this range of depth, the imaging of biological samples is degraded by scattering effects. Both raw volumes (i.e., SHG and TPF) dimension have 2048ˆ2048ˆ100 voxels, from which we extract a subpart with size 256 ˆ256 ˆ100 voxels for the purpose of our study.

We follow the computational pipeline FAMOUS previously introduced in [START_REF] Lefort | FAMOUS: a fast instrumental and computational pipeline for multiphoton microscopy applied to 3d imaging of muscle ultrastructure[END_REF]. We estimate a depth-variant Gaussian PSF field within the 3D microscopic volume by applying the 3D Gaussian fitting algorithm FIGARO from [START_REF] Chouzenoux | Optimal multivariate Gaussian fitting with applications to PSF modeling in two-photon microscopy imaging[END_REF] to volume of interests extracted from the second image channel, displaying fluorescence of calibrated microbeads. Each volume of interest is selected through an automatic search of connected components within a filtered and binarized version of the observed volume. Then, FIGARO method is ran, yielding parameters (i.e., mean, covariance, scaling, shift) of a 3D Gaussian shape. This allows to build, through a simple interpolation strategy, a model for a depth-variant PSF (see more details in [START_REF] Lefort | FAMOUS: a fast instrumental and computational pipeline for multiphoton microscopy applied to 3d imaging of muscle ultrastructure[END_REF]Sec.2.4]). Since no ground truth is available, the regularization parameters pλ, δ, κ, ηq " p10 2 , 2, 10, 10 ´3q, are selected by retrospective visual inspection. The reconstruction shown in Figure 6 exhibits clear contrasts and sharpness properties. Comparative videos of the original and restored volume are available at https://github.com/mathieuchal/BD3MG. The native signal from the raw image was presenting a high level of noise and blur due to the presence of scattering elements all along the 140 µm of sample depth. Thanks to the proposed restoration strategy, the localisation of the myosin in the muscle sample is made possible, and the spatial organization of this protein into the down side of the EDL is revealed. The volume restoration took 305 seconds and " 2000 iterations on a 12 cores setting, when using stopping criterion ε " 10 ´3. This represents a much faster convergence than the competing approaches, and, in particular, greatly increases computational tractability of manual hyperparameters finetuning.

Conclusion

In this paper, we have presented a new block distributed Majorize-Minimize algorithm, BD3MG, devised to tackle large-size differentiable optimization problems met in a wide range of applications. Our main contribution lies in a distributed asynchronous formulation that allows for delays in the current solution computed between workers, while securing convergence guarantees under mild assumptions. Our new algorithm BD3MG has been tested in the context of 3D image restoration with depth-variant blur.

Experimental results underlined the speedup potential of this method and its concrete applicability in the field of fluorescence microscopy. Future work will be dedicated to extension to more general distributed graph topologies.

Declarations

Figure 1 :

 1 Figure 1: Examples of graph topologies. The graph in (c) is encompassed by our framework.

) Algorithm 2 .

 2 BD3MG algorithm -Worker loop (Wait for a feedback from the master) (a) Receive px, S, dq from Master. (b) Dpx, S, dq " r´∇ pSq f pxq | ds. (c) Compute ∇ pSq f pxq and A pSq pxq.(d) u " ´`Dpx, S, dq J A pSq pxq Dpx, S, dq ˘: Dpx, S, dq ∇ pSq f pxq.

Case 1 :

 1 ι k " 0. By definition, c 0 , . . . , c k´1 ‰ c k and an easy induction gives S 0 c k " . . . " S k c k . Non-overlap rule (5) with c 0 , . . . , c k´1 ‰ c k yields

Lemma 4 . 2 .

 42 [Uniform KL property [8, Lemma 6]] Let C a compact set of R d and g : R d Ñ R a continuous function satisfying KL property on C and constant on the latter. Then, there exist ϵ, ζ ą 0 and φ P Φ ζ such that every x P C and all x P R d satisfying both dpx, Cq ă ϵ, 0 ă gpxq ´gpxq ă ζ, we have }∇gpxq}φ 1 p|gpxq ´gpxq|q ě 1.

Figure 2 :

 2 Figure 2: Evolution of quantitative metrics along time (in seconds), for algorithms 3MG (blue), ABGD (orange), BP3MG (green) and BD3MG (red), for FlyBrain restoration. Evolution of }x k ´x} (left), }x k`1 ´xk }{}x k } (middle), and SNR in dB (right).

Figure 3 :

 3 Figure 3: Restoration results of Flybrain: ground truth volume (top), degraded version (middle), and results of BD3MG restoration (bottom). Visual comparison along the X ´Z axis (left) the X ´Y axis (middle) and zoomed details (right). The optimization process recovers fine details of the original volume that were lost in its degraded version.

Figure 4 :

 4 Figure 4: Numerical comparisons between BD3MG and BP3MG for FlyBrain restoration under imbalanced computing power: evolution of }x k`1 ´xk }{}x k } along time (in sec.) for each of the three experimental settings in log-log scale (left), and averaged ratio of workers CPU idle time over the entire optimization process for each scenario (right).

Figure 5 :

 5 Figure 5: Speed-up ratio of the computation time for 1 to 30 cores for BD3MG for the restoration of Aneurysm.

Figure 6 :

 6 Figure 6: Slices (12, 5µm ˆ12, 5µm) for depths z " 5, 25 and 70 (from top to bottom) of the original acquisition (left) and after restoration (right). The comparison shows that the definition of the muscular structure has been enhanced by the reconstruction.

Table 1 :

 1 Characteristics and performance of compared algorithms on the Flybrain restoration task.

	Version Block-parallel Asynchrony MM SNR (dB)	f px f q	Error Time (s.)
	3MG	✗	✗	✓	18.58	1266.04 76.89	1683.79
	ABGD	✓	✓	✗	18.14	1268.80 76.73	305.76
	BP3MG	✓	✗	✓	18.11	1264.08 75.33	489.99
	BD3MG	✓	✓	✓	18.33	1261.59 74.72 147.16

Table 2 :

 2 Performance of BP3MG and BD3MG under imbalanced computed power, for a stopping criterion ε " 5 ˆ10 ´4 for Flybrain restoration. We additionally provide results for the vanilla 3MG algorithm for sake of comparison.

	3 cores hold a delay following a uniform distribution Upr0, 0.5sq.
	3 cores hold a delay following a uniform distribution Upr0, .25sq.
	' Type III: All worker cores are affected by a delay that follows a uniform
	distribution Upr0, 1sq.
	Method (Scenario)	SNR (dB)	f px f q	Time (s.)
	3MG (no delay)	18.132	1247.01	1683.79
	BP3MG (Type I)	17.941	1247.14	623.07
	BD3MG (Type I)	18.679	1246.04	211.34
	BP3MG (Type II)	17.941	1247.14	707.92
	BD3MG (Type II)	18.681	1246.03	220.65
	BP3MG (Type III)	17.941	1247.14	752.83
	BD3MG (Type III)	18.670	1246.02	219.90

This research work received funding support from the European Research Council Starting Grant MAJORIS ERC-2019-STG-850925.

Funding: This research work received funding support from the European Research Council Starting Grant MAJORIS ERC-2019-STG-850925.

Proof. Let k ě 2τ . Assumption 2 allows us to bound the gradient of f at x k , as

Let us extract the root of the above terms, and use triangular and gradient-Lipschitz inequalities, leading to

For all ℓ P k ´τ, k ´1 , by Proposition 4.1, ι ℓ ě ℓ ´τ `1 ě k ´2τ . Thus,

The right term of [START_REF] Jezierska | Spatially variant psf modeling in confocal macroscopy[END_REF] does not depend on index ℓ. Using [START_REF] Jezierska | Spatially variant psf modeling in confocal macroscopy[END_REF] and inequality [START_REF] Florescu | A majorize-minimize memory gradient method for complex-valued inverse problems[END_REF] finally proves the result.

Lemma 3.3 is useful as it provides a bound on the gradient at step x k only depending on the 2τ `1 past iterates x k , . . . , x k´2τ . Lemma 3.4. Let pu k q kPN , pv k q kPN two sequences of positive real. If there exists P P N and k ˚ě P such that

with r ă 1{P and

with

Plugging ( 47) into (46), yields

Proof. Function L is differentiable. The expression of its gradient is

´@Z P R pτ `1qN ¯∇LpZq " g 0 `L?

where g 0 " ˜∇f pZ 0 q

Let us apply twice the Jensen inequality for the square of the norm and then the majoration τ ´ℓ `1 ď τ for 1 ď ℓ ď τ . This yields

Using ? a 2 `b2 ď a `b for the two quantities at the right of ( 63) and then standard norm majoration inequalities, we get:

´@Z P R pτ `1qN ¯}∇LpZq} ď ?

The application of (64) to sequence pZ k q kPN leads to p@k ě τ q }∇LpZ k q} ď

By Lemma 3.3 and (65), we finally deduce that p@k ě 2τ q }∇LpZ k q} ď

which concludes the proof taking ρ " ? 2Lτ `?2ν `2Lτ 2 .

The following analysis makes use of recent theoretical results around the KL inequality [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] that we recall hereafter. For every ζ ą 0, we denote by Φ ζ the set of concave functions φ : r0, ζq Þ Ñ R `verifying : Competing interests: The authors have no relevant financial or non-financial interests to disclose.