Thomas Moreau

Alexandre Gramfort

DiCoDiLe: Distributed Convolutional Dictionary Learning

Keywords: Convolutional Dictionary Learning, Distributed System, Coordinate Descent. !

Convolutional dictionary learning (CDL) estimates shift invariant basis adapted to represent signals or images. CDL has proven useful for image denoising or inpainting, as well as for pattern discovery on multivariate signals. Contrarily to standard patch-based dictionary learning, patterns estimated by CDL can be positioned anywhere in signals or images. Optimization techniques consequently face the difficulty of working with extremely large inputs with millions of pixels or time samples. To address this optimization problem, we propose a distributed and asynchronous algorithm, employing locally greedy coordinate descent and a soft-locking mechanism that does not require a central server. Computation can be distributed on a number of workers which scales linearly with the size of the data. The parallel computation accelerates the parameter estimation and the distributed setting allows our algorithm to be used with data that do not fit into a single computer's RAM. Experiments confirm the theoretical scaling properties of the algorithm. This allows to demonstrate an improved pattern recovery as images grow in size, and to learn patterns on images from the Hubble Space Telescope containing tens of millions of pixels.

INTRODUCTION

C ONVOLUTIONAL sparse linear models consider that signals or images are obtained by linear combination of a few patterns localized in time or space. Patterns, which typically have a finite support, encode local variations such as edges in images or transient events in signals. The benefit of this model over traditional sparse coding technique based on local patches is that it is shift-invariant. While sparse coding encodes every local patch in the data, the convolutional linear model selects and encodes few patches leading to sparser and sharper reconstruction [START_REF] Bergstra | The Statistical Inefficiency of Sparse Coding for Images (or One Gabor to Rule them All)[END_REF]. Over the last ten to fifteen years, these models have been used successfully for various signal and image processing applications. It was first proposed by Grosse et al. [START_REF] Grosse | Shift-Invariant Sparse Coding for Audio Classification[END_REF] for music recognition. It was then used extensively for denoising and inpainting of images [START_REF] Bristow | Fast convolutional sparse coding[END_REF][START_REF] Kavukcuoglu | Learning Convolutional Feature Hierarchies for Visual Recognition[END_REF][START_REF] Papyan | Convolutional Dictionary Learning via Local Processing[END_REF][START_REF] Wohlberg | Efficient Algorithms for Convolutional Sparse Representations[END_REF]. Beyond denoising or inpainting, recent applications used this shiftinvariant model as a way to learn and localize relevant patterns in signals or images. Yellin et al. [START_REF] Yellin | Blood cell detection and counting in holographic lens-free imaging by convolutional sparse dictionary learning and coding[END_REF] used it to count the number of blood-cells present in holographic lens-free images. Jas et al. [START_REF] Jas | Learning the Morphology of Brain Signals Using Alpha-Stable Convolutional Sparse Coding[END_REF] and Dupré la Tour et al. [START_REF] Dupré La Tour | Multivariate Convolutional Sparse Coding for Electromagnetic Brain Signals[END_REF] learned recurring patterns in univariate and multivariate signals from neuroscience. This model has also been used for denoising astronomical data [START_REF] Starck | Image Decomposition Via the Combination of Sparse Representation and a Variational Approach[END_REF] and localizing predefined patterns from space satellite images [START_REF] Del Aguila Pla | Convolutional Group-Sparse Coding and Source Localization[END_REF]. For both denoising and pattern recovery tasks, there is a clear advantage in using large images. It is an expected statistical effect that when more realisations of the patterns are present, it is easier to disentangle them from other patterns and from the noise. Thus, for all aforementioned applications, there is a need for scalable and efficient algorithm that can handle large images.

Convolutional dictionary learning (CDL) is the technique commonly employed to estimate the parameters of such shiftinvariant models. The dictionary is here formed by the set of patterns (or atoms), whose positions in time or space are encoded by some sparse activation vectors. CDL then boils down to an optimization problem, where the variables are the dictionary as well as the activations. This problem is nonconvex and approximate solutions can be found by using alternating minimization over each variable. The estimation of • T. Moreau and A. Gramfort are with Université Paris-Saclay, Inria, CEA, Palaiseau, 91120, France.

the activation -which is commonly referred to as convolutional sparse coding (CSC) -is the critical bottleneck when it comes to applying CDL to large data.

As CSC is a convex problem, classical optimization tools can be employed. Chalasani et al. [START_REF] Chalasani | A fast proximal method for convolutional sparse coding[END_REF] used an accelerated proximal gradient method based on Fast Iterative Soft-Thresholding Algorithm (FISTA; [2]). Bristow et al. [START_REF] Bristow | Fast convolutional sparse coding[END_REF] proposed the Fast Convolutional Sparse Coding (FCSC) method, a splitting algorithm based on Alternating Direction Method of Multipliers (ADMM; [START_REF] Gabay | A dual algorithm for the solution of nonlinear variational problems via finite element approximation[END_REF]). While it is possible to leverage the Fast Fourier Transform to reduce the complexity of these approaches [START_REF] Wohlberg | Efficient Algorithms for Convolutional Sparse Representations[END_REF], both algorithms require to compute Fourier transforms of the entire data at each iteration. The cost can therefore be prohibitive for large signals or images. To avoid such costly global operations, several works have proposed to use local computations [START_REF] Grosse | Shift-Invariant Sparse Coding for Audio Classification[END_REF][START_REF] Kavukcuoglu | Learning Convolutional Feature Hierarchies for Visual Recognition[END_REF][START_REF] Moreau | DICOD: Distributed Convolutional Sparse Coding[END_REF][START_REF] Papyan | Convolutional Dictionary Learning via Local Processing[END_REF][START_REF] Zisselman | A Local Block Coordinate Descent Algorithm for the Convolutional Sparse Coding Model[END_REF]. A first line of work proposed by Kavukcuoglu et al. [START_REF] Kavukcuoglu | Learning Convolutional Feature Hierarchies for Visual Recognition[END_REF] uses greedy coordinate descent (GCD; [START_REF] Osher | Coordinate descent optimization for $\ell 1$ minimization with application to compressed sensing; a greedy algorithm[END_REF]) with an efficient mechanism to compute the update value. In their seminal work, Grosse et al. [START_REF] Grosse | Shift-Invariant Sparse Coding for Audio Classification[END_REF] propose to adapt the Feature Sign Search algorithm to only process small window of the signal at a time. Similarly, Papyan et al. [START_REF] Papyan | Convolutional Dictionary Learning via Local Processing[END_REF] proposed to solve the problem locally over smaller continuous chunks of signal and then aggregate these local results by solving a global optimization problem. In our recent work [START_REF] Moreau | DICOD: Distributed Convolutional Sparse Coding[END_REF], we proposed to use locally greedy coordinate descent (LGCD).

Based on updates similar to GCD, this method has a lower periteration complexity as it relies on local computation in the data. It does so like a patch-based technique would do, yet solving the non-local problem.

To further reduce the computational burden of CDL, recent studies consider parallelization strategies. Skau and Wohlberg [START_REF] Skau | A Fast Parallel Algorithm for Convolutional Sparse Coding[END_REF] proposed a parallel method based on consensus ADMM. It does so by splitting the computation across the different atoms. This limits the number of cores that can be used to the number of atoms in the dictionary, and it still requires Fourier transforms on the entire data. Concurrently, we proposed an alternative parallel algorithm for CSC, called DICOD [START_REF] Moreau | DICOD: Distributed Convolutional Sparse Coding[END_REF], that distributes GCD updates when working with one dimensional signals. DICOD could be used with multidimensional data such as images, yet only by splitting the data along one dimension, hence limiting its interest in this setting.

In the present paper, we propose a distributed solver for CDL, which can be used to learn recurring patterns in large multidimensional data. For the CSC step, we extend DICOD to use local updates with lower complexity while still working in the asynchronous setting without a central server. We detail a soft-lock mechanism, necessary to ensure the convergence of the algorithm with a grid of workers while remaining asynchronous. Then, we explain how to also leverage the distributed set of workers to further reduce the cost of the dictionary update. This is done by precomputing in parallel sufficient statistics used in the gradient operator. Extensive numerical experiments confirm the agreement between the theoretical complexity of our algorithm and its running time, allowing us to learn patterns from an astronomical image formed by tens of millions of pixels. A fully asynchronous implementation based on MPI and non-blocking communication primitives is made available to replicate our results.

In Section 2, we describe the convolutional linear model and the estimation of its parameters with CDL.

to Ti -1, |Ω| = d
i=1 Ti is its size. We will denote T = d i=1 Ti whenever the link between T and Ti is unambiguous. We denote X P Ω (resp. X P ×Q

Ω

) the space of observations defined on Ω with values in R P (resp. R P ×Q ). For instance, X 3 Ω with d = 2 is the space of RGB-images with height T1 and width T2. The value of an observation X ∈ X P Ω at position ω ∈ Ω is given by X[ω] = X[ω1, . . . , ω d ] ∈ R P , and its restriction to the p th coordinate at each position is denoted Xp ∈ X 1 Ω . All signals are 0-padded, i.e. for ω ∈ Ω, X[ω] = 0. The convolution operator is denoted * and is defined for

X ∈ X P Ω and Y Y Y ∈ X P ×Q Ω as the signal X * Y Y Y ∈ X Q Ω with (X * Y Y Y )[ω] = τ ∈Ω X[ω -τ ] • Y Y Y [τ ] , ∀ω ∈ Ω . (1)
For any signal X ∈ X P Ω , the reversed signal is defined as

X [ω] = X[T1 -ω1, . . . , T d -ω d ] .
Computing the convolution with X is equivalent to computing the cross-correlation with X. The p-norm is defined as X p = ω∈Ω X[ω] p p 1/p . We denote F the Discrete Fourier transform (DFT) of a signal and for a multivariate signal X ∈ X P Ω , X = F(X) denotes the DFT of each of its coordinates. Finally, we will also denote ST the soft-thresholding operator defined as

ST(u, λ) = sign(u) max(|u| -λ, 0) .

CONVOLUTIONAL DICTIONARY LEARNING

In this section, we recall the convolutional linear model and the classical way to estimate its parameters by alternating minimization.

Convolutional Linear Model

For an observed X ∈ X P Ω , signal or image, the convolutional linear model reads

X = Z * D D D + ξ , (2)
with As patterns typically occur at few locations in each observation, the coding signal Z is considered to be sparse. The aim of this model is truly to separate the shape of the local variation in the observationencoded in the dictionary -from their localization in space or time. Figure 1 illustrates the decomposition of a temporal signal with one pattern.

D D D ∈ X K×P Θ a dictionary of K atoms D D D k in dimension P defined on a sub-domain Θ = d i=1 [0, Li[⊂ Ω,
The parameters of this model are estimated by solving the following CDL optimization problem:

min Z,D D D D D D k ≤1 1 2 X -Z * D D D 2 2 F (Z,D D D) + λ Z 1 G(Z) , (3)
where the first term F (Z, D D D) measures how well the model fits the data and the second term G(Z) promotes sparse activations. The regularization parameter λ > 0 balances these two terms. The constraint D k 2 ≤ 1 is used to avoid the scale ambiguity in our model due to the 1 constraint. Indeed, the 1-norm of Z can be made arbitrarily small by rescaling D by α ∈ R * + and Z by 1 α . This problem is bi-convex, i.e. it is not jointly convex in (Z, D D D) but it is convex in each of its coordinates. It is usually solved with an alternating minimization algorithm, updating at each iteration one of the block of coordinates, Z or D D D. As each of the sub-problems are convex, the alternated descent is guaranteed to decrease the cost function, converging to a critical point of (3). It is however not guaranteed to be a global minimizer.

Convolutional Sparse Coding (CSC)

Given a dictionary of patterns D D D, CSC aims to retrieve the sparse decomposition Z * associated to the signal X. When minimizing only over Z, [START_REF] Bergstra | The Statistical Inefficiency of Sparse Coding for Images (or One Gabor to Rule them All)[END_REF] becomes

Z = argmin Z 1 2 X -Z * D D D 2 2 + λ Z 1 . (4)
This boils down to a Lasso problem with a Toeplitz design matrix. Therefore, classical optimization techniques for the Lasso can be applied with the same convergence guarantees. In practice, these methods are adapted to leverage the convolutional structure of the problem, reducing the computational complexity to solve (4).

Locally Greedy Coordinate Descent (LGCD)

Coordinate descent (CD) is an algorithm which updates a single coefficient at each iteration. Following Moreau et al. [START_REF] Moreau | DICOD: Distributed Convolutional Sparse Coding[END_REF], we use locally greedy coordinate selection to solve efficiently the CSC (4). Algorithm 1 summarizes the computations for LGCD.

Coordinate update For (4), it is possible to compute efficiently the optimal update for a given coordinate when all the others are fixed. Denoting q the iteration number, and

Z (q)
k 0 [ω0] (atom k0 at position ω0) the optimal value of the activation Z (q) k 0 [ω0], the update reads

Z (q) k 0 [ω0] = 1 D D D k 0 2 2 ST(β (q) k 0 [ω0], λ) , (5)
where β (q) is an auxiliary variable in

X K Ω defined for (k, ω) ∈ [1, K] × Ω as β (q) k [ω] = X -Z (q) * D D D + Z (q) k [ω]eω * D D D k * D D D k [ω] , (6)
and where eω is a dirac (canonical vector) with value 1 in ω and 0 elsewhere. The value of the auxillary variable in (k0, ω0) corresponds to the gradient of F for a point Z (q) with the coordinate (k0, ω0) set to 0. To enable efficient coordinate update, an important property is that Z (q) [ω] only depends on a small neighborhood around ω.

Definition 2.1 (Θ-neighborhood). Let us consider a support

Θ = d i=1 [0, Li[. For Ω0 ∈ Ω, the Θ-neighborhood VΘ(ω0) is defined as VΘ(ω0) = d i=1 [ω0,i -Li + 1, ω0,i + Li[ .

Proposition 2.2 (Weak dependence).

If the coordinate (k0, ω0) is updated in Z (q) , then the value of Z (q+1) differs from Z (q) only for ω in a Θ-neighborhood VΘ(ω0) of ω0.

Proof.

If the coordinate Z k 0 [ω0] is updated with an additive update ∆Z k 0 [ω0] = Z (q) k 0 [ω0] -Z (q)
k 0 [ω0], Kavukcuoglu et al. [START_REF] Kavukcuoglu | Learning Convolutional Feature Hierarchies for Visual Recognition[END_REF] showed that β (q+1) can be obtained using the relation

β (q+1) k [ω] = β (q) k [ω] -(D D D k 0 * D D D k )[ω -ω0]∆Z (q) k 0 [ω0] , (7) for all (k, ω) = (k0, ω0). As the support of (D D D k 0 * D D D k )[ω] is VΘ(0), this implies that the value β (q+1) k only changes in the neighborhood VΘ(ω0) = ω0 + VΘ(0).
This proposition shows that the optimal value for a coordinate only depends on the current values of the activations in a neighborhood of size |VΘ(0)| = 2 d L. This corresponds to twice the size of the dictionary in each dimension. Using this property, the computational complexity of updating a coordinate is O(2 d KL), as it is only necessary to update β and Z (q) on VΘ(ω0).

Coordinate selection

The selection of the updated coordinate (k (q) 0 , ω (q) 0 ) can follow different strategies. Following our work in [START_REF] Dupré La Tour | Multivariate Convolutional Sparse Coding for Electromagnetic Brain Signals[END_REF][START_REF] Moreau | DICOD: Distributed Convolutional Sparse Coding[END_REF], we propose to select the coordinates with a locally greedy rule (LGCD). Definition 2.3 (Disjoint and contiguous partition). Let us consider a domain Ω. For M ∈ N, a M disjoint and contiguous partition {Cm} M m=1 is a set of subsets of Ω such that

d i=1 [ωi, ω i ] ⊂ Cm for ω, ω ∈ Cm , m Cm = Ω , Cm ∩ C m = ∅ for m = m .
At iteration q, the coordinate to update is selected greedily among the m-th sub-domain Cm with (k

(q) 0 , ω (q) 0 ) = argmax (k,ω)∈[1,K]×Cm Z (q) k [ω] -Z (q) k [ω] ∆Z (q) k [ω] , (8)
with m = q mod M . This strategy cycles through each of the M sub-domains and select the coordinate in Cm which is the farthest away from its optimal value when all the other coordinates are fixed. The quantity ∆Z k [ω] acts as a proxy for the reduction of the cost function obtained with the update of coordinate (k, ω) to the value [START_REF] Osher | Coordinate descent optimization for $\ell 1$ minimization with application to compressed sensing; a greedy algorithm[END_REF]. The complexity of the locally greedy selection rule is linear with the size of the sub-segments O(K|Cm|). In the case where M = KT , this rule boils down to the cyclic CD [START_REF] Friedman | Pathwise coordinate optimization[END_REF] with a coordinate selection complexity in O(1), while setting M = 1 boils down to greedy CD [START_REF] Osher | Coordinate descent optimization for $\ell 1$ minimization with application to compressed sensing; a greedy algorithm[END_REF] where the complexity of the selection is linear with the size of the signal O(KT ). The selection of M in LGCD acts as a tradeoff between the computational complexity and the convergence speed of these two methods.

Z (q) k [ω]
As it was shown in recent studies [START_REF] Karimireddy | Efficient Greedy Coordinate Descent for Composite Problems[END_REF][START_REF] Nutini | Coordinate Descent Converges Faster with the Gauss-Southwell Rule Than Random Selection[END_REF], updating in priority the important coordinates by greedy selection leads to faster convergence rates, although it increases the cost of one iteration. The LGCD strategy proposes to use sub-domains Cm such that the computational complexities of selecting and updating the coordinate match. By using sub-domains of the VΘ(0)'s size i.e. |Cm| = 2 d L, the global iteration complexity of LGCD scales with 2 d times the size of the dictionary D D D i.e. O(2 d KL). By doing so LGCD achieves faster convergence rates than nongreedy coordinate selection rule while limiting the per-iteration cost.

Stopping criterion

The stopping criterion used for this algorithm follows [START_REF] Osher | Coordinate descent optimization for $\ell 1$ minimization with application to compressed sensing; a greedy algorithm[END_REF]. The algorithm is stopped once the magnitude of all possible coordinate updates are below a threshold > 0 specified by the user, i.e.

max (k,ω)∈[1,K]×Ω |∆Z (q) k [ω]| < . (9)

Global solvers

Chalasani et al. [START_REF] Chalasani | A fast proximal method for convolutional sparse coding[END_REF] propose to solve (4) by proximal gradient methods exploiting the convolution in the gradient computation to save time. The core idea behind this method is to use a proximal splitting of the cost function between its smooth part F (•, D D D) and a proximable part G. Proximal gradient descent (ISTA) reads

Z (q+1) = prox γG Z (q+1) -γ∇Z F (Z (q+1) , D D D) , ( 10
)
where γ is the step size of the algorithm, usually chosen as the inverse of the smoothness constant of F , and prox γG is the proximal operator of γG. For G(•) = λ • 1, this corresponds to the soft-thresholding operator ST(•, γλ). To make this algorithm efficient the key is to rely on the FFT for the gradient computations, as described in [START_REF] Wohlberg | Efficient Algorithms for Convolutional Sparse Representations[END_REF]. Doing the computation in the

Algorithm 1 Locally Greedy Coordinate Descent Require: D D D, X, parameter > 0, sub-domains Cm, 1: Initialization: ∀(k, ω) ∈ [1, K] × Ω, 2: Z k [ω] = 0, β k [ω] = D D D k * X [ω] 3: repeat 4:
for m = 1 . . . M do 5:

∀(k, ω) ∈ [1, K] × Cm, Z (q) k [ω] = 1 D D D k 2 2 ST(β (q) k [ω], λ) , 6:
Choose (k0, ω0) = arg max

(k,ω)∈[1,K]×Cm |∆Z k [ω]| 7:
Update β using [START_REF] Combettes | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF] and

Z k 0 [ω0] ← Z (q) k 0 [ω0] 8:
end for 9: until ∆Z (q) ∞ <

Fourier domain, the computational complexity of each iteration of this algorithm is O(KT log T ). It is also possible here to use an accelerated version FISTA to solve the problem [2,[START_REF] Chalasani | A fast proximal method for convolutional sparse coding[END_REF].

An alternative global algorithm proposed by Bristow et al. [START_REF] Bristow | Fast convolutional sparse coding[END_REF] is based on the ADMM algorithm [START_REF] Gabay | A dual algorithm for the solution of nonlinear variational problems via finite element approximation[END_REF]. The core idea in this line of work is to introduce an auxiliary variable and split the two part of the objective with an equivalent problem

min Y ,Z 1 2 X -Y D D D 2 F + λ Z 1 s.t. F(Z) = Y . ( 11
)
This problem is then solved using ADMM in its scaled form [4, Section 3.1.1]. The critical part in this algorithm is the minimization of the augmented Lagrangian relatively to Y which necessitates multiple FFT as well as the inversion of T linear systems of size K × K. The overall complexity of this method is thus O(KT (log

T + K 2 )).
For these global solvers, the complexity of each iteration scales with the size of the data in O(T log T ) and the convergence rate is independent of T . In contrast, the iteration of LGCD are independent of T but the number of iterations needed to reach convergence scales with the number of non-zero coefficients which also depends on T . For large and sparse data, LGCD tends to be more efficient as the overall complexity depends linearly in T and not in T log T . However, the choice between a local vs a global algorithm depends on the sparsity level of the desired solution, controlled by the value of the regularization parameter λ [START_REF] Dupré La Tour | Multivariate Convolutional Sparse Coding for Electromagnetic Brain Signals[END_REF].

Dictionary Update

Given an activation signal Z (q) , the dictionary is updated to improve how the model reconstructs the signal. The problem reads

D D D * = argmin D D D k 2 ≤1 1 2 X -Z (q) * D D D 2 2 . ( 12
)
This problem is smooth and convex and can be solved using algorithms designed for smooth constrained optimization.

Projected Gradient Descent

The projected gradient descent (PGD; [START_REF] Combettes | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF]) and its accelerated version (APGD) are possible strategies. These two algorithms are equivalent to ISTA and FISTA iterations where the proximal operator is a simple projection [START_REF] Combettes | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF]. The efficiency of these two methods relies on pre-computations of sufficient statistics for fast evaluation of the gradient of the objective ∇ D D D F [START_REF] Dupré La Tour | Multivariate Convolutional Sparse Coding for Electromagnetic Brain Signals[END_REF]. The computations can be factorized as

∇ D D D F (Z, D D D) = Z * (X -Z * D D D) = Ψ Ψ Ψ -Φ Φ Φ * D D D , (13)
where

Φ Φ Φ ∈ X K×K V Θ (0) and Ψ Ψ Ψ ∈ X K×P Θ
are respectively defined as the convolution Z * Z restricted to VΘ(0) and the convolution Z * X restricted to Θ, i.e.

Φ Φ Φ[ω] = Z * Z = τ ∈Ω Z[ω + τ ]Z [τ ]
for ω ∈ VΘ(0) , ( 14)

Ψ Ψ Ψ[ω] = Z * X = τ ∈Ω Z[ω + τ ]X [τ ] for ω ∈ Θ . ( 15
)
The complexity to obtain the two constants Φ Φ Φ and Ψ Ψ Ψ scales as O(K(K+P )LT ). This is comparable with the cost of computing the gradient directly O(KP T log T ). However, once the two constants are pre-computed, the gradient of F is obtained using [START_REF] Giavalisco | The Great Observatories Origins Deep Survey: Initial Results From Optical and Near-Infrared Imaging[END_REF] with complexity O(2 d K 2 P L log(L)). This is independent of the data size T , making each iteration efficient.

Related solvers

The block coordinate descent for dictionary learning [START_REF] Mairal | Online Learning for Matrix Factorization and Sparse Coding[END_REF], which consists in updating one atom at a time, can also be extended to the convolutional case. As PGD for full dictionary updates, this technique benefits from sufficient statistics pre-computations and the complexity is similar. Another possibility is to rely on ADMM to update the dictionary [START_REF] Bristow | Fast convolutional sparse coding[END_REF][START_REF] Wohlberg | Efficient Algorithms for Convolutional Sparse Representations[END_REF]. However, ADMM algorithms can be slow to reach a convergence, yielding suboptimal reconstruction performances as the feasible point is not directly optimized but results from a consensus. Yellin et al. [START_REF] Yellin | Blood cell detection and counting in holographic lens-free imaging by convolutional sparse dictionary learning and coding[END_REF] proposed to adapt the K-SVD dictionary update method [START_REF] Aharon | K-SVD: An Algorithm for Designing Overcomplete Dictionaries for Sparse Representation[END_REF] to the convolutional setting but this technique is adapted to an 0 penalty for the activation Z. It also requires to compute costly singular value decompositions (SVD).

Choice of initialization and regularization parameter

In practice, the quality of the solution obtained by CDL is impacted by the choice of the regularization parameter λ. Also, due to the iterative nature of the solvers, the initialization of the dictionary D D D 0 can play an important role.

Regularization parameter selection

The choice of a regularization parameter λ in the context of dictionary learning is often selected by cross-validation on a subsequent task such as denoising or classification. However in the context of purely unsupervised learning, e.g. matrix factorization, this is still an open problem. This problem has recently received some attention using Bayesian methods [START_REF] Pereyra | Maximum-A-Posteriori Estimation with Unknown Regularization Parameters[END_REF][START_REF] Vidal | Maximum Likelihood Estimation of Regularization Parameters[END_REF]. For the specific case of CDL using an 1 norm, it is useful to note that there exists a value λmax such that for λ ≥ λmax, 0 is solution of (4). The first order condition for (4) reads

∇Z F (0, D D D (0) ) = X * D D D (0) ∈ ∂G(0) = [-λ, λ] KT . ( 16
)
This condition is always verified if

λ > λmax = X * D D D (0) ∞.
In the following, the regularization parameter λ is chosen as a fraction of λmax. Using this data dependent value, the ratio λ/λmax is chosen in [0, 1], making the definition of the regularization parameter common to all input signals X.

Dictionary initialization

The standard strategy to initialize D D D (0) for CDL methods is to generate random atoms with Gaussian white noise. However, as these atoms generally poorly correlate with the data, the initial value of λmax is low compared to the following ones, leading to dense activation Z. For dictionary learning, a classical way to initialize the dictionary is to take samples from the input data [START_REF] Mairal | Online Learning for Matrix Factorization and Sparse Coding[END_REF]. In the convolutional setting, this boils down Compute Z (q+1) with DiCoDiLeZ (X, D D D (q) , W ).

3:

Compute φ and ψ with ( 14) and W workers.

4:

Update D D D (q+1) with PGD and Armijo backtracking linesearch. 5: until Cost variation is smaller than ν to taking random patches of the size of the dictionary directly from the data. This can be theoretically motivated by the work of Zhang et al. [START_REF] Zhang | Structured local minima in sparse blind deconvolution[END_REF] that demonstrate improvements in pattern recovery under suitable assumptions.

DISTRIBUTED CONVOLUTIONAL DICTIONARY LEARN-ING (DICODILE)

In this section, we introduce DiCoDiLe, a distributed algorithm for convolutional sparse coding. Its steps are summarized in Algorithm 2. The CSC part of this algorithm will be referred to as DiCoDiLeZ . It extends our previous work DICOD [START_REF] Moreau | DICOD: Distributed Convolutional Sparse Coding[END_REF] by handling multidimensional data, in particular images. The algorithm proposes to distribute the computation on a grid of workers. To guarantee convergence with multidimensional data in an asynchronous setting, we introduce the notion of soft-locks described in Subsection 3.1. Details on how the dictionary updates can also benefit from the distributed setting are presented in Subsection 3.2.

DiCoDiLeZ : Distributed CSC with coordinate descent

In the convolutional setting, coordinates that are far enough -relatively to the size of the dictionary -are only weakly dependent (Proposition 2.2). It is thus natural to distribute the computation of CD by splitting the data in continuous subdomains [START_REF] Moreau | DICOD: Distributed Convolutional Sparse Coding[END_REF].

The DiCoDiLeZ algorithm, described in Algorithm 3, is distributed over W workers which update asynchronously the coordinates of Z. We consider {Sw} W w=1 , a disjoint and contiguous partition of Ω. Each worker w ∈ [1, W ] is in charge of updating the coordinates of the sub-domain Sw. These sub-domains are also partitioned with disjoint and contiguous sub-domains C (w) m of size |VΘ(0)|. Then, the worker w chooses at each iteration q an update candidate (k0, ω0) which corresponds to a greedy update in one of its sub-domain

C (w) m , i.e. (k0, ω0) = argmax (k,ω)∈[1,K]×C (w) m |∆Z k [ω]| , (17)
with m = q mod M . The sub-domains C (w) m are selected in a cyclic order but empirical test showed no significant differences for random sub-domain selection. Then, this candidate update is accepted or rejected using the soft-lock mechanism described in Subsection 3.1.1. If the candidate is rejected, the worker moves on to the next sub-domain C w m+1 without updating a coordinate. If the candidate is accepted, the value of the coordinate

Z k 0 [ω0] is replaced by Z (q)
k 0 [ω0] and β is updated with [START_REF] Combettes | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF]. To ensure the convergence, it is sufficient to notify the neighboring workers if a local update changes the value of β outside of Sw, i.e. if VΘ(ω0) ⊂ Sw. Receive messages and update Z and β with [START_REF] Combettes | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF] 8:

Choose (k0, ω0) = argmax (k,ω)∈[1,K]×C (w) m |∆Z k [ω]| 9: if |∆Z k 0 [ω0]| < max (k,ω)∈[1,K]×V Θ (ω 0 ) |∆Z k [ω] then 10:
The coordinate is soft-locked, goto 6

11:

end if 12:
Update β with [START_REF] Combettes | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF] and

Z k 0 [ω0] ← Z (q) k 0 [ω0] 13:
if ω0 ∈ B2L(Sw) then end for 17: until global convergence ∆Z ∞ < sub-domain for dimension i. For Θ = d i=1 [0, Li[, the Θ-border is defined as the set of points ω ∈ Sw such that VΘ(ω) ⊂ Sw, i.e.

BΘ(Sw) = d i=1 [li, li + Li[∪[ui -Li, ui[ . (18)
In the case where the updated coordinate (k, ω0) is in Sw \ BΘ(Sw), then VΘ(ω0) is included in the set of coordinates Sw updated by the worker w. In this case, the update is not at the border of the sub-domain, so the values of β and Z (q) need no update for the other workers w . Thus, there is no need to exchange information between the workers. However, for ω0 ∈ BΘ(Sw), it is necessary to notify other workers w whose domain overlap with the Θ-neighborhood of ω0 i.e. VΘ(ω0) ∩ S w = ∅. The notification can be done by sending the triplet (k0, ω0, ∆Z k 0 [ω0]) to these workers. These information are sufficient to keep β up to date for all workers using [START_REF] Combettes | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF]. Figure 2 illustrates this communication process on a 2D example. With this distributed algorithm, the inter-processes communications are minimal, as only d + 2 values need to be exchanged. Moreover, as it does not rely on centralized communication, the necessary bandwidth between the workers is also reduced as a worker only communicates with its direct neighbors. This contrasts with other distributed algorithms for CDL such as [START_REF] Skau | A Fast Parallel Algorithm for Convolutional Sparse Coding[END_REF] that requires centralized communication of full support values between the workers. Finally, the communication only occurs when ω0 ∈ BΘ(Sw), which happens rarely especially when working with large data and patterns of small support. This reduces even further the need of bandwidth between workers.

Interferences and Soft-Locks in asynchronous setting

While it is now clear how communication should be done between workers, one need to clarify how to proceed to ensure convergence in the case of asynchronous updates. When W coordinates (kw, ωw) W w=1 of Z are updated simultaneously by respectively ∆Z kw [ωw], the updates might not be independent. The local version of β used for the update does not account for the other updates. The cost difference resulting from all these simultaneous updates is denoted ∆E, and the cost reduction induced by only one update w is denoted

∆E kw [ωw]. Simple S w S w-1 S w+1 S w-W1-1 S w-W1 S w-W1+1 S w+W1 S w+W1+1 T 1 /W 1 T 2 /W 2 L 1 L 1 ω 0 • ω 1 • k 1 , ω 1 , ∆Z k1 [ω 1 ] ω 2 • k 2 , ω 2 , ∆Z k2 [ω 2 ] S w V Θ (ω i ) B Θ (S w ) E Θ (S w ) Soft-
-

w =w (D D D kw * D D D k w )[ω w -ωw]∆Z kw [ωw]∆Z k w [ω w ]) interference ,
If for all w, all other updates w are such that if ω w / ∈ VΘ(ωw), then (D D D kw * D D D k w )[ω w -ωw] = 0. The updates can be considered to be sequential as the interference term is zero. This also directly results from Proposition 2.2. When {ωw} w ∩ VΘ(ω0) = I0 > 1, the interference term does not vanish. Moreau et al. [START_REF] Moreau | DICOD: Distributed Convolutional Sparse Coding[END_REF] show that if I0 < 3, then, under mild assumption, the interference term can be controlled and DICOD converges. While this is sufficient for 1D partitioning as I0 cannot be larger than 2 (if Sw is larger than VΘ(0)), this analysis cannot be extended to I0 ≥ 3. This limits the partitioning of Ω along one single dimension.

To avoid this limitation, we propose the soft-lock, a mechanism to avoid interfering updates of higher order. The idea is to avoid updating concurrent coordinates simultaneously. A classical tool in computer science to address such concurrency issue is to rely on a synchronization primitive, called lock. This primitive can ensure that at most one worker enters a protected block of code simultaneously. In our case, it would be possible to use this primitive to restrict updates in BΘ(Sw) only when none of its neighbors is updating BΘ(S w ). This would avoid interfering updates. However, the drawback of this method is that it prevents the algorithm to be asynchronous, as typical lock implementations rely on centralized communications. Instead, we propose to use an implicit locking mechanism which keeps our algorithm asynchronous.

Definition 3.2 (Θ-extension). We consider

Sw = d i=1 [li, ui[ a sub-domain of Ω. For Θ = d i=1 [0, Li[, the Θ-extension of Sw is defined as EΘ(Sw) = (Sw + VΘ(0)) \ Sw, i.e. EΘ(Sw) = d i=1 [li -Li, li[∪[ui, ui + Li[ . ( 20
)
The idea of the soft-lock is the following. Additionally to the sub-domain Sw, each worker also maintains the value of Z (q) and β on the Θ-extension EΘ(Sw) (See Figure 2). These extra coordinates are referred to as ghost cells in the computer science literature [START_REF] Kjolstad | Ghost Cell Pattern[END_REF]. At each iteration, the worker w considers a candidate coordinate (k0, ω0) to update from [START_REF] Kavukcuoglu | Learning Convolutional Feature Hierarchies for Visual Recognition[END_REF]. If ω0 ∈ BΘ(Sw), the coordinate is updated like in LGCD. When ω0 ∈ BΘ(Sw), the candidate is accepted if there is no better coordinate update in VΘ(ω0) ∩ EΘ(Sw), i.e. if

|∆Z k 0 [ω0]| > max k,ω∈[1,K]×V Θ (ω 0 )∩E Θ (Sw ) |∆Z k [ω]| . (21)
In case of equality, the update in the sub-domain S w with the minimal index w is preferred. Using this mechanism effectively prevents any interfering update. If two workers w < w have update candidates ω0 and ω 0 such that ω 0 ∈ VΘ(ω0), then with ( 21), only one of the two updates will be accepted. It will be either the largest one if

|∆Z k 0 [ω0]| = |∆Z k 0 [ω 0 ]
| or the one from w if both updates have the same magnitude. The other update will be rejected and there will be no interference. By doing so this mechanism is fully asynchronous. The updates on the borders are implicitly ordered, preventing interference. The update are guaranteed to be performed as if they were run sequentially, ensuring convergence of the algorithm. Unlike classical ghost cells [START_REF] Kjolstad | Ghost Cell Pattern[END_REF], the soft-lock mechanism for DiCoDiLe does not require a synchronization step after each update in a worker. The update of the ghost cells in neighboring workers are handled using the asynchronous messages used to maintain β up to date (see Algorithm 3 l.14). The soft-lock mechanism requires to increase slightly the communication between the workers and to double the size of the Θ-border in each direction. This is however not a bottleneck as demonstrated in the following experiments.

Speed-up analysis of DiCoDiLeZ

In this section, we consider the speed-up provided by this distributed algorithm. As each worker runs LGCD locally, the computational complexity of a local iteration in a worker w is O(2 d KL). This complexity is independent of the number of workers used to run the algorithm. Thus, the speed-up obtained is equal to the number of updates that are performed in parallel. If we consider W update candidates (kw, ωw) chosen by the workers, the number of updates that are performed corresponds to the number of updates that are not soft-locked.

To evaluate this, we make the assumption that the updates are uniformly distributed on each Sw. In this case, the probability that an update candidate located in ω ∈ Sw is not soft-locked can be lower bounded by

P (ω ∈ SL) ≥ d i=1 1 - WiLi Ti (22)
Computations are detailed in Proposition B.1. When the ratios T i W i L i are large for all i, i.e. when the size of the workers subdomains |Sw| are large compared to the dictionary size L, then P (ω ∈ SL) In comparison, we report in [START_REF] Moreau | DICOD: Distributed Convolutional Sparse Coding[END_REF] a super-linear speed-up for DICOD on 1D signals. This is due to the fact that we are using GCD locally in each worker. GCD has a very high iteration complexity when W is small, as it scales linearly with the size of the considered sub-domain. Thus, when sub-dividing the signal between more workers, the iteration complexity is decreasing and more updates are performed in parallel. This explains why the speed-up is almost quadratic for low W . As the complexity of iterative GCD are worse than LGCD, DiCoDiLeZ has better performance than DICOD in low W regime. Furthermore, in the 1D setting, DICOD becomes similar to DiCoDiLeZ once the size of the workers sub-domain becomes too small. This happens when the local domain of LGCD becomes comparable to the size of the sub-domains Sw. This implies that DiCoDiLeZ always outperforms DICOD.

Stopping criterion

In this distributed setting, it is non-trivial to know when the algorithm reaches convergence. To match the sequential stopping criterion of GCD ( 9), the convergence is considered to be reached once no worker can perform a coordinate update |∆Z (q) | that is higher than > 0. However, as a worker can be affected by its neighboring workers, the stopping criterion of CD cannot be applied independently in each worker. Once a worker reaches this state, the potential update can still increase above because of updates in EΘ(Sw) by another worker. To make sure the correct criterion is achieved, workers that reach the local criterion are paused. They are restarted when a neighboring workers triggers an update on EΘ(Sw). The global convergence can be decided when all workers are paused. While this process is based on centralized communication, it can still be performed asynchronously using message passing to a designated master process that will notify the end of the algorithm to all workers.

Dictionary updates in a distributed setting

For the dictionary, a change in a coordinate in Z at any point in Ω can impact the solution of [START_REF] Gabay | A dual algorithm for the solution of nonlinear variational problems via finite element approximation[END_REF]. Moreover, the atoms are interdependent of each other. It is thus complicated to parallelize the dictionary update as it is centralized by nature and there is no natural division of the computation. However, the size of D D D is way smaller than Z and ( 12) can be solved efficiently when applying the PGD algorithm described in Subsection 2.3. The main computational bottleneck is the computation of the constants Φ Φ Φ and Ψ Ψ Ψ. For these, it is possible to leverage the distributed setting as the computation are separable between the different workers. The computations can be made embarrassingly parallel among the workers as

Φ Φ Φ[ω] = W w=1 τ ∈Sw Z[τ ]Z [τ + ω] for ω ∈ VΘ(0) , (23)
Ψ Ψ Ψ[ω] = W w=1 τ ∈Sw Z[τ ]X [τ + ω] for ω ∈ Θ . ( 24
)
For all (ω, τ ) ∈ VΘ(0) × Sw, we have ω + τ ∈ Sw ∪ EΘ(Sw).

Thus, the computations of the inner sum Φ Φ Φ w can be made independently by each worker. These partial results are gathered and reduced with Φ Φ Φ = W w=1 Φ Φ Φ w . By making use of these distributed computation, Φ Φ Φ and Ψ Ψ Ψ can be computed with respectively O(K 2 |Sw| log(|Sw ∪ EΘ(Sw)|)) and O(KP |Sw| log(|Sw ∪ EΘ(Sw)|)) operations. The values of Ψ Ψ Ψ and Φ Φ Φ are sufficient statistics to compute ∇ D D D F and F . The dictionary can thus be updated without having to gather the full activation Z on one machine. This allows to learn a dictionary for a signal X that does not fit in memory, as each worker only needs to store a part of the full signal and the dictionary updates only require storing these sufficient statistics.

NUMERICAL EXPERIMENTS

The numerical experiments are performed on a SLURM cluster with 30 nodes. Each node has 20 physical cores, 40 logical cores and 250 GB of RAM. A worker is typically one logical core on one of the machine. In practice we use up to 400 workers. DiCoDiLe is implemented in Python [START_REF]Python Software Fundation[END_REF] using

mpi4py [8] 1 .
Experiments make use of simulated one dimensional signals, generated following the sparse convolutional linear model ( 2) with d = 1 and P = 7. The dictionary is composed of K = 25 atoms D D D k of length L = 250. Each atom is sampled from a standard Gaussian distribution and then normalized. The sparse code entries are drawn from a Bernoulli-Gaussian distribution, with Bernoulli parameter ρ = 0.007, mean 0 and standard variation 10. The noise term ξ is sampled from a standard Gaussian distribution with variance 1. The length of the signals X is denoted T . For experiments with images, we will use the standard colored image Mandrill (P = 3) at full resolution (512 × 512) as well as text images generated with L A T E X and ImageMagick [32]. Finally, for large scale experiments, we used an image of the Great Observatories Origins Deep Survey (GOODS) South field, acquired by the Hubble Space Telescope [START_REF] Giavalisco | The Great Observatories Origins Deep Survey: Initial Results From Optical and Near-Infrared Imaging[END_REF]. We used the STScI-H-2016-39 image 2 with resolution 6000 × 3600.

In all experiments, we used λ = 0.1λmax as a default regularization parameter, unless specified otherwise.

Performance of distributed sparse coding

The performances of the distributed sparse coding step in DiCoDiLe is impacted by the choice of coordinate descent strategy, as well as soft-lock mechanism and the data partitioning scheme. The following experiments discuss their effects on the scalability on DiCoDiLeZ .

Coordinate Selection Strategy

Using one worker, Figure 3 compares the average running times of three CD strategies to solve (4): Cyclic, Greedy and LGCD.

LGCD consistently outperforms the other strategies for both T = 150L and T = 750L. For LGCD and Cyclic, the iteration complexity is constant compared to the size of the signal, whereas the complexity of each Greedy iteration scales linearly with T. This makes running time of Greedy explode when T grows. Besides, as LGCD also prioritizes the more important coordinates, it is more efficient than Cyclic that selects the coordinate to update independently of the data. This explains the performances observed in Figure 3 and justifies the use of LGCD as the base of the DiCoDiLeZ distributed algorithm. In the following, we will not include Cyclic coordinate selection in our distributed algorithm as it is less efficient. We include Greedy as its complexity depends on the size of the signal, meaning that it can be competitive when the workload is distributed on many workers.

Need for soft-lock mechanism

The reconstruction of the Mandrill with and without soft-locks is shown in Figure 4. It uses a dictionary D D D consisting of K = 25 patches of size 16×16 extracted randomly from the original image. The reconstructed image Z * D D D obtained without soft-lock shows artifacts around the edges of some of the sub-domains Sw, as the activations diverge in these locations. This visually demonstrates the need for controlling the interfering updates when more than two workers can interfere. When using the soft-lock mechanism, the algorithm returns the expected result.

Impact of data partitioning

To demonstrate the impact of the domain partitioning strategy on images, Figure 5 shows the average runtime for DiCoDiLeZ using either a line of worker (unidirectional partitioning) or a grid of workers (bidirectional partitioning). In this experiment K = 5 and the atom size is 8 × 8. Both strategies scale similarly in the regime with low number of workers. But when W reaches T1/3L1, the performances of DiCoDiLeZ stops improving when using the unidirectional partitioning of Ω. Indeed in this case many update candidates are selected at the border of the sub-domains Sw, therefore increasing the chance to reject an update. Moreover, the scaling of the unidirectional partitioning is limited to W = T1/2L1 = 32, whereas the bidirectional partitioning can be used with W up to 1024 workers. Linear Split Grid Split Fig. 5. Scaling on 2D images of DiCoDiLe Z with the number of workers for two partitioning strategies of Ω: (blue) Ω is split only along one direction, as in DICOD [START_REF] Moreau | DICOD: Distributed Convolutional Sparse Coding[END_REF], (orange) Ω is split along both directions, on a grid. The running times are similar for low number of workers but the performances with the linear splits stop improving once W reaches the scaling limit T 1 /4L 1 (green). With the grid of workers adapted to images, the performance improves further. The linear split stops when W reached T 1 /L 1 because no more coordinate in Sw are independent from other neighbors.

T = 150L

Scaling of DiCoDiLeZ

Scaling for 1D signals Figure 6 investigates the role of CD strategies in the distributed setting with 1D signals. The average runtime of DICOD (with GCD on each worker), and DiCoDiLeZ (with LGCD) are reported in Figure 6. DiCoDiLeZ with LGCD scales sub-linearly with the number of workers used, whereas DICOD scales almost quadratically. However, DiCoDiLeZ outperforms DICOD consistently, as the performance of GCD is poor for low number of workers using large sub-domains. The two coordinate selection strategies become equivalent when W reaches T /4L, as in this case DiCoDiLe 0:1¸m ax 0:2¸m ax 0:5¸m ax

LGCD GCD Fig. 7. Scaling of DiCoDiLe Z with the number of workers W for different values of λ and for two strategies of coordinate selection: Greedy and Locally Greedy. The convergence is faster when λ is large because the activation becomes sparser and less coordinates need to be updated. Also, the locally greedy coordinate selection outperforms the greedy selection up to a certain point where the performances become similar as the sub-domain Sw becomes too small to be further partitioned with sub-domains of size Θ.

has only one sub-domain

C (w)
1 . This is consistent for larger signals. Indeed, in the case where T = 750, DiCoDiLeZ scales linearly with the number of workers. We do not see the drop in performance, as in this experiments W does not reach the scaling limit T /4L. Importantly, DiCoDiLeZ still outperforms DICOD for all regimes.

Scaling for 2D signals Figure 7 illustrates the scaling performance of DiCoDiLeZ on images. The average running times of DiCoDiLeZ are reported as a function of the number of workers W for different regularization parameters λ and for greedy and locally greedy selection. As for the 1D CSC problem, the locally greedy selection consistently outperforms the greedy selection. This is particularly the case when using low number of workers W . As for signals, once the size of the sub-domains becomes smaller, the performances of both coordinate selection methods become closer as both algorithms become equivalent. One can also notice that the convergence of DiCoDiLeZ is faster for large λ. This is expected, as in this case, the solution Z * is sparser and less coordinates need to be updated.

Comparison with Skau and Wohlberg [28]

Figure 8 displays the evolution of the objective function (3) as a function of time for DiCoDiLe and the Consensus ADMM algorithm 3 proposed by Skau and Wohlberg [START_REF] Skau | A Fast Parallel Algorithm for Convolutional Sparse Coding[END_REF]. To the best of our knowledge, the work of Skau and Wohlberg [START_REF] Skau | A Fast Parallel Algorithm for Convolutional Sparse Coding[END_REF] is the only approach in the literature that proposes a parallel algorithm for CDL. Both algorithms were run on a single node with 36 workers. We used the Hubble Space Telescope GOODS South image as input data. As the algorithm by Skau and Wohlberg [START_REF] Skau | A Fast Parallel Algorithm for Convolutional Sparse Coding[END_REF] raised a memory error with the full image on a machine with 250 GB of RAM, we used a random-patch of size 512 × 512 as input data. Due to the non-convexity of the problem, the algorithms converge to different local minima, even when initialized with the same dictionary. To deal with this difficulty, the Figure 8 shows 5 restarts of both algorithms, as well as the median curve, obtained by interpolation. We can see that DiCoDiLe outperforms the Consensus ADMM as it converges faster and to a better local minima. To ensure a fair comparison, the objective function is computed after projecting back the atoms of the dictionary D D D k to the the inverse of the dictionary norm i.e.

1 D D D k 2 2
. This projection and scaling step is necessary for the solver as ADMM iterations do not guarantee that all iterates are feasible (i.e., satisfy the constraints). This also explains the presence of several bumps in the evolution curve of the objective function.

Pattern recovery with large images

To illustrate the benefit of working with large images, we propose two pattern recovery experiments on text images and on an image of the Hubble Space Telescope.

Text images Figure 9 presents the patterns that are recovered using CDL and Online Dictionary Learning (Online DL 4 , Mairal et al. 20) from a text image composed of the letters PAMI. The text is generated by drawing uniformly 5000 characters from a set composed of these K * = 4 letters and 2 spaces. Text images X, mimicking documents obtained with a typewritter, are generated from this text using L A T E X and font beramono. They are then converted to PNG using ImageMagick [32]. We add a Gaussian white noise with standard deviation σ std(X) and σ = 3 to X (see Figure 9 bottom). The figure presents the qualitative results obtained by drawing uniformly K = 10 patches (dashed green), or learning K patterns with DiCoDiLe (orange; λ = 0.2) and Online DL (blue; λ = 0.1). The quality of the recovered atoms is evaluated using ρ, the convolutional cosine similarity following optimal assignment defined as

min σ∈S K K * l=1 ρ(d * l , d σ(l) ); ρ(d, d) = max t∈V Θ (0) W τ ∈Θ d[t+τ ] d[τ ] . ( 25
)
The atoms that are presented for each line are selected with the argmin permutation in [START_REF] Pereyra | Maximum-A-Posteriori Estimation with Unknown Regularization Parameters[END_REF]. This argmin over the set of permutations is computed using the so-called Hungarian algorithm [START_REF] Kuhn | The Hungarian method for the assignment problem[END_REF]. While DiCoDiLe is able to successfully recover the original patterns, Online DL learns blurred atoms that are unrecognizable. This highlights one of the main interest of CDL over ordinary DL, which is its ability to capture the patterns in the data in a shift invariant way. Indeed, as DL encodes all patches in a image, the recovered patterns are inevitably blurred as all phases of a pattern are of equal importance. In contrast, CDL selects specific shifts that are used to compute the atoms, allowing for sharper pattern recovery.

To assess the impact of the image size on the performance of CDL with DiCoDiLe, we generated images following the same protocol as before but drawing uniformly K * = 5 letters from a-z and then generating text with these letters for varying text-length before converting them to images X. A Gaussian white noise with standard deviation sigma * std(X) is then added to the image. Figure 10 presents for varying text sizes the recovery score ρ averaged on 10 repetition for Online DL (blue) and DiCoDiLe (orange) both with λ = 0.2. We also report the performances obtained with random patches (dashed green) extracted from the noisy input image X and with random normal patches (dashed black) obtained by sampling each pixel from a normal distribution N (.5, 1). For low level of noise σ = .1 (left), the recovery of the patterns with CDL is already very close to 1 with only small images (0.2 Mpx) with 150 characters. We see here that CDL consistently outperforms DL which have a performance a bit better than extracting random patches from the image. For large level of noise σ = 3 (right), we see that the random patches are very noisy and perform worst than normal patches that are centered around 0.5. Both Online DL and CDL with the size of the image as they take advantage of the repetition of the patterns to separate them from the noise. DiCoDiLe is able to consistently recover the shape of the original patterns, up to a level close to the one obtained in noiseless setting for very large images (4.8 Mpx) with 5000 characters. This experiment clearly demonstrates the need for CDL algorithms that are able to scale to large input data, as it is a necessary condition to design efficient pattern recovery system. We then report running times of DiCoDiLe and Online DL for this use case. In this figure, we can see that for small images, the running time of one iteration of DiCoDiLe is similar for any number of workers. This can be explained by the fact that the cost of the dictionary update step is dominating the computations. Also using a lot of workers for small images increases the number of interfering updates -thus reducing the benefit of parallelism. However, when the size of the images increases, the cost of the Z step increases and using DiCoDiLe with more than W = 1 worker reduces massively the running time -up to 10 times for larger images with 16 workers. In comparison, the benefit of using more workers is constant for Online DL. For large images where one epoch amounts to many patches, the running time of Online DL with 16 workers for one epoch becomes similar to the one of DiCoDiLe with the same number of workers. This experiment demonstrates that DiCoDiLe is able to significantly accelerate the resolution of CDL for large images by leveraging distributed computations, and can be competitive with Online DL in term of computational performances.

Hubble Telescope images

We present here results using Di-CoDiLe to learn patterns from an image of the GOODS South field, acquired by the Hubble Space Telescope [START_REF] Giavalisco | The Great Observatories Origins Deep Survey: Initial Results From Optical and Near-Infrared Imaging[END_REF]. We used W = 400 workers to learn 25 atoms of size 32 × 32 with regularization parameter λ = 0.1λmax. The learned patterns are presented in Figure 12. This unsupervised algorithm is able to highlight structured patterns in the image, such as small stars. Patterns 1 and 7 are very fuzzy. This is probably due to the presence of very large object in the foreground. As this algorithm is not scale invariant, the objects larger than our patterns are encoded using these low frequency atoms.

CONCLUSION

This work introduces DiCoDiLe, an asynchronous distributed algorithm to apply convolutional dictionary learning to very large data such as signals and images. The number of workers that can be used to solve the CSC scales linearly with the data size, allowing to adapt the computational resources to the problem dimensions. Using careful distributed pre-computations, the complexity of the dictionary update step is independent of the data size. Experiments show that the proposed algorithm outperforms the current state-of-theart parallel algorithms for CDL. Thanks to its good scaling properties, DiCoDiLe is able to process large data such as astronomical images, improving the unsupervised recovery of recurrent patterns in large datasets. Finally, recent development are pushing to refine CDL models such as the multi-layer CSC [START_REF] Sulam | Multi-Layer Convolutional Sparse Modeling: Pursuit and Dictionary Learning[END_REF][START_REF] Sulam | On Multi-Layer Basis Pursuit, Efficient Algorithms and Convolutional Neural Networks[END_REF], CSC with strided convolution [START_REF] Simon | Rethinking the CSC Model for Natural Images[END_REF] or rank-1 CDL [START_REF] Dupré La Tour | Multivariate Convolutional Sparse Coding for Electromagnetic Brain Signals[END_REF]. These models rely on a similar encoding step that can be computed efficiently with our distributed method, opening for further investigations and improved implementations.

APPENDIX A COMPUTATION FOR THE COST UPDATES

When a coordinate Z k [ω] is updated to Z k 0 [ω0], then the change in the cost function ∆E k 0 [ω0] is given by:

∆E k 0 [ω0] = D D D k 0 2 2 2 (Z k 0 [ω0] 2 -Z (q) k 0 [ω0] 2 ) -β k 0 [ω0](Z k 0 [ω0] -Z (q) k 0 [ω0]) + λ(|Z k 0 [ω0]| -| Z (q) k 0 [ω0]|).
Proof. We denote Z (1)

k [ω] = Z (q) k 0 [ω0], if (k, ω) = (k0, ω0) Z k [ω],
elsewhere .

Let α0 be the residual when coordinate (k0, ω0) of Z is set to 0.

(α k 0 [ω] -D D D k 0 [ω -ω0]Z k 0 [ω0]) 2 - 1 2 ω∈Ω α k 0 [ω] -D D D k 0 [ω -ω0] Z (q) k 0 [ω0] 2 + λ(|Z k 0 [ω0]| -| Z (q) k 0 [ω0]|) = 1 2 ω∈Ω D D D k 0 [ω -ω0] 2 (Z k 0 [ω0] 2 -Z (q) k 0 [ω0] 2 ) - ω∈Ω α k 0 [ω]D D D k 0 [ω -ω0](Z k 0 [ω0] -Z (q) k 0 [ω0]) + λ(|Z k 0 [ω0]| -| Z (q) k 0 [ω0]|) = D D D k 0 2 2 2 (Z k 0 [ω0] 2 -Z (q) k 0 [ω0] 2 ) -(α k 0 * D D D k 0 )[ω0] β k 0 [ω 0 ] (Z k 0 [ω0] -Z (q) k 0 [ω0]) + λ(|Z k 0 [ω0]| -| Z (q) k 0 [ω0]|)
This concludes our proof.

Using this result, we can derive the optimal value Z (q) k 0 [ω0] to update the coordinate (k0, ω0) as the solution of the following optimization problem: In the case where multiple coordinates (kw, ωw) are updated in the same iteration to values Z (q) kw [ωw], we obtain the following cost variation. Assuming that the updates are spread uniformly on Ω, the probability that an update candidate at random localization ω is accepted is lower bounded by

P (ω ∈ SL) ≥ d i=1 1 - WiLi Ti .
Proof. We denote ω0 the position of an update candidate, and Wω 0 = {w ∈ [1, W ] s.t. Sw ∩ VΘ(ω) = ∅} the set of workers that are impacted by this update. We denote jω 0 = |Wω 0 | the number of workers in this set. With the definition of the softlock, there exists at least one worker w in Wω 0 which has a coordinate that is not soft-locked. As the updates are uniformly spread across the workers, the probability of an update in ω0 being accepted satisfies P (ω ∈ SL | ω = ω0) ≥ 1 jω 0 .

Since one assumes uniform update distribution, the probability to update a coordinate at location ω0 is P (ω0) = 1 |Ω| . The probability for an update candidate at a random location ω to be accepted satisfies

Fig. 1 .

 1 Fig. 1. Decomposition of a noiseless univariate signal X (blue) as the convolution Z * D D D between a temporal pattern D D D (orange) and a sparse activation signal Z (green).

Algorithm 2

 2 DiCoDiLe with W workersRequire: D D D (0) , X, stopping criterion ν.1: repeat 2:

Definition 3 . 1

 31 (Θ-Border). Let us consider a sub-domain, Sw = d i=1 [li, ui[, where (li, ui) are the lower and upper limits of this Algorithm 3 DiCoDiLeZ with W workers 1: Input: D D D, X, parameter > 0, sub-domains Sw, 2: In parallel for w = 1 • • • W 3: Compute a partition {C (w) m } M m=1 of the worker domain Sw with sub-domains of size 2 d L. 4: Initialize β k [ω] and Z k [ω] ∀(k, ω) ∈ [1, K] × Sw, 5: repeat 6: for m=1. . . M do 7:

Fig. 3 .

 3 Fig. 3. Average running time of CD with three coordinate selection schemes -(green) Cyclic, (orange) Locally Greedy and (blue) Greedy -for three signal lengths and λ = 0.5λmax.LGCD consistently outperforms the two other strategies.

Fig. 4 .

 4 Fig. 4. Reconstruction of the image Mandrill using DiCoDiLe Z (left) without and (right) with soft-locks for a grid of 7 × 7 workers. The dashed lines show the partitions of the domain Ω. The algorithm diverges at the edges of some of the sub-domains due to interfering updates between more than two workers.

Fig. 6 .

 6 Fig.6. Scaling of DICOD (orange) and DiCoDiLe Z (blue) with the number of workers W . DiCoDiLe Z scales sub-linearly while DICOD scales super-linearly. However, DiCoDiLe Z is more efficient than DICOD in a regime with a low number of workers W , leading to improved performances in all regimes. The green line denotes the number of cores where DiCoDiLe Z and DICOD become the same as DiCoDiLe Z only has one sub-domain in each worker. The results are consistent for both small (top, T = 150L) and large (bottom T = 750L) signals sizes.

4 .Fig. 9 .

 49 Fig.9. Letter recovery experiment using a noisy image of a text with 5000 characters drawn uniformly from PAMI (top). For each row, the atoms that are displayed have been selected among K = 10 with the best assignment ρ with the original patterns (dashed red). Random patches drawn uniformly from the original image (dashed green) are used as initial estimates D 0 . We use DiCoDiLe (orange) and Online DL (blue,[START_REF] Mairal | Online Learning for Matrix Factorization and Sparse Coding[END_REF]) to learn K = 10 atoms. While the atoms recovered with Online DL are blurred, CDL is able to correcly recover the 4 original letters.

Fig. 10 .Fig. 11 .

 1011 Fig.10. Effect of the image size (in Megapixels) on letter recovery for two noise levels. The quality of the recovered patterns is evaluated with best assignment ρ. One observes that increasing image size improves recovery of the patterns which are here letters, as the number of characters present in the image increases.

Figure 11 displays

 11 Figure 11 displays the running time per iteration for DiCoDiLe and per epoch for Online DL with W in {1, 4, 16} for the text images used in the previous experiment. We consider data with σ = 3 and different image sizes. The reported running times are the median time over 10 repetitions of the experiment.In this figure, we can see that for small images, the running time of one iteration of DiCoDiLe is similar for any number of workers. This can be explained by the fact that the cost of the dictionary update step is dominating the computations. Also using a lot of workers for small images increases the number of interfering updates -thus reducing the benefit of parallelism. However, when the size of the images increases, the cost of the Z step increases and using DiCoDiLe with more than W = 1 worker reduces massively the running time -up to 10 times for larger images with 16 workers. In comparison, the benefit of using more workers is constant for Online DL. For large images where one epoch amounts to many patches, the running time of Online DL with 16 workers for one epoch becomes similar to the one of DiCoDiLe with the same number of workers. This

Fig. 12 .

 12 Fig. 12. 25 atoms 32 × 32 learned from the Hubble Space Telescope images (STScI-H-2016-39-a) using DiCoDiLe with 400 workers and a regularization parameter λ = 0.1λmax. The atoms are sorted based on Z k 1 from the top-left to the bottom-right. Most atoms displays a structure which corresponds to spatial objects. The atoms 1 and 7 have fuzzier structure because they are mainly used to encode larger scale objects.

k

  [ω], the cost update is a simple function ofZ k [ω] and Z (q) k [ω]. Proposition A.1. If the coordinate (k0, ω0) is updated from value Z k 0 [ω0] to value Z (q)

1 - 1 2 ω∈ΩX 2 [

 1122 For ω ∈ Ω,α k 0 [ω] = (X -Z * D D D)[ω] + D D D k 0 [ω -ω0]Z k 0 [ω0] = (X -Z * D D D)[ω] + Z k 0 [ω0](eω 0 * D D D k 0 )[ω0] .We also haveα k 0 [ω] = (X -Z (1) * D D D)[ω] + D D D k 0 [ω -ω0] Z Z * D D D) 2 [ω] + λ Z -Z (1) * D D D

  β k 0 [ω0]) 2 + λ|u| .

Proposition A. 2 . 2 + λ Z 1 ( 2 -(αProposition B. 1 .

 22121 The update of the W coordinates (kw, ωw) W w=1 with additive update ∆Z kw[ωw] changes the cost by:D kw * D D D k w )[ω w -ωw]∆Z kw [ωw]∆Z k w [ω w ] interference , Proof. We define Z (1) k [ω] = Z (q) kw [ωw], if (k, ω) = (kw, ωw) Z k [ω],otherwise .Letγ[ω] =(X -Z * D D D)[ω] + K w=1 D D D kw [ω -ωw]Z kw [ωw] (A.2) =(X -Z (1) * D D D)[ω] + K w=1 D D D kw [ω -ωw] Z (q) kw [ωw] , (A.3)andαw[ω] =(X -Z * D D D)[ω] + D D D kw [ω -ωw]Z kw [ωw] (A.4) =γ[ω]w =w D D D k w [ω -ω w ]Z k w [ω w ] . ] -Z * D D D[ω] ] -Z (1) * D D D[ω] kw [ω -ωw]Z kw [ωw] kw [ω -ωw] 2 (Z kw [ωw] 2 -Z |Z kw [ωw]| -| Z (q) kw [ωw]|) -W w=1 ω∈Ω γ[ω]D D D kw [ω -ωw] ∆Z kw [ωw] (A.8) w =w ω∈Ω D D D kw [ω -ωw]D D D k w [ω -ω w ] × (Z kw [ωw]Z k w [ω w ] -Z kw [ω]D D D kw [ω -ωw] α kw * D D D kw [ωw ] ∆Z kw [ωw] + λ(|Z kw [ωw]| -| Z D k w * D D D kw )[ωw -ω w ] × ∆Z kw [ωw] Z (q) k w [ω w ] + ∆Z k w [ω w ] Z (q) kw [ωw] -Z kw [ωw]Z k w [ω w ] -Z D k w * D D D kw )[ωw -ω w ]∆Z kw [ωw]∆Z k w [ω w ]APPENDIX B PROBABILITY OF A COORDINATE TO BE SOFT-LOCKEDWe consider DiCoDiLeZ with W workers (Wi in each direction), run on a signal with d dimensions. We denote Ti the length of the signal and Li the length of the atoms in each dimension i.

P

  (ω ∈ SL) = W w=1 ω 0 ∈Sw P (ω ∈ SL | ω = ω0)P (ω0)

  Our distributed algorithm DiCoDiLe is introduced in Section 3. Section 4 presents results on simulations, signals and images, demonstrating improvements in both scalability and pattern recovery.

Notation For a vector U ∈ R P we denote Up ∈ R its p th coordinate. For a finite d-dimensional domain Ω = d i=1 [0, Ti[, where [0, Ti[ are the integers from 0

  lock area Messages Fig.2. Communication process in DiCoDiLe Z with d = 2 and 9 workers centered around worker w. The update in ω 0 is independent of the other workers. The update in ω 1 is performed if no better coordinate update is possible in the soft-lock area V Θ (ω 1 ) ∩ S w+1 (red hatched area). If accepted, the worker w + 1 needs to be notified. The update in ω 2 changes the value of the optimal updates in the Θ-extension of the other workers sub-domains. Thus it needs to notify all the neighboring workers.

	iterative steps
	W
	∆E =	∆E kw [ωw]
	i=1

derivations, detailed in Proposition A.2, show that

  close to W and DiCoDiLeZ scales almost linearly. When Wi reaches T i L i ( 2 1/d 2 1/d -1 ), then P (ω ∈ SL)

	1 2 ,
	and the acceleration is reduced to W 2 . Numerical experiments
	in Section 4 will evaluate if this result holds in practice, al-
	though the updates are not distributed uniformly when run-
	ning LGCD on real data.

1.

In this case, the expected number of accepted update candidates is

  2-ball by scaling it with D D D k 2 2 . This implies that we also needed to scale Z by 3. Code available at https://sporco.readthedocs.io/

	10 4			Skau et al. (2018)
	Loss F(Z; D)			DiCoDiLe
	10 3	0	500	1000 1500 2000 2500 3000
				Time [sec]
	Fig. 8. Comparison of DiCoDiLe algorithm with Consensus ADMM
	dictionary learning proposed by Skau and Wohlberg [28]. Algorithms
	were run 5 times with 36 workers on a of size 512 × 512
	taken from the Hubble Space Telescope GOODS South image. The
	solid lines denote the median curves obtained through interpolation.

ACKNOWLEDGMENTS

This work was supported by the ERC Starting Grant SLAB ERC-StG-676943.

where (B.3) results from the summation on Ω. Basically, in each worker w, there is d i=1 (Ti/Wi -2Li) coordinates for which jω 0 = 1. This corresponds to the coordinate in the worker that are not on the border i.e. Sw \ BΘ(Sw). For a coordinate at the border with jω 0 = k > 1, it suffices to notice that it will appear k times in the sum due to symmetry considerations. This is illustrated in Figure 13.


Farah Abdel Samad

Frederic Medail

Alex Baumel

Sonja Siljak Yakovlev

Anthony Roukoz

Georges Tohme

Henriette Tohme

Magda Bou

Dagher Kharrat

New insights on the conservation status of the Endangered coastal endemic plant Astragalus berytheus (Fabaceae) in Lebanon

Keywords: Astragalus berytheus, coastal ecosystem, conservation, Critically Endangered, Eastern Mediterranean, endemic, extinction, habitat degradation

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Supplementary material for this article is available at doi.org/./S C oastal vegetation includes unique species with adap- tations to the specific conditions of this environment (Voronkova et al., ). In the Mediterranean region urban development of coastal areas has severely affected coastal ecosystems and led to loss of biodiversity (Doxa et al., ). In Lebanon, one of the  regional biodiversity hotspots in this region (Médail & Myers, ), the increase in the human population and intensive urbanization has destroyed one-third of the  km coastline (MOE, ).

Astragalus berytheus Boiss. & Blanche (Plate ) is a coastal species endemic to Lebanon and Palestine/Israel (Tohmé & Tohmé, ), and one of the most threatened taxa in Lebanon (MOE, ). It occurs in prime tourist areas, such as sandy beaches, which are targeted for the establishment of resorts. These habitats are categorized as highly threatened because of extensive tourism throughout the Mediterranean (IUCN, ). Here we assess the conservation status of A. berytheus, using the IUCN Red List criteria (IUCN, ), and make appropriate recommendations for its conservation.

Information on the historical distribution of A. berytheus was gathered from regional floras (Mouterde, ; Tohmé & Tohmé, ) and herbaria (Royal Botanic Garden Edinburgh Herbarium, Herbarium WU at University of Vienna, and MNHN Paris). During -, all sandy shores along the coast of Lebanon were visited by GT and HT to determine the species' occurrence. Historically the species was known from five localities in Lebanon (Fig. , Supplementary Table ).

Three of the five Lebanese populations (Beirut, Ouzaï and Khaldeh) were already extinct when we began this study in . The coastal strip was severely damaged by the construction of the highway connecting Beirut to southern Lebanon and the construction of Beirut airport in . The last population in the Beirut area (Bir Hassan) became extinct following construction of a residential building in . This population was in an area of , m  and was characterized by annual fluctuation in the number of individuals: eight in ,  in  and five in .

The only extant population of A. berytheus in Lebanon is in an area of c. , m  in Tyre Coast Nature Reserve. In  there were .  individuals, bearing healthy fruits and seeds. These were, however, affected by the invasion of camphor-weed Heterotheca subaxillaris (Lam.) Britton & Rusby We calculated the AOO (area of occupancy) and the EOO (extent of occurrence) of A. berytheus using the Geospatial Conservation Assessment Tool (Bachman et al., ), considering the Lebanese population and the populations in Palestine/Israel (Danin & Fragman-Sapir, ; Supplementary Table ). With an EOO of c.  km  and an AOO of c.  km  , A. berytheus should be categorized as Endangered based on criteria Bab(ii); i.e. with AOO ,  km  (B), severely fragmented (a) and projected decline of quality of habitat b(ii).

Nationally in Lebanon, A. berytheus should be categorized as Critically Endangered based on criteria Aac; BBab(i, ii, iii, iv); Cb; i.e. with a population reduction observed of . % and where the causes of reduction may not be reversible (A) and direct observation (a) of decline in AOO and EOO (c), with an EOO ,  km² (B) and an AOO ,  km² (B) showing only one location (a) and a continuing decline of EOO (i), AOO (ii), quality of habitat (iii), and number of subpopulations (iv), and with ,  mature individuals in the wild (C) and extreme fluctuations in the number of mature individuals (b).

Passive habitat protection is not a guarantee of longterm plant conservation (Martinell et al., ), and the last A. berytheus population in Lebanon would probably have been extirpated if encroachment of the invasive H. subaxillarisi had not been halted. Plant species known from only one or two populations, growing in areas with some protection, are known to have gone extinct (Bañares et al., ). Appropriate planning of local management actions, with recovery plans for individual species (Fenu et al., ; Ballesteros et al., ; Cambrollé et al., ), are required for successful plant conservation. Environmental monitoring, population viability analysis and population genetics and studies of reproductive biology are required to support the conservation of A. berytheus. The microreserve concept (Laguna et al., ), for the monitoring and in situ conservation of threatened plants in small protected areas, could be appropriate for A. berytheus and for other narrow endemic species of Lebanon, such as Matthiola crassifolia Boiss. & Gaill. (Bou Dagher Kharrat et al., ).

In ,  pods of A. berytheus were collected from the Bir Hassan population (before it vanished), containing a total of  seeds, but only nine seeds escaped infestation by a chalcid wasp Bruchophagus sp. (Hymenoptera: Eurytomidae) often associated with the seeds of legumes (Burks, ). Theses seeds, together with some seeds collected from the Tyre Coast Nature Reserve in  were stored in the Laboratory of Seed Germination and Conservation-Jouzour Loubnan at Saint-Joseph University, Beirut, Lebanon.

In ,  plantlets grown from ex situ germination of these seeds will be introduced to the Tyre archaeological site,  km from Tyre Nature Reserve, with similar ecological conditions. As the archaeological site is a UNESCO world heritage site this will facilitate protection and monitoring. The presence of endemic plants on archaeological sites is an opportunity to raise the conservation awareness of visitors (Siljak-Yakovlev et al., ).

A national legal framework is required to halt the loss of biodiversity in Lebanon, especially of narrow endemics plant species (Médail & Baumel, ), working from the principle that the narrower the range of a taxon the greater the local responsibility (Gauthier et al., ). An additional imperative to conserve A. berytheus is that it bears the name of Beirut, the nation's capital. If the results of the pilot introduction of the species to the Tyre archaeological site are favourable, future introductions could be considered in other suitable coastal locations.

(

  Family Asteraceae), which originated from the USA; it is an invader of coastal sand dunes (Sternberg, ) and was introduced into the Middle East in  (di Castri et al., ). It was first observed in the Reserve in  (Tohmé & Tohmé, ), and in  it invaded A. berytheus habitat. Eradication has, however, allowed A. berytheus to recover (Bou Dagher Kharrat et al., ), with  individuals recorded in  and  individuals bearing fruits in .

PLATE 1

 1 FIG. 1 Present and historical distribution of Astragalus berytheus on the eastern Mediterranean coast (see Supplementary Tablefor additional details). The only extant population in Lebanon is in Tyre Coast Nature Reserve (see text for details).

  FIG. 1 Present and historical distribution of Astragalus berytheus on the eastern Mediterranean coast (see Supplementary Tablefor additional details). The only extant population in Lebanon is in Tyre Coast Nature Reserve (see text for details).

  FIG. 1 Present and historical distribution of Astragalus berytheus on the eastern Mediterranean coast (see Supplementary Tablefor additional details). The only extant population in Lebanon is in Tyre Coast Nature Reserve (see text for details).

Acknowledgements This research was supported by a PhD grant from the National Council of Scientific Research of Lebanon and from Critical Ecosystem Partnership Fund grant #108497.

Author contributions Fieldwork: FAS, GT, HT; writing: FAS, FM, MBDK, AB, SSY.

Conflicts of interest None.