Tung Duy Luu
email: duy-tung.luu@ensicaen.fr.

Jalal Fadili
email: jalal.fadili@ensicaen.fr.

Christophe Chesneau
email: christophe.chesneau@unicaen.fr.

Sampling from non-smooth distributions through Langevin diffusion

In this paper, we propose proximal splitting-type algorithms for sampling from distributions whose densities are not necessarily smooth nor log-concave. Our approach brings together tools from, on the one hand, variational analysis and non-smooth optimization, and on the other hand, stochastic diffusion equations, and in particular the Langevin diffusion. We establish in particular consistency guarantees of our algorithms seen as discretization schemes in this context. These algorithms are then applied to compute the exponentially weighted aggregates for regression problems involving non-smooth priors encouraging some notion of simplicity/complexity. Some popular priors are detailed and implemented on some numerical experiments.

Introduction

Problem statement

We consider the linear regression problem

y = Xθ 0 + ξ, (1.1)
where y ∈ R n is the vector of observations, X ∈ R n×p is the design matrix, ξ is the vector of errors, and θ 0 ∈ R p is the unknown regression vector we wish to estimate. X ∈ R n×p can be seen as the sensing or degradation operator in inverse problems raising in, e.g., signal and image processing, or the design matrix for a regression problem in statistics and machine learning. Generally, problem (1.1) is either under-determined (p < n), or determined (p = n) but X is ill-conditioned. In both cases, (1.1) is ill-posed. The idea of aggregating elements in a dictionary has been introduced in machine learning to combine different techniques (see [START_REF] Littlestone | The weighted majority algorithm[END_REF][START_REF] Vovk | Aggregating strategies[END_REF]) with some procedures such as bagging [START_REF] Breiman | Bagging predictors[END_REF], boosting [START_REF] Freund | Boosting a weak learning algorithm by majority[END_REF][START_REF] Schapire | The strength of weak learnability[END_REF] and random forests [1, 8-10, 13, 36]. In the recent years, there has been a flurry of research on the use of low-complexity regularization (among which sparsity and low-rank are the most popular) in various areas including statistics and machine learning in high dimension. The idea is that θ 0 generally conforms to some notions of sparsity/low-complexity. Namely, it has either a simple structure or a small intrinsic dimension. This makes it possible to build an estimate X θ with good provable performance guarantees under appropriate conditions. In literature, the information of sparsity/low-complexity has been taken into account through two families of estimators: Penalized Estimators and Exponentially Weighted Aggregates (EWA).

Variational/Penalized Estimators

The penalized approach consists in imposing on the set of candidate solutions some prior structure on the object to be estimated. The class of estimators are obtained by solving the convex optimization problem

θ PEN n ∈ Argmin θ∈R p V (θ) def = F (Xθ, y) + J λ (θ) , (1.2)
where F : R n × R n → R is a general loss function assumed to be differentiable, J λ : R p → R is the regularizing penalty promoting some specific notion of simplicity/low-complexity which depends on a vector of parameters λ. Regularization is now a central theme in many fields including statistics, machine learning and inverse problems. A prominent member covered by (1.2) is the Lasso [START_REF] Bickel | Simultaneous analysis of lasso and Dantzig selector[END_REF][START_REF] Bühlmann | Statistics for High-Dimensional Data: Methods, Theory and Applications[END_REF][START_REF] Candès | Near-ideal model selection by 1 minimization[END_REF][START_REF] Chen | Atomic decomposition by basis pursuit[END_REF][START_REF] Donoho | For most large underdetermined systems of linear equations the minimal 1 -norm solution is also the sparsest solution[END_REF][START_REF] Osborne | A new approach to variable selection in least squares problems[END_REF][START_REF] Tibshirani | Regression shrinkage and selection via the Lasso[END_REF] and its variants such the analysis/fused Lasso [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF][START_REF] Tibshirani | Sparsity and smoothness via the fused Lasso[END_REF] or group Lasso [START_REF] Bach | Consistency of the group lasso and multiple kernel learning[END_REF][START_REF] Bakin | Adaptive regression and model selection in data mining problems[END_REF][START_REF] Chesneau | Some theoretical results on the grouped variables lasso[END_REF][START_REF] Wei | Consistent group selection in high-dimensional linear regression[END_REF][START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF]. Another example is the nuclear norm minimization for low rank matrix recovery motivated by various applications including robust PCA, phase retrieval, control and computer vision [START_REF] Candès | Exact matrix completion via convex optimization[END_REF][START_REF] Candès | Phaselift: Exact and stable signal recovery from magnitude measurements via convex programming[END_REF][START_REF] Fazel | A rank minimization heuristic with application to minimum order system approximation[END_REF][START_REF] Recht | Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization[END_REF]. See [START_REF] Bühlmann | Statistics for High-Dimensional Data: Methods, Theory and Applications[END_REF][START_REF] Negahban | A unified framework for high-dimensional analysis of M-estimators with decomposable regularizers[END_REF][START_REF] Vaiter | Low complexity regularization of linear inverse problems[END_REF][START_REF] Van De Geer | Weakly decomposable regularization penalties and structured sparsity[END_REF] for generalizations and comprehensive reviews.

Exponential Weighted Aggregation (EWA)

An alternative to the variational estimator (1.2) is the aggregation by exponential weighting which combines all of candidate solutions with the aggregators promoting the prior information. The aggregators are defined via the probability density function

µ(θ) = exp (-V (θ)/β) Θ exp (-V (ω)/β)dω , (1.3)
where β > 0 is called temperature parameter. If all θ are candidates to estimate the true vector θ 0 , then Θ = R p . The aggregate is thus defined by

θ EWA n = R p
θ µ(θ)dθ.

(1.4)

Aggregation by exponential weighting has been widely considered in the statistical and machine learning literatures, see e.g. [START_REF] Dalalyan | Pac-bayesian bounds for the expected error of aggregation by exponential weights[END_REF][START_REF] Dalalyan | Aggregation by exponential weighting, sharp pac-bayesian bounds and sparsity[END_REF][START_REF] Dalalyan | Aggregation by exponential weighting and sharp oracle inequalities[END_REF][START_REF] Dalalyan | Sparse regression learning by aggregation and langevin montecarlo[END_REF][START_REF] Luu | PAC-Bayesian risk bounds for group-analysis sparse regression by exponential weighting[END_REF][START_REF] Guedj | Pac-bayesian estimation and prediction in sparse additive models[END_REF][START_REF] Lecué | Simultaneous adaptation to the margin and to complexity in classification[END_REF][START_REF] Nemirovski | Topics in non-parametric statistics[END_REF][START_REF] Rigollet | Linear and convex aggregation of density estimators[END_REF][START_REF] Yang | Aggregating regression procedures to improve performance[END_REF] to name a few.

The Langevin diffusion

In this paper, we focus on the computation of EWA. Computing θ EWA n in (1.4) corresponds to an integration problem which becomes very involved to solve analytically or even numerically in high-dimension. A classical alternative is to approximate it via a Markov chain Monte-Carlo (MCMC) method which consists in sampling from µ by constructing an appropriate Markov chain whose stationary distribution is µ, and to compute sample path averages based on the output of the Markov chain. The theory of MCMC methods is based on that of Markov chains on continuous state space. As in [START_REF] Dalalyan | Sparse regression learning by aggregation and langevin montecarlo[END_REF], we here use the Langevin diffusion process; see [START_REF] Roberts | Exponential Convergence of Langevin Distributions and Their Discrete Approximations[END_REF].

Continuous dynamics

A Langevin diffusion L in R p , p ≥ 1 is a homogeneous Markov process defined by the stochastic differential equation (SDE)

dL(t) = 1 2 ρ(L(t))dt + dW (t), t > 0, L(0) = l 0 , (1.5)
where ρ = ∇ log µ, µ is everywhere non-zero and suitably smooth target density function on R p , W is a p-dimensional Brownian process and l 0 ∈ R p is the initial value. Under mild assumptions, the SDE (1.5) has a unique strong solution and, L(t) has a stationary distribution with density precisely µ [54, Theorem 2.1]. L(t) is therefore interesting for sampling from µ. In particular, this opens the door to approximating integrals R p f (θ)µ(θ)dθ, where f : R p → R, by the average value of a Langevin diffusion, i.e.,1 T T 0 f (L(t))dt for a large enough T . Under additional assumptions on µ, the expected squared error of the approximation can be controlled [START_REF] Xuerong | Stochastic differential equations and applications[END_REF].

Forward Euler discretization

In practice, in simulating the diffusion sample path, we cannot follow exactly the dynamic defined by the SDE (1.5). Instead, we must discretize it. A popular discretization is given by the forward (Euler) scheme, which reads

L k+1 = L k + δ 2 ρ(L k ) + √ δZ k , t > 0, L 0 = l 0 ,
where δ > 0 is a sufficiently small constant discretization step-size and {Z k } k are iid ∼ N (0, I p ). The average value 1 T T 0 L(t)dt can then be naturally approximated via the Riemann sum

δ T T /δ -1 k=0 L k , (1.6)
where T /δ denotes the interger part of T /δ. It is then natural to approximate θ by applying this discretization strategy to the Langevin diffusion with µ as the target density. However, quantitative consistency guarantees of this discretization require µ (hence ρ) to be sufficiently smooth. For a comprehensive review of sampling by Langevin diffusion from smooth and log-concave densities, we refer the reader to e.g. [START_REF] Dalalyan | Theoretical guarantees for approximate sampling from a smooth and log-concave density[END_REF]. To cope with non-smooth densities, several works have proposed to replace log µ with a smoothed version (typically involving the Moreau-Yosida regularization/envelope, see Definition 2.2) [START_REF] Dalalyan | Sparse regression learning by aggregation and langevin montecarlo[END_REF][START_REF] Durmus | Non-asymptotic convergence analysis for the Unadjusted Langevin Algorithm[END_REF][START_REF] Durmus | Sampling from convex non continuously differentiable functions, when Moreau meets Langevin[END_REF][START_REF] Pereyra | Proximal markov chain monte carlo algorithms[END_REF]. In [START_REF] Durmus | Sampling from convex non continuously differentiable functions, when Moreau meets Langevin[END_REF][START_REF] Pereyra | Proximal markov chain monte carlo algorithms[END_REF] for instance, the authors proposed proximal-type algorithms to sample from possibly non-smooth log-concave densities µ using the forward Euler discretization and the Moreau-Yosida regularization. In [START_REF] Pereyra | Proximal markov chain monte carlo algorithms[END_REF] 1 , -log µ is replaced with its Moreau envelope, while in [START_REF] Durmus | Sampling from convex non continuously differentiable functions, when Moreau meets Langevin[END_REF], it is assumed that -log µ = L + H, L is convex Lipschitz continuously differentiable, and H is a proper closed convex function replaced by its Moreau envelope. In both these works, convexity plays a crucial role to get quantitative convergence guarantees. Proximal steps within MCMC methods have been recently proposed for some simple (convex) signal processing problems [START_REF] Chaari | A hamiltonian monte carlo method for non-smooth energy sampling[END_REF], though without any guarantees.

Contributions

Our main contributions are summarized as follows.

• We aim to enlarge the family of µ covered by [START_REF] Dalalyan | Sparse regression learning by aggregation and langevin montecarlo[END_REF][START_REF] Durmus | Non-asymptotic convergence analysis for the Unadjusted Langevin Algorithm[END_REF][START_REF] Durmus | Sampling from convex non continuously differentiable functions, when Moreau meets Langevin[END_REF][START_REF] Pereyra | Proximal markov chain monte carlo algorithms[END_REF] by relaxing some underlying conditions. Especially, in our study, µ is structured as µ in (1.3), and it is not necessarily differentiable nor logconcave.

• We propose two algorithms based on forward-backward proximal splitting for which we prove theoretical consistency guarantees.

• These algorithms are applied to compute EWA estimators with several popular penalties in the literature, and illustrated on some numerical problems.

Paper organization

Some preliminaries, definitions and notations are introduced in Section 2. Section 3 establishes key properties of a Moreau-Yosida regularized version of µ under mild assumptions of the latter. In turn we will consider the SDE (1.5) with such the smoothed density. Well-posedness of this SDE and consistency guarantees for its discrete approximations are proven in Section 4. Section 5 provides a large class of functions, namely prox-regular functions, for which the previous theoretical analysis applies. From this analysis, two algorithms are derived in Section 6 and applied in Section 7 to compute the EWA estimator with several penalties. The numerical experiments are described in Section 8. The proofs of all results are collected in Section 9.

Notations and Preliminaries

Before proceeding, let us introduce some notations and definitions.

Vectors and matrices F or a d-dimensional Euclidean space R d , we endow it with its usual inner product •, • and associated norm • 2 . I d is the identity matrix on R d . For r ≥ 1, • r will denote the r norm of a vector with the usual adaptation for r = +∞. Let M ∈ R d×d symmetric positive definite, we denote •, • M = •, M • and • M its associated norm. For a matrix M , we denote σ min (M ) its smallest singular value and |||M ||| its spectral norm. Of course,

• M and • 2 are equivalent. Let x ∈ R d and the subset of indices I ⊂ {1, . . . , d}. We denote x I the subvector whose entries are those of x indexed by I. For any matrix M , M denotes its transpose.

Sets For a set C, denote I C its characteristic function, i.e., 1 if the argument is in C and 0 otherwise, and ι C its its indicator function, i.e., 0 if the argument is in C and +∞ otherwise. For an index set I, I is its cardinality.

Functions

We will denote (•) + = max(•, 0) the positive part of a real number. For a function f :

R d → R∪ {-∞, +∞}, its effective domain is dom(f ) = x ∈ R d : f (x) < +∞ and f is proper if f (x) > -∞ for all x and dom(f ) = ∅ as is the case when it is finite-valued. A function f : R d → R ∪ {-∞, +∞} is lower semi continuous (lsc) at x 0 if lim inf x→x 0 f (x) ≥ f (x 0 ). For a differentiable function f , ∇f is its (Euclidean) gradient. Define C 1,+ (R d ) (resp. C 1,1 (R d ))
the set of differentiable functions in R d whose gradient is locally (resp. globally) Lipschitz continuous. We also define

C 1,+ (R d ) def = f ∈ C 1,+ (R d ) : ∃K > 0, ∀x ∈ R d , x, ∇f (x) ≤ K(1 + x 2 2 ) . The following lemma shows that C 1,1 (R d ) ⊂ C 1,+ (R d ). Lemma 2.1. Assume that f : R d → R d is Lipschitz continuous, then there exists K > 0 such that f (x), x ≤ K(1 + x 2 2 ), ∀x ∈ R d .
Let us also consider some definitions and properties of variational analysis. A more comprehensive account on variational analysis in finite-dimensional Eudlidean spaces can be found in [START_REF] Rockafellar | Variational analysis[END_REF]. Definition 2.1 (Subdifferential). Given a point x ∈ R d where a function f : R d → R ∪ {+∞} is finite, the subdifferential of f at x is defined as

∂f (x) = v ∈ R d : ∃x k → x, f (x k ) → f (x), v ← v k ∈ ∂ F f (x k ) , where the Fréchet subdifferential ∂ F f (x) of f at x, is the set of vectors v such that f (w) ≥ f (x) + v, w -x + o w -x 2 .
We say that f is subdifferentially regular at x if and only if f is locally lsc there with ∂f (x) = ∂ F f (x).

Let us note that ∂f (x) and

∂ F f (x) are closed, with ∂ F f (x) convex and ∂ F f (x) ⊂ ∂f (x) [55, The- orem 8.6]. In particular, if f is a proper lsc convex function, ∂ F f (x) = ∂f (x)
prox M γf (x) def = Argmin w∈R d 1 2γ w -x 2 M + f (w) , M ,γ f (x) def = inf w∈R d 1 2γ w -x 2 M + f (w) ,
prox M γf here is a set-valued operator since the minimizer, if it exists, is not necessarily unique. When M = I p , we simply write prox γf and γ f . Operators For a set-valued operator S : R d ⇒ R d , its graph is gph(S) = (x, v) : v ∈ S(x) .

Definition 2.3 (Hypomonotone and monotone operators). A set-valued operator

S : R d ⇒ R d is hypomono- tone of modulus r > 0 if x -x, v -v ≥ -r x -x 2 2 , ∀(x, v) ∈ gph(S), (x , v ) ∈ gph(S).
It is monotone if the inequality holds with r = 0.

Moreau-Yosida regularization

In our framework, the target distribution µ is defined as

µ(θ) ∝ exp -L(θ) + H • D (θ) , (3.1)
where L ∈ C 1,+ (R p ), D ∈ R p×q and H : R q → R. Moreover, H is assumed neither differentiable nor convex. To overcome these difficulties, we invoke arguments from variational analysis [START_REF] Rockafellar | Variational analysis[END_REF]. Namely, we will replace H by its Moreau envelope and state the following assumptions to exploit some key properties of the latter. To avoid trivialities, from now on, we assume that Argmin(H) = ∅.

(H.1) H : R q → R is lsc and bounded from below.

(H.2) prox M γH is single valued.

Let us start with some key properties of the Moreau envelope.

Lemma 3.1. Let M ∈ R q×q depending on γ ∈]0, γ 0 [ with γ 0 > 0, we denote it M γ , such that M γ is symmetric positive definite for any γ ∈]0, γ 0 [, and γ → θ M γ , ∀θ ∈ R q , is a decreasing mapping on ]0, γ 0 [. Assume that (H.1) holds. (i) prox M γ
γH (x) are non-empty compact sets for any x, and

x ∈ Argmin(H) ⇒ x ∈ prox M γ γH (x).
(ii) M γ ,γ H(θ) is finite and depends continuously on

(x, γ) ∈ R q ×]0, γ 0 [, and M γ ,γ H(x) γ∈]0,γ 0 [ is a decreasing net. More precisely, M γ ,γ H(x) H(x) for all x as γ 0.
The fixed points of this proximal mapping include minimizers of H. They are not equal however in general, unless for instance H is convex.

Lemma 3.2. Let M γ ∈ R q×q symmetric positive definite, assume that (H.1) and (H.2) hold. Then prox M γ γH is continuous on (x, γ) ∈ R q ×]0, γ 0 [, and M γ ,γ H ∈ C 1 (R q ) with gradient ∇ M γ ,γ H = γ -1 M γ I q -prox M γ γH .
In plain words, Lemma 3.2 tells us that under (H.1)-(H.2), the Moreau envelope is a smooth function, hence the name Moreau-Yosida regularization. Moreover, the action of the operator prox M γ γH is equivalent to a gradient descent on the Moreau envelope of H in the metric M γ with step-size γ. Remark 3.1. When the metric matrix does not depend on γ, Lemmas 3.1 and 3.2 hold with γ 0 = +∞.

Let us now define the smoothed density

µ γ (θ) = exp -L(θ) + ( M ,γ H) • D (θ) Z γ , (3.2)
where

Z γ = R p exp -L(θ ) + ( M ,γ H) • D (θ ) dθ .
The following proposition answers the natural question on the behaviour of µ γ -µ as a function of γ.

Proposition 3.1. Assume that (H.1) holds. Then, µ γ converges to µ in total variation as γ → 0.

Langevin diffusion with Moreau-Yosida regularization

Let us define the following SDE with the Moreau-Yosida regularized version of H dL(t) = ψ(L(t))dt + dW (t), t > 0,

where ψ : θ ∈ R p → - 1 2 ∇ L + ( M ,γ H) • D (θ), (4.1)
ψ is the drift coefficient.

Recall that (H.1) and (H.2) were mild assumptions required to establish key properties of Moreau-Yosida regularization, which in turn allow to compute ∇ M ,γ H by exploiting its the relation between ∇ M ,γ H and prox M γH as stated in Lemma 3.2. Now, to guarantee well-posed (existence and uniqueness) and discretization consistency of the SDE (4.1), we will also need the following assumptions.

(H.3) prox M γH is locally Lipschitz continuous. (H.4) There exists C > 0 such that D θ, prox M γH (D θ) M ≤ C(1 + θ 2 2 ), ∀θ ∈ R p .

Well-posedness

We start with the following characterization of the drift ψ.

ψ(θ), θ ≤ K(1 + θ 2 2 ), for some K > 0, (4.2

Discretization

Approach 1 Inserting the identities of Lemma 3.2 into (4.1), we get the SDE

dL(t) = - 1 2 ∇L + γ -1 DM I q -prox M γH • D (L(t))dt + dW (t), L(0) = l 0 , t > 0. (4.4)
Consider now the forward Euler discretization of (4.4) with step-size δ > 0, which can be rearranged as

L k+1 = L k - δ 2 ∇L(L k ) - δ 2γ DM D L k -prox M γH (D L k ) + √ δZ k , t > 0, L 0 = l 0 . (4.5)
Note that by Lemma 3.2, and without the stochastic term √ δZ k , (4.5) amounts to a relaxed form of gradient descent on L and the Moreau envelope of H in the metric M with step-size δ.

From (4.5), an Euler approximate solution is defined as

L δ (t) def = L 0 - 1 2 t 0 ∇L(L(s)) -γ -1 DM D L(s) -prox M γH (D L(s)) ds + t 0 dW (s),
where

L(t) = L k for t ∈ [kδ, (k + 1)δ[. Observe that L δ (kδ) = L(kδ) = L k , hence L δ (t)
[ L(0) r 2 ] < ∞ for any r ≥ 2. Then E L δ (T ) -E [L(T )] 2 ≤ E sup 0≤t≤T L δ (t) -L(t) 2 -→ δ→0 0. (4.6)
The convergence rate is of order δ 1/2 when prox M γH is globally Lipschitz continuous.

Approach 2 Assume now that the metric also depends on γ ∈ (]0, γ 0 [ with γ 0 > 0, and we emphasize this by denoting it M γ such that M γ is symmetric positive definite for any γ ∈]0, γ 0 [, γ → θ M γ , for any θ ∈ R q , is a decreasing mapping on ]0, γ 0 [, and M γ → γ→0 I q (such a choice is motivated by the scheme described in Section 6.1). One can consider an alternative version of the SDE (4.1), i.e.,

dL(t) = - 1 2 ∇ L + ( M γ ,γ H) • D • M γ -1/2 (L(t))dt + M γ 1/2 dW (t), t > 0. (4.7)
Denote the drift coefficient of (4.7) by φ, we get that

φ(θ), θ = ψ(u), u ,
where u = M γ -1/2 θ. Therefore, it is easily seen that φ also satisfies (4.2) and (4.3) under assumptions (H.1)-(H.4). In turn, Proposition 4.2 applies to (4.7) the diffusion L is unique, non explosive and admits an unique invariant measure µ γ having density

θ → exp -L + ( M γ ,γ H) • D • M γ -1/2 (θ) /Z γ where Z γ = det(M γ ) R p exp -L + ( M γ ,γ H) • D (u) du. Since det(M γ ) → γ→0
1, applying the reasoning in the proof of Proposition 3.1, we also deduce that µ γ converges to µ in total variation as γ → 0.

By the change of variable U (t) = M γ -1/2 L(t), we get the following SDE

dU (t) = - 1 2 M γ -1 ∇ L + ( M γ ,γ H) • D (U (t))dt + dW (t), t > 0. (4.8)
In an analogous way to (4.5), the forward Euler discretization of (4.8) has a deterministic part which is a relaxed gradient descent in the metric M γ -1 . In turn, mean square convergence of the Euler discretizations of both (4.7) and (4.8) and of their first-order moments can be established exactly in the same way as in Theorem 4.1. We omit the details here for the sake of brevity.

Prox-regular penalties

We now present a large class of penalties, namely prox-regular functions, which satisfy the key assumptions (H.2) and (H.3).

Roughly speaking, a lsc function f is prox-regular at x ∈ dom(f ) if it has a "local quadratic support" at x for all (x, v) ∈ gph(∂f ) close enough to (x, v) ∈ gph(∂f ) with f (x) nearby f (x). This is formalized in the following definition. Definition 5.1 (Prox-regularity). Let f : R d → R ∪ {+∞}, given a point x ∈ dom(f ). f is prox-regular at x for v, with v ∈ ∂f (x) if f is locally lsc at x, there exist > 0 and r > 0 such that

f (x ) > f (x) + (x -x) v - 1 2r x -x 2 2 , when x -x 2 < and x -x 2 < with x = x and f (x) -f (x) 2 < while v -v 2 < with v ∈ ∂f (x).
When this holds for all v ∈ ∂f (x), f is said prox-regular at x. When f is prox-regular at every x ∈ dom(f ), f is said prox-regular.

Example 5.1. The class of prox-regular functions is large enough to include many of those used in statistics.

For instance, here examples where prox-regularity is fullfilled (see [START_REF] Rockafellar | Variational analysis[END_REF]Chapter 13, Section F] and [START_REF] Poliquin | Local differentiability of distance functions[END_REF]):

(i) Proper lsc convex functions.

(ii) Proper lsc lower-C 2 (or semi-convex) functions, i.e., f is such that f

+ 1 2r • 2 2 is convex, r > 0. (iii) Strongly amenable functions, i.e., f = g • R, R : R d → R q ∈ C 2 (R d
) and g : R q → R ∪ {+∞} proper lsc convex.

(iv) A closed set C ⊂ R d is prox-regular if, and only if, ι C is a prox-regular function. This is also equivalent to: for any x ∈ R d and for any γ > 0,

P C (x) = Argmin v∈R d 1 γ x -v 2 2 + ι C (v) = prox γι C (x)
is single valued and continuous, or equivalently, to

d 2 C = min v∈R d 1 γ • -v 2 2 + ι C (v) = γ ι C ∈ C 1,+ (R d ).
The following lemma summarizes a fundamental property of prox-regular functions.

Lemma 5.1 ([50, Theorem 3.2]). When f : R d → R ∪ {+∞} is locally lsc at x ∈ R d , the following are equivalent (i) f is prox-regular at x for v ∈ ∂f (x).
(ii) v is a proximal subgradient to f at x, i.e., there exist r > 0 and > 0 such that

f (x) ≥ f (x) + v, x -x - r 2 x -x 2 2 , ∀x such that x -x 2 < .
Moreover, there exist r > 0 and an f -attentive -localization (with > 0) of ∂f around (x, v) defined by

T f ,x,v (x) = v ∈ ∂f (x) : v -v 2 < if x -x 2 < and f (x) -f (x) 2 < , ∅ otherwise, such that T f ,x,v + rI d is monotone.
Let us consider a prox-regular function satisfying (H.1). Owing to the following lemma, such type of functions also fullfills (H. .

When prox M

γH is globally Lipschitz continuous, the optimal convergence rate in (4.6) is of order δ 1/2 in view of Theorem 4.1.

Forward-Backward type LMC algorithms

Let us now deal with our main goal: computing the EWA estimator in (1.4) by sampling from µ. Remind that

µ(θ) ∝ exp - F (Xθ, y) + J λ (θ) β ,
where F : R n × R n → R is a general loss and J λ : R p → R is the penalty. Assume that F (X•, y) ∈ C 1,+ (R p ) and the penalty takes the form J λ = W λ • D . Let us impose the following assumptions on W λ .

(H.1') W λ : R q → R is lsc and bounded from below.

(H.2') prox γW λ is single valued.

(H.3') prox γW λ is locally Lipschitz continuous.

To lighten notation, we will write F β def = F (X•, y)/β. This section aims to describe our Forward-Backward type Langevin Monte-Carlo (LMC) algorithms to implement (1.4). These algorithms are based on wise specializations of the results reported in Section 4.

Forward-backward LMC (FBLMC)

In (3.1), we set D = I p (hence J λ = W λ ), L ≡ 0, and H = F β + J λ /β, where F is a quadratic loss, i.e., F β (θ) = y -Xθ

prox γW λ /β • (I p -γ∇F β )(θ), θ M γ ≤ C FB (1 + θ 2 2 ), ∀θ ∈ R p .
By Lemma 2.1, a sufficient condition for (H.4'-FB) to hold is that the proximal mapping of W λ is Lipschitz continuous. From Lemmas 3.2 and 6.1, we get

∇ M γ ,γ H = γ -1 M γ I p -prox M γ γH = γ -1 M γ I p -prox γJ λ /β (I p -γ∇F β ) .
With this expression at hand, the forward Euler discretization of the SDE (4.1), specialized to the current case, reads

L k+1 = L k - δ 2γ M γ L k -prox γJ λ /β (L k -γ∇F β (L k )) + √ δZ k , t > 0, L 0 = l 0 , (6.2)
Similarly, the forward Euler discretization of the SDE (4.8) is given by

U k+1 = (1 -δ 2γ )U k + δ 2γ prox γJ λ /β (U k -γ∇F β (U k )) + √ δZ k , t > 0, U 0 = l 0 . (6.
3)

The familiar reader may have recognized that the deterministic part of (6.3) is nothing but the relaxed form of the so-called Forward-Backward proximal splitting algorithm [START_REF] Bauschke | Convex analysis and monotone operator theory in Hilbert spaces[END_REF]. This terminology reflects that there is a forward Euler discretization on F β and a Euler backward discretization on J λ .

Semi-Forward-Backward LMC (Semi-FBLMC)

The main limitation of (6.2) is that the proximal mapping of J λ must be easy to compute. This may not be true even if the proximal mapping of W λ is accessible as, for for example, when D does not have orthogonal rows [START_REF] Bauschke | Convex analysis and monotone operator theory in Hilbert spaces[END_REF]. Our goal is to circumvent this difficulty. Toward this goal, in (3.1), consider now L = F β , H = W λ /β and M = I q . Owing to (H.1')-(H.3'), one can check that (H.1)-(H.3) are fulfilled. Assumption (H.4) is verified by imposing the following on W λ .

(H.4'-SFB) There exists

C SFB > 0 such that prox γW λ /β (u), u ≤ C SFB (1 + u 2 2 ), ∀u ∈ R q .
From Lemma 3.2, we obtain

∇ ( γ H) • D (θ) = γ -1 D(D θ -prox γW λ /β (D θ)).
Thus, the forward Euler discretization of SDE (4.1) now reads

L k+1 = L k - δ 2 ∇F β (L k ) - δ 2γ D D L k -prox γW λ /β (D L k ) + √ δZ k , t > 0, L 0 = l 0 . (6.4)
In the case where D = I p , F β and W λ are convex, we recover the scheme studied in [START_REF] Durmus | Sampling from convex non continuously differentiable functions, when Moreau meets Langevin[END_REF].

Applications to penalties in statistics

In this section, we exemplify our LMC sampling algorithms for some popular penalties in the statistical and machine learning literature. Our goal is by no means to be exhaustive, but rather to be illustrative and show the versatility of our framework. For each penalty, we aim at checking that assumptions (H.1')-(H.3'), (H.4'-FB) and (H.4'-SFB) hold, and to compute prox γW λ /β . In turn, this allows to apply our algorithms (6.3) and (6.4) to compute EWA with such penalties.

Analysis group-separable penalties

We first focus on a class of penalties where J λ is analysis group-separable, i.e.,

J λ (θ) = W λ (D θ) where W λ (u) = L l=1 w λ u G l 2 , (7.1)
for w λ : R + → R, and some uniform partition (G l ) l∈{1,...,L} of {1, . . . , q}, i.e., ∪ L l=1 G l = {1, . . . , q} and G l ∩ G l , ∀l = l . Remark 7.1. It is worth mentioning that separability of W λ does not entail that of J λ . In fact, overlapping groups can be easily taken intro account as any overlapping-group penalty can be written as the composition of W λ with a linear operator, say B, such that B B is diagonal, and B acts as a group extractor, see [START_REF] Chen | An efficient proximal-gradient method for general structured sparse learning[END_REF][START_REF] Peyré | Group sparsity with overlapping partition functions[END_REF].

A first consequence of separability is that prox γW λ /β is also separable under the following mild assumptions on w λ .

(W.1) w λ is bounded from below on ]0, +∞[. (W.2) w λ is non-decreasing functions on ]0, +∞[. Lemma 7.1. Assume that Assumptions (W.1) and (W.2) hold, and w λ is continuous on ]0, +∞[. For any u ∈ R q and γ > 0, we have

prox γW λ /β (u) =      prox γw λ /β ( u G 1 2 ) u G 1 u G 1 2 . . . prox γw λ /β u G L 2 u G L u G L 2      .

SCAD penalty

The SCAD penalty, proposed in [START_REF] Fan | Variable selection via nonconcave penalized likelihood and its oracle properties[END_REF] is parameterized by λ = (λ, a) ∈]0, +∞[×]2, +∞[ as

w λ (t) =      λt if 0 ≤ t ≤ λ, -t 2 -2aλt+λ 2 2(a-1) if λ < t ≤ aλ, (a+1)λ 2 2 if t > aλ, (7.5)
The following lemma establishes the validity of w λ and computes prox γw λ /β . Lemma 7.4. Let w λ defined in (7.5), and κ = γ/β. For any γ < (a -1)β,

(i) w λ satisfies (W.1) -(W.3), (
ii) The proximal mapping of the SCAD penalty is given by the shrinkage

prox γw λ /β (t) =      (t -κλ) + if 0 ≤ t ≤ (κ + 1)λ, (a-1)t-kaλ a-1-κ if (κ + 1)λ < t ≤ aλ, t if t > aλ. (7.6)
Since a > 2, one can set κ = 1. In this case, (7.6) specializes to [START_REF] Fan | Variable selection via nonconcave penalized likelihood and its oracle properties[END_REF]Equation (2.8)].

∞ penalty The ∞ norm penalty is convex and continuous but is not separable, unlike the previous ones. It is a suitable prior to promote flat vectors, and has found applications in several fields [START_REF] Jégou | Anti-sparse coding for approximate nearest neighbor search[END_REF][START_REF] Lyubarskii | Uncertainty principles and vector quantization[END_REF][START_REF] Studer | Signal representations with minimum ∞ -norm[END_REF]. It entails the following penalty W λ :

J λ (θ) = W λ (D θ) where W λ (u) = λ max l∈{1,...,L} [u] G l 2 , (7.7)
where λ = λ > 0. Since W λ is not separable, Lemma 7.1 is not applicable. Nevertheless, the proximal mapping of W λ can still be obtained easily from the projector on the 1 unit ball, i.e.,

prox γW λ /β (u) = u -P x : l x G l 2 ≤ β λγ (u). (7.8)
This projector can be obtained from [31, Proposition 2] (see also references therein). One can see that (H.1')-(H.3') hold. We report the verification of (H.4'-FB) and (H.4'-SFB) in the proof of the following lemma.

Lemma 7.5. Let W λ in (7.7). Then (H.4'-FB) and (H.4'-SFB) hold.

Numerical experiments

In this section, some numerical experiments are conducted to illustrate and validate our LMC algorithms.

Image processing experiments

Let θ 0 is a 2-D image which is a matrix in R 128×128 . Let us denote vec : R √ p× √ p → R p the vectorization operator, i.e. the operator which stacks the columns of its arguments. We then consider the following linear regression problem y = X vec(θ 0 ) + ξ.

Here p = 1282 and ξ ∼ N (0, σ 2 I n ). The noise level σ is chosen according to the simulated θ 0 through the signal-to-noise ratio SNR, i.e. σ = n -1/2 Xθ 0 2 /10 SNR/10 . In our experiments, we take SNR = 5.

The goal is estimating θ 0 by computing the EWA estimators via the penalties proposed in Section 7. Three types of problems are considered: compressed sensing, inpainting and deconvolution whose regression function described in what follows.

• Compressed sensing: in this case X is drawn from a random ensemble. In our experiments, X is drawn uniformly at random from the Rademacher ensemble, i.e., its entries are iid Rademacher random variables. We also set n = 9p/16.

• Inpainting In this case, X acts as a masking operator. Let M ⊂ {1, . . . , p} be the set indexing masked pixels. Thus X vec(θ 0 ) = vec(θ) j∈{1,...,p}\M .

In our numerical experiments, we mask out 20% of the pixels, and thus n = p-20%p where p the stands of integer part of p. The masked positions are chosen randomly from the uniform distribution.

• Deconvolution In this case X is the convolution operator with a Gaussian kernel with periodic boundary conditions, such that X is diagonalized in the discrete Fourier basis. In this experiment, the standard deviation of the kernel is set to 1.

Assuming that the targeted image is piecewise smooth, a popular prior is the so-called isotropic total variation [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF]. To cas this into our framework, define D c : R

√ p× √ p → R √ p× √ p and D r : R √ p× √ p → R √ p×
√ p the finite difference operators along, respectively, the columns and rows of an image, with Neumann boundary conditions. We define D TV as

D TV : θ ∈ R √ p× √ p → vec (vec(D r (θ)), vec(D c (θ))) ∈ R 2p .
With this notation, our prior penalty J λ reads

J λ (θ) = p l=1 w λ vec(D r (θ)) 2 l + vec(D c (θ)) 2 l = W λ (D TV θ), (8.2)
which clearly has the form (7.1) with p blocks of equal size 2.

To estimate θ 0 from (8.1), we employ the EWA estimator (1.4) with F (Xθ, y) = y -X vec(θ) 2 2 and J λ in (8.2). For each problem instance (compressed sensing, inpainting or deconvolution), we tested w λ as the 1 , SCAD and FIRM penalties. The corresponding EWA estimators are denoted respectively EWA-1 , EWA-SCAD and EWA-FIRM. Because of the presence of the analysis operator D TV , which is not unitary, we applied Semi-FBLMC scheme (6.4) to compute EWA with β = 1/(pn), γ = β, and δ = 5β/10 3 , 5β/10 2 , 5β/10 6 respectively associated to inpainting, deconvolution and compressed sensing problems. The results are depicted in Figure 1.

Signal processing experiments

Here we consider reconstructing a piecewise flat 1D signal from compressed sensing measurements using EWA. For this, we create a p = 128 sample signal whose coordinates are valued in {-1, 1} and compute the observations (8.1) where X is drawn from the Rademacher ensemble with n > p 2 . We set F (Xθ, y) = y -Xθ 2 2 , J λ (θ) = θ ∞ , i.e. D = I p and the size of groups is 1. We can then use the FBLMC scheme (6.3), where we choose β = 1/(pn), γ = β, and δ = 5/10 2 . The results are shown in Figure 2.

f (x), x ≤ f (x) 2 2 /2 + x 2 2 /2 ≤ f (x) -f (x * ) 2 2 + f (x * ) 2 2 + x 2 2 /2 ≤ K x -x * 2 2 + f (x * ) 2 2 + x 2 2 /2 ≤ 2 K + 1/2 x 2 2 + 2 K x * 2 2 + f (x * ) 2 2 ≤ K(1 + x 2 2 ), with K ≥ max 2 K + 1/2, 2 K x * 2 2 + f (x * ) 2 2 .
Recalling that f is bounded on bounded sets concludes the proof.

Proof of Lemma 3.1 (i) In view of (H.1), H is prox-bounded by [START_REF] Rockafellar | Variational analysis[END_REF]Exercise 1.24] for any γ ∈]0, γ 0 [, and then for any x, (ii) The continuity and finiteness properties follow from [55, Theorem 1.17(c)] (see also [START_REF] Rockafellar | Variational analysis[END_REF]Theorem 1.25]), where we use (H.1). For the second claim, we have ∀x

∈ R q -∞ < inf H ≤ M γ ,γ H(x) ≤ H(x).
Moreover, let p ∈ prox

M γ γH (x). Then, ∀δ > γ, M δ ,δ H(x) = inf w∈R q 1 2δ w -x 2 M δ + H(w) ≤ 1 2δ p -x 2 M δ + H(p) ≤ 1 2γ p -x 2 M γ + H(p) = M γ ,γ H(x).
This together with continuity concludes the proof of Assertion (ii).

Proof of Lemma 3.2 By virtue of Lemma 3.1-(i) and (H.2), prox

M γ
γH is clearly non-empty and single valued. The continuity property follows from [55, Theorem 1.17(b)] (see also [START_REF] Rockafellar | Variational analysis[END_REF]Theorem 1.25]) and single-valuedness. By Lemma 3.1-(ii), M γ ,γ H(θ) is finite. Since (H.1) holds, H is prox-bounded with threshold ∞ by [START_REF] Rockafellar | Variational analysis[END_REF]Exercise 1.24]. Invoking [START_REF] Rockafellar | Variational analysis[END_REF]Example 10.32], we get that -M γ ,γ H is locally Lipschitz continuous, subdifferentially regular and

∂ -M γ ,γ H (θ) = γ -1 M γ prox M γ γH (θ) -θ , ∀θ ∈ R p .
Combining this with [START_REF] Rockafellar | Variational analysis[END_REF]Theorem 9.18] applied to -M γ ,γ H, we obtain that M γ ,γ H is differentiable and its gradient is precisely as given.

Proof of Proposition 3.1 Recall that we denote with the same symbol the measure and its density with respect to the Lebesgue measure. Thus

µ γ -µ TV = R M µ γ (θ) -µ(θ) dθ, where µ γ (θ) = exp -L(θ) + ( M ,γ H) • D (θ) /Z γ and µ(θ) = exp -L(θ) + H • D (θ) /Z, and Z = R M exp -(L(θ ) + H • D (θ )) dθ .
In view of Lemma 3.1(ii), applying the monotone convergence theorem, we conclude that Z γ → Z when γ → 0. This together with Lemma 3.1(ii) again yield that µ γ converges to µ pointwise. We conclude using Scheffé(-Riesz) theorem [START_REF] Kusolitsch | Why the theorem of scheffé should be rather called a theorem of riesz[END_REF][START_REF] Scheffe | A useful convergence theorem for probability distributions[END_REF].

Proof of Proposition 4.1 In view of Lemma 3.2, the drift term reads

ψ(θ) = - 1 2 ∇(L + ( M ,γ H) • D )(θ) = - 1 2 ∇L(θ) - 1 2γ DM D θ + 1 2γ DM prox M γH (D θ).
Since L ∈ C 1,+ (R p ) and (H.4) holds, there exist K 1 > 0 and K 2 > 0 such that

ψ(θ), θ = - 1 2 ∇L(θ), θ - 1 2γ D θ 2 M + 1 2 prox M γH (D θ), D θ M ≤ K 1 (1 + θ 2 2 ) + |||D||| 2 |||M |||/(2γ) θ 2 2 + K 2 (1 + θ 2 2 ) ≤ K(1 + θ 2 2 ),
where

K ≥ K 1 + K 2 + |||D||| 2 |||M |||/(2γ
> 0 such that H(x ) > H(x) + v, x -x - 1 2λ 0 x -x 2 2 > H(x) + v, x -x - 1 2λ 0 σ min (M ) x -x 2 M , (9.1)
for any x = x and (

x, v) ∈ gph T H ,x,v . Let γ 0 = λ 0 σ min (M ), γ ∈]0, γ 0 [ and u = x + γM -1 v, (9.1) becomes H(x ) + 1 2γ x -u 2 M > H(x) + 1 2γ x -u 2 M .
Therefore, prox M γH (u) = x where (x, v) ∈ gph T H ,x,v . That yields prox M γH (x + γM -1 v) = x. Since H is lsc, proper and prox-bounded, from [55, Theorem 1.17(c)] (see also [START_REF] Rockafellar | Variational analysis[END_REF]Theorem 1.25]), we have We conclude the proof of Lemma 6.1.

x ∈ prox M γH (u), u → x + γM -1 v =⇒ x → prox M γH (x + γM -1 v) = x, H(x) = M ,γ H(u) -1 2γ x -u 2 2 → H(x).
Proof of Lemma 7.1 This is a probably known result, for which we provide a simple proof. Since W λ is separable and w λ is continuous and lower-bounded, we have (i) Observe that w λ is continuously differentiable on ]0, +∞[ with w λ (t) = κλ I(t ≤ λ) + (aλ -t) + (a -1)λ I(t > λ) ≥ 0, w λ is then non decreasing and bounded from below by w λ (0) = 0 on ]0, +∞[. Thus, w λ satisfies (W.1) and (W.2). Let us check (W.3). Let u(t) = t + κw λ (t), we obtain that

• u(0) = κλ, • if 0 < t ≤ λ, u(t) = t + κλ > κλ,

• if λ < t ≤ aλ, since a -1 > κ > 0, u(t) = t + κ(aλ-t) a-1 = κλ + a-1-κ a-1 t + κλ a-1 > κλ, • if t > aλ, since a -1 > κ, u(t) = t > aλ > κλ.

Thus, t = 0 is the unique mimimum in [0, +∞[ of t + p λ (t). In other words, w λ satisfies (W.3).

(ii) For the sake of simplified notation, we denote p = prox γw λ /β (t). Owing to Lemma 7.2, we obtain that p = 0 if t ≤ κλ, t -κλ I(p ≤ λ) + (aλ-p) + (a-1)λ I(p > λ) otherwise. (9.7)

From (9.7), we get the following assertions when t > κλ,

• if p ≤ λ, p = t -κλ, and t = p + κλ ≤ (κ + 1)λ,

• if λ < p ≤ aλ, p = t -κ(aλ -p)/(a -1) implies that p = (a-1)t-κaλ a-1-κ . Since λ < p ≤ aλ, κ < a -1 and a > 2, we also get that

(1 + κ)λ < t = a -1 -κ a -1 p + κaλ a -1 ≤ aλ,
• if p > aλ, p = t, and t > aλ.

That concludes the proof of (ii), Lemma 7.4.

Proof of Lemma 7.5 Set u = D θ, α = γλ/β and p u = P Thus, Assumption (H.4'-FB) holds and we conclude the proof of Lemma 7.5.

Definition 2 . 2 (

 22 and f is subdifferentially regular at any point x where ∂f (x) = ∅. Proximal mapping and Moreau envelope). Let M ∈ R d×d symmetric positive definite. For a proper lsc function f and γ > 0, the proximal mapping and Moreau envelope in the metric M are defined respectively by

Proposition 4 . 1 .

 41 Assume that (H.1)-(H.4) hold. Then,

Figure 1 :

 1 Figure 1: (a): Original image. (b,c) Observed masked and blurry images. (d, e, f): EWA-1 estimated images from masked image, compressed sensing measurements, and blurry image. (g, h, i): EWA-FIRM estimated images from masked image, compressed sensing measurements, and blurry image. (j, k, l): EWA-SCAD estimated images from masked image, compressed sensing measurements, and blurry image. 16

Figure 2 :

 2 Figure 2: Compressed sensing with EWA using the ∞ penalty. * is the original signal and • is the the estimated one.

1

 1 Let x * ∈ C, a bounded subset of R d . Using Young and Jensen inequalities as well as K-Lipschitz continuity of f , we obtain

1 2γ x -• 2 M

 12 γ + H is proper lsc and level-bounded uniformly in (x, γ) ∈ R q ×]0, γ 0 [. This entails that the set of minimizers of this function, i.e. prox M γ γH (x), is a non-empty compact set. For the last claim, suppose that x ∈ Argmin(H) = ∅ and bounded but x / ∈ prox M γ γH (x). Thus, for any p ∈ prox M γ γH (x), we have p = x andH(p) < 1 2γ pxleading to a contradiction with x is a minimizer of H.

(9. 2 ) 2 M 2 .

 222 and the latter operator is well-defined as a single-valued operator since S is maximal monotone; see[5, Proposition 3.22 (ii)(d)]. Let p = prox M γH (θ) and p = prox M γH (θ ). It then follows that M θ -γδM p ∈ γS(p) and M θ -γδM p ∈ γS(p ), and monotonicity of S yields pp, M (θ -θ) ≥ γδ pp ≥ γδσ min (M ) pp 2 Using Cauchy-Schwartz inequality, we then obtain pp 2 ≤ κ θθ 2 ,where κ -1 = γδσ min (M )

2 2 -

 2 X(θw), Xθy + γJ λ /β (θ -γ∇F (θ)).

  v 2 ),and thus, ∀l ∈ {1, . . . , L},prox γW λ /β (u)

  and L(t) are continuous-time extensions to the discrete-time chain {L k } k .Mean square convergence of the pathwise approximation (4.5) and of its first-order moment can be established as follows. Assume that (4.2) and (4.3) hold and E

	Theorem 4.1.

  2) and (H.3). Let M ∈ R p×p symmetric positive definite and γ small enough, assume that H : R p → R is prox-regular and satisfies (H.1). Then prox M γH is single-valued and locally Lipschitz continuous. Lower-C 2 (or semi-convex) functions, see Example 5.1-(ii), satisfy the global counterpart of Lemma 5.1-(ii). For a lower-C 2 penalty H satisfying (H.1), the following lemma shows that prox M γH is globally Lipschitz continuous with a proper choice of γ which in turn implies directly (H.4) according to Lemma 2.1. Assume that H is lower-C 2 (with constant r) satisfying (H.1) and γ ∈]0, rσ min (M )[, prox M

	Lemma 5.2. Lemma 5.3.

γH is single-valued and Lipschitz continuous with constant |||M ||| σ min (M ) 1 -γ rσ min (M ) -1

2 2

 2 /β. Observe that H satisfies (H.1) by assumption (H.1'). To check (H.2)-(H.4), we need to design a metric in which prox M γH is expressed as a function of prox γJ λ /β . This idea is formalized in the following lemma. = prox γJ λ /β • (I p -γ∇F β ).

	Lemma 6.1. Assume that 0 < γ < β/(2|||X||| 2 ) and (H.1') holds. Define M γ	def = I p -(2γ/β)X X, which
	is symmetric positive definite. Then
	prox	M γ	(6.1)
	In view of Lemma 6.1, (H.2') and (H.3'), it is immediate to check that (H.2) and (H.3) are satisfied.
	It remains now to verify (H.4) which is fulfilled by imposing the following assumption on W λ (or J λ ).

γH

(H.4'-FB) There exists C FB > 0 such that

  ). Moreover, under (H.3), ( M ,γ H) • D is locally Lipschitz continuous by virtue of Lemma 3.2, which applies thanks to assumptions (H.1)-(H.2). Clearly( M ,γ H) • D ∈ C 1,+ (R p ). Since C 1,+ (R p) is closed under addition, we conclude the proof. The proof of Lemma 5.3 is based on the one of [55,Proposition 13.37] and generalizes to the proximal mapping in metric M for any M ∈ R p×p symmetric positive definite.Without loss of generality, we prove the claim on a neighbourhood of x where H is lsc. Let x ∈ R p , v ∈ ∂H(x), since H is prox-regular at x for v and H is prox-bounded, owing to [7, Lemma 4.1], there exist > 0 and λ 0

	Proof of Proposition 4.2 Claim (i) follows by combining Proposition 4.1 and [69, Theorem 3.6, Chapter
	II]. Claim (ii) is a consequence of Proposition 4.1 and [54, Theorem 2.1].

Proof of Theorem 4.1 Owing to Proposition 4.1 and [69, Theorem 4.1, Chapter II], we get that the r-th moments of L(t) are bounded for any r ≥ 2 and t ≥ 0. A similar reasoning also entails that the r-th moments of the continuous-time extension L δ are also bounded. Moreover, according to Proposition 4.1, the drift ψ is locally Lipschitz continuous. The claim then follows from [38, Theorem 2.2] and Jensen's inequality. In the globally Lipschitz continuous case, we get the claimed rate by putting together Lemma 2.1, Jensen's inequality and [69, Theorem 7.3, Chapter II] or [40, Theorem 10.2.2 and Remark 10.2.3]. Proof of Lemma 5.2

  Owing to(7.8) and Young's inequality, we obtain thatu, prox γW λ /β (u) = u, up u ≤ u 2 2 + u 2 p u 2 ≤ ≤ |||D||| 2 θ 2 2 , (H.4'-SFB) holds. Set v = 2γX y/β, t θ = θ -γ∇F β (θ) = M γ θ + v and p t θ = P x : α l x G l 2 ≤1(t θ ). By Young's inequality, we obtain thatprox γW λ /β (t θ ), θ M γ = M γ prox γW λ /β (t θ ), θ ≤ 1 2 |||M γ ||| 2 prox γW λ /β (t θ )

							3 2	u 2 2 +	1 2	p u	2 2 ≤	3 2	u 2 2 +	1 2α 2 .
	According to the fact that u 2 2 = D θ	2 2 2 2	+	1 2	θ 2 2 .
	Moreover, owing to (7.8), we get that
	prox γW λ /β (t θ )	2 2	= t θ -p t θ	2 2 ≤ 2 t θ	2 2 + 2 p t θ	2 2 ≤ 4|||M γ ||| 2 θ 2 2 + 4 v 2 2 +	2 α 2 .

x : α l x G l 2 ≤1 (u).

The author however applied it to problems where -log µ = L + H. But the gradient of the Moreau envelope of a sum, which amounts to computing the proximity operator of -log µ does not have an easily implementable expression even if those of L and H do.

The overdetermined regime is known to yield good performance for the ∞ penalty[START_REF] Vaiter | Model selection with low complexity priors[END_REF].

M γ + H(p) ≤ H(x),

Acknowledgement. This work was supported by Conseil Régional de Basse-Normandie and partly by Institut Universitaire de France.

Our aim is now to design a family of penalties that will allow to establish (H.1')-(H.3'), (H.4'-FB) and (H.4'-SFB), while involving a form of shrinkage which is ubiquitous in low-complexity regularization. Inspired by the work of [START_REF] Antoniadis | Regularization of Wavelet Approximations[END_REF], we make the following assumptions on w λ .

(W.3) w λ is continuously differentiable on ]0, +∞[ and the problem min t∈[0,+∞[ {t + γ β w λ (t)} has a unique solution at 0 for a given γ.

Under these assumptions, prox γw λ /β has a indeed convenient shrinkage-type form. Lemma 7.2 ([2, Theorem 1]). Assume that (W.2) and (W.3) hold for some γ > 0. Then, prox γw λ /β are the single-valued continuous mappings, and satisfy, for t ∈ [0, +∞[,

Let us turn to check our assumptions. (H.1')-(H.3') are fulfilled thanks to (W.1)-(W.3). It remains to check (H.4'-FB) and (H.4'-SFB). This is the subject of the following lemma. Lemma 7.3. Assume that (W.2) and (W.3) hold for some γ > 0, then (H.4'-FB) and (H.4'-SFB) also hold.

We now discuss some popular penalties w λ that satisfy (W.1)-(W.3) for some γ > 0.

Examples

1 penalty Take w λ : t ∈ R + → λt. This entails the analysis group Lasso penalty

Clearly, w λ is a continuous positive convex function which verifies (W.1)-(W.3) for any γ > 0, and its proximal mapping corresponds to soft-thresholding, i.e., prox γw λ /β (t) = (t -γλ/β) + , ∀t ≥ 0.

FIRM penalty

The FIRM penalty is given by [START_REF] Gao | Waveshrink with firm shrinkage[END_REF] w

which entails the corresponding analysis group FIRM penalty

≥ 0, w λ is non-decreasing and bounded from below by w λ (0) = 0 on ]0, +∞[. Thus, w λ satisfies (W.1) and (W.2). Assumption (W.3) also holds for any γ < βµ/λ. The operator prox γw λ /β can be constructed from [START_REF] Woodworth | Compressed sensing recovery via nonconvex shrinkage penalties[END_REF]Definition II.3]. Its formula is defined as

where α = γλ/β. The formula (7.4) can also be found using Lemma 7.2. Observe that the FIRM shrinkage (7.4) interpolates between hard-(see [START_REF] Woodworth | Compressed sensing recovery via nonconvex shrinkage penalties[END_REF]Definition II.2]) and soft-thresholding. In particular, (7.4) coincides with soft-thresholding when µ → ∞.

For any x ∈ prox M γH (u), by Fermat rules we get

For any γ ∈]0, γ 0 [, owing to (9.2) and (9.3), there exists

So that

where δ = 1/γ -1/γ 0 , S = T H ,x,v + 1/γ 0 M . From (9.1), S is maximal monotone, the latter operator is well defined as a single valued operator (see [

Using Cauchy-Schwarz's inequality, we obtain

where

Let us note that when γ decrease, Inequality (9.1) can be hold for a larger that enlarges N γ,x,v and x + γM -1 v concentrate to x for any v. Thus, when γ is small enough, there exists a neighbourhood x that includes in N γ,x,v for any v ∈ ∂H(x). That concludes the proof of Lemma 5.3.

Proof of Lemma 5.3 From [55, Example 12.28(b)], ∂H is hypomonotone of modulus 1 r . In turn

Thus

is strongly convex uniformly in θ and γ complying with γ < γ 0 . It then follows that prox M γH is single-valued. We have

By Fermat's rule, we then get

For u G l = 0, by isotropy of problem (9.4), we can write

The inner minimization problem amounts to solving for the orthogonal projector on the 2 sphere in R G of radius t, which is obviously

since u G l = 0. Inserting this into (9.5) and rearranging the terms, (9.4) becomes

where we used even-symmetry of w λ .

Proof of Lemma 7.3 Before proceeding, let us discuss about the term prox γw λ /β . In view of (W.2), w λ /β is positive on ]0, +∞[. According to Lemma 7.2 we get that, for any t ≥ 0, prox

Set u = D θ, from Lemma 7.1 and (9.6), we get that

According to the fact that u 2 2 = D θ

Thus, Assumption (H.4'-FB) holds and we conclude the proof of Lemma 7.3.


Martin Rosalie
email: martin.rosalie@uni.lujan.dentler@uni.lusomasundar.kannan@uni.lumiguel.olivaresmendez@uni.lu

Jan E Dentler

Grégoire Danoy
email: gregoire.danoy@uni.lupascal.bouvry@uni.lu

Pascal Bouvry

Somasundar Kannan

Miguel A Olivares-Mendez

Holger Voos
email: holger.voos@uni.lu

Area exploration with a swarm of UAVs combining deterministic chaotic ant colony mobility with position MPC

published or not. The documents may come

I. INTRODUCTION

The extensive operational space and sensing capabilities of UAVs facilitate the exploration of areas. By using multiple autonomously operating UAVs, the response time in search and rescue scenarios within large areas can be reduced significantly. However, the coverage of the surveyed area is dependent on the coordination of the UAVs. A comprehensive overview on mechanisms to coordinate and control such swarms are given in [START_REF] Gazi | Swarm Coordination and Control Problems[END_REF] and [START_REF] Bertrand | MPC Strategies for Cooperative Guidance of Autonomous Vehicles[END_REF]. A detailed recent state of art of the multi-UAV area exploration problem from a perspective of optimization and artificial intelligence is summarized in [START_REF] Nigam | The multiple unmanned air vehicle persistent surveillance problem: A review[END_REF]. A typical approach which is covered in all of these contributions is a cooperative area exploration control with potential fields over a receding horizon. The nature inspired interpretation of this method has been first introduced in [START_REF] Gaudiano | Swarm intelligence: a new C2 paradigm with an application to control of swarms of UAVs[END_REF] as Ant Colony Optimization with pheromones to coordinate a swarm of UAVs. This deterministic path planning is in contrast to Dynamic Data Driven Application Systems (DDDAS, see [START_REF] Mccune | Swarm control of UAVs for cooperative hunting with DDDAS[END_REF] for an example applied on swarm of UAVs), where the UAVs are coordinated according to the data they collected.

Kuiper and Nadjm-Tehrani [START_REF] Kuiper | Mobility models for UAV group reconnaissance applications[END_REF] define an area exploration scenario with ten UAVs with a pheromone based approach for a path planning. Their study compares a random mobility model to a pheromone-based mobility model using few metrics for coverage. Some recent studies propose to replace the random part of such metaheuristics to increase their performance. For instance Gandomi et al. use this principle for Particle Swarm Optimization method [START_REF] Gandomi | Chaosenhanced accelerated particle swarm optimization[END_REF] and Firefly Algorithm [START_REF] Gandomi | Firefly algorithm with chaos[END_REF]. Thus, in our previous work [START_REF] Rosalie | From random process to chaotic behavior in swarms of UAVs[END_REF] we use chaotic dynamic to create the Chaotic Ant Colony Optimization to Coverage (CACOC) mobility model to increase the area coverage performance of an Ant Colony Algorithm. However, the high-level simulations used to evaluate CACOC do not take into account the dynamics of a real quadrotor, which are limiting the unpredictability of a real UAV trajectory. To analyze the influence of the UAV dynamics, the major contribution of this paper is a stochastic evaluation of the performance of CACOC in terms of coverage in a realistic scenario.

The structure of this use case exploration scenario is illustrated in Fig. 1. The latter consists of ten UAV models which are approximating a real quadrotor motion as introduced in some of the authors previous work [START_REF] Dentler | A realtime model predictive position control with collision avoidance for commercial low-cost quadrotors[END_REF]. The positions of these UAVs are controlled with model predictive control (MPC) as presented in the same paper. The system dynamics are equal to the MPC inherent prediction model which has been applied in [START_REF] Dentler | A realtime model predictive position control with collision avoidance for commercial low-cost quadrotors[END_REF] on real quadrotors and therefore justifies its representative usage. To finally assess the influence of the closed-loop UAV dynamics, the waypoints which are generated with the CACOC mobility model (indicated as mobility model) are given as target position for the MPC controlled UAV (indicated as control). The remainder of this paper is organized as follows. We first provide the description of the CACOC UAV mobility model and the control model of the UAVs in section II. Then the experimentation environment, the metrics used to compare the model and statistical analysis are given in section III. In section IV we finally provide conclusions and future work perspectives that this approach provides.

II. UAV MOBILITY MODEL

This section first provides a brief description of the CA-COC mobility model used to generate the waypoints to be reached by the UAVs. The second part introduces the model predictive position control (MPC) of a real quadrotor used in the experiments.

A. CACOC mobility model (Waypoints)

In order to generate the different waypoints to be followed by the UAVs, the CACOC mobility model introduced in one of our previous works [START_REF] Rosalie | From random process to chaotic behavior in swarms of UAVs[END_REF] has been used. In a high level simulation process, CACOC is used to describe the mobility of a swarm of ten UAVs. Its performance was compared to existing models using random, chaos and classical ant colony optimization algorithms [START_REF] Rosalie | From random process to chaotic behavior in swarms of UAVs[END_REF]. We empirically demonstrated that CACOC permits an increase of performance for several aspects of coverage (overall coverage, recent coverage and fairness of the coverage).

This section provides a brief description of this mobility model, for more details the interested reader should refer to [START_REF] Rosalie | From random process to chaotic behavior in swarms of UAVs[END_REF]. CACOC is based on the idea that the chaotic behavior of a dynamical system can enhance the performance of optimization algorithms. More precisely, CACOC replaces the random parts of the ant colony optimization algorithm [START_REF] Dorigo | Optimization, learning and natural algorithms[END_REF] by a chaotic behavior in such a way that the exploration capabilities of the UAV swarm is improved. The CA-COC performance on the coverage problem was empirically demonstrated using metrics as described below in Sec. III.

For the chaotic dynamics, we used the Rössler system [START_REF] Rössler | An equation for continuous chaos[END_REF] because it is a well-known dynamical system that can exhibit various chaotic mechanisms [START_REF] Rosalie | Templates and subtemplates of rössler attractors from a bifurcation diagram[END_REF]. CACOC works as follows: UAVs have a constant speed and choose at each step a direction: A for ahead, R for 45 • right and L for 45 • left. If there is no virtual pheromone to guide the UAV (pheromones are deposited by each UAV to indicate areas they already visited), then the next choice is given by the first return map (Fig. 2). This map underline the dynamical signature of the Rössler system giving ρ n+1 versus ρ n . Thus, the next action depends on the previous one:

• if ρ n < 1/3 then direction is right; • if 1/3 < ρ n < 2/3 then direction is left;
• else the direction is ahead. In that case, the good exploration performance of the UAVs is due to the periodic orbits of the system that lead to patterns (Fig. 2). If there are pheromones, they have repulsive properties and the next choice depends on the total amount of pheromones sensed around the UAV. Such a repulsive effect prevents the UAVs to reach recently visited areas since pheromones evaporate in time. We also used the first return map values to choose the next direction with the pheromones perception instead of a random number. Consider that p L , p A and p R are inversely proportional to the total amount of pheromones sensed respectively to the left, ahead and right of the UAV and that p R + p L + p A = 1. Thus, with ρ n taken from the first return map (Fig. 2), the next direction is chosen according to these rules:

• if ρ n < p R next direction is right; • if p R < ρ n < p R + p L next direction is left;
• else the direction is ahead. We demonstrated that the performance of the mobility model is increased with this chaotic behavior instead of a random one (see [START_REF] Rosalie | From random process to chaotic behavior in swarms of UAVs[END_REF] for details). CACOC relies on a simplified quadrotor model, this article therefore proposes to use a realistic one to evaluate CACOC performance in a physical environment. In this work, the CACOC mobility model provides waypoints that UAVs will try to reach. For the sake of clarity, the mobility model is referred to as Waypoints in the remainder of this article. As representation of a real quadrotor system, the motion model of an AR.Drone 2.0 is used in this work. Its system dynamics according to Fig. 3 and model predictive control (MPC) have been presented in some of the authors previous work [START_REF] Dentler | A realtime model predictive position control with collision avoidance for commercial low-cost quadrotors[END_REF]. Considering multiple quadrotors, the dynamics of quadrotor i can be described, using state x i and control vector u i :

B. UAV control (Controlled)

x i = x i , y i , z i , Ψ i , v f,i , v s,i (1)
u i = u f,i , u s,i , u z,i , u Ψ,i . (2)
For means of visualization, the time dependency of states and controls (e.g., x i (t)) is not explicitly shown here. The system state (1) consists of the quadrotor's x i , y i , z i position, its yaw angle Ψ i and the forward v f,i and sidewards velocities v s,i . Its control vector (2) concatenates forward u f,i , sidewards u s,i , upward u z,i and yaw u Ψ,i velocity. The approximation of the real quadrotor's dynamics by means of system identification, results in:

ẋi = f i (x i , u i ) =         v f,i cos (Ψ i ) -v s,i sin (Ψ i ) v f,i sin (Ψ i ) + v s,i cos (Ψ i ) u z,i 1.6u Ψ,i -0.5092v f,i + 1.458u f,i -0.5092v s,i + 1.458u s,i         . ( 3
)
The scalar parameters in (3) have been identified from a real AR.Drone 2.0 system by means of a motion capture system. This system model is used in an MPC approach, which computes optimal controls by solving the following optimal control problem (OCP):

min u i Ji = t f t 0 (x * i -xi) Q (x * i -xi) + u i Rui dτ (4) s.t. ẋi = fi (xi, ui) (5) c ≤ u f (t) 2 -1 us (t) 2 -1 uz (t) 2 -1 u 2 Ψ -1 (6)
xi (0) = 50, 0, i, 0, 0, 0

Q = D{ 1.5, 1.5, 1.6, 0.1, 0, 0 }, (8) R = D{ 1, 1, 1, 1 } (7)
for each drone i = 1...10 over a receding horizon. The system dynamics (3) are represented as constraints (5). The system inputs are defined between -1 ≤ u ≤ 1, which is taken into account with constraint (6). To avoid collisions, the quadrotors initial positions x i (0) differ in their height [START_REF] Gandomi | Chaosenhanced accelerated particle swarm optimization[END_REF]. The tracking of the desired state x * i is achieved by a quadratic state penalty with Q (8) in cost-function J i (4). For means of visibility of the area coverage, the desired trajectory waypoints in the given scenario are limited to changes in the xy-plane, while the initial z i = i and ψ i = 0 are tracked. Furthermore the quadratic control penalty with R (9) is minimizing the control effort and thus the energy consumption.

Analog to the previous work [START_REF] Dentler | A realtime model predictive position control with collision avoidance for commercial low-cost quadrotors[END_REF], a Condensed Multiple Shooting Continuation Generalized Minimal Residual (CM-SCGMRES) [START_REF] Shimizu | A real-time algorithm for nonlinear receding horizon control using multiple shooting and continuation/krylov method[END_REF] approach is used to solve the MPC with OCP (4)-( 9). The CMSCGMRES is parameterized with a maximum number of iterations k max = 6, horizon length T = 1s, control update interval of ∆T = 0.1s, n hor = 10 steps within the horizon, a forward difference step of h = 0.001s, solution tolerance of = 10 -4 , continuation factor ζ = 10 and the horizon expansion factor α = 2. More details about the parameters are given in [START_REF] Ohtsuka | A continuation/gmres method for fast computation of nonlinear receding horizon control[END_REF].

III. EXPERIMENTATION

The first part of this section provides information on the simulation framework used for the experiments. Then the experimental setup, including all parameters and metrics used are detailed. Finally the last part contains the experimental results and their analysis.

A. Simulation framework

For the experimental validation, the quadrotor dynamics (3) are implemented in the V-Rep simulation environment [START_REF] Rohmer | V-rep: a versatile and scalable robot simulation framework[END_REF]. The MPC is interfaced via ROS. An example of a scenario with three controlled UAVs is shown in Fig. 4. For the evaluation of the large area coverage with ten UAVs, as illustrated in Fig. 5, the 3D environment data is mapped onto the xy-plane as shown in Fig. 6 for 10 UAVs1 .

B. Experimental setup

The simulation is performed using the environment shown in Tab. I. The same three area coverage metrics as introduced in [9] are used. A short description is provided hereinafter. Fig. 7 underlines how the metrics are calculated for one simulation. The coverage is the percentage of the total area visited during the whole simulation. The environment is a 100 m by 100 m square area. The coverage value varies during the whole simulation. To have a representative value of the coverage, we choose to compare the coverage value after 4000 steps for each model. This indicates the efficiency of the models to visit the total area. On the other hand, we want to evaluate the first steps of each model to compare their initial behavior. We extract the slope of a linear regression a × x considering the 500 first steps (Fig. 7). 2) Fairness: The fairness measures if all cells are regularly and equally scanned. This is computed as the standard deviation of their respective number of scans [START_REF] Schleich | UAV fleet area coverage with network connectivity constraint[END_REF]. To evaluate the fairness during the whole simulation, we perform a linear regression a × x + b using the last 3500 steps (Fig. 7). This measure is complementary to the coverage initial slope that only evaluates the initial UAV trajectories. Here, the slope value is considered as a measure to evaluate the fairness of the models without this initialization part.

3) Recent coverage ratio: This metric introduced in [18] represents the percentage of coverage during the last 100 iterations. These 100 steps correspond to the pheromones' evaporation time. We exclude the 100 first iterations of our simulation to compute the mean value of the recent coverage (Fig. 7).

C. Experimental results

The previously described metrics are obtained with the same program to avoid bias. First, they are computed from the waypoints and secondly using the synchronized positions of the UAVs from the V-Rep simulation using ROS. As CACOC is a deterministic mobility model, we will only vary the initial conditions of the dynamical system to create a wide range of simulations with the same mobility model without changing its parametrization as it has already been done previously [START_REF] Rosalie | From random process to chaotic behavior in swarms of UAVs[END_REF]. Fig. 8 gives the box plot for each metric and method with a smooth density estimation. It should be noted that some of the results presented here are not comparable to our previous simulations [START_REF] Rosalie | From random process to chaotic behavior in swarms of UAVs[END_REF] because their duration is shorter due to the time necessary to compute the controlled model. Associated with the Fig. 8 we perform pairwise statistical tests. The purpose of the Shapiro-Wilks test is to evaluate if the values are normally distributed. With 95% confidence, this hypothesis is satisfied for our results (except the recent coverage of the controlled model) because the p-value ≤ 0.05.

The F-test evaluates the equality of the variance with the same mechanism for the hypothesis for a given p-value. The variances are similar for all the metrics except for the recent coverage one with 95% confidence. Finally, if the data are normally distributed with equality of the variance, we can perform a Student's test to check if the data follow the same distribution. Else, a Wilcoxon signed rank test is performed for not normally distributed data. With 95% of confidence we can say that the metric data of the Waypoints and Controlled model do not follow the same distribution. Even though, Tab. II indicates that the slope of fairness is the metric that is the most similar: the p-value of the Student's test is equal to 0.02913. We assume that the values for this metric are close because of the use of a linear regression over 87% of the simulation time. For this metric, the average value is higher for the Controlled model than for the Waypoints the average value is higher (Fig. 7). This means that the coverage is more fair for the Waypoints method than for the Controlled one. This indicates that despite some small change in the trajectories of the UAVs, the models remain globally efficient considering this metric. This is not the case for the recent coverage because the gap between the two average values represent 10% of the optimal value (see [START_REF] Rosalie | From random process to chaotic behavior in swarms of UAVs[END_REF] for details). This metric evaluates the performance on a short term giving the percentage of area covered during the last 100 steps. If we only consider this metric we can erroneously consider that the controlled method does not perform as well as the Waypoints model. Considering the whole statistics, we can conclude that the global performance of the controlled method is very satisfactory. As previously mentioned, the most noticeable performance degradation is for the recent coverage performance. Conversely, the quality of slope of fairness is very similar between the two methods. In addition, it seems that the shift between the expected position and the real one are compensated on a long term for the coverage. Due the position tracking, the quadrotor asymptotically approaches the desired target points. For moving targets as in the present work, this delay leads to a position tracking error. This explains the difference in the slope of coverage. Furthermore this position tracking error results in less observed boundary areas, as all target points lie within the given rectangular region. Accordingly, the coverage in the studied scenario is slightly lower for the controlled system. As this difference in coverage is the direct result of the closed loop dynamics of the controlled quadrotor system, this difference is supposed to decrease with an increasing size of the observed area.

IV. CONCLUSION

In this paper, we compared the coverage performance of a state-of-the-art UAV swarm mobility model (CACOC) against its physical simulation in a robotic dedicated environment. An empirical evaluation of both the initial Waypoints model and of the Controlled model was conducted on a 10 UAVs scenario. The difference between these two simulations was evaluated using four area coverage metrics and assessed with statistical tests. The obtained results show that the presented CACOC-MPC approach implemented in the robotics simulator was able to accomplish the covering task efficiently. Indeed, the high similarity of the results between the Waypoints and Controlled demonstrates that it can be expected that it will perform well in a real environment with the only consideration of a fine tuning of the MPC approach.

Future work on the UAV control part will treat the reduction of the position tracking error, for instance, with the introduction of a disturbance model, the adaptation of the MPC control law, or a target position controller as already presented in [START_REF] Dentler | A realtime model predictive position control with collision avoidance for commercial low-cost quadrotors[END_REF]. This two sided simulator will also permit to work on the CACOC parameters (pheromone evaporation time, covered area by step, etc.) to increase the performance of the mobility model by comparing the theoretical and realistic results.

Finally, we also plan to work on collision avoidance for the UAVs. In this work we consider that each UAV has its own altitude but we want to compare the performance of CACOC waypoints with controlled UAVs at the same altitude. The main question we would like to address is the following: is it statistically efficient to implement a collision avoidance mechanism for the Waypoints or the collision avoidance included in the Controlled method is sufficient?

Fig. 1 .

 1 Fig. 1. Use case: implementation structure of area exploration scenario.

Fig. 2 .

 2 Fig. 2. Details of the CACOC mobility model with the first return maps from the Rössler system and the pattern it generates when there is no pheromone to guide the UAVs. The next choice depends on the previous one with: A for ahead, L for left and R for right. The period 1 orbits leads to symbols AAAAA. . . (strong violet straight line) and the period 2 leads to the pattern ARARA. . . (lime green large circle).

Fig. 3 .

 3 Fig. 3. Quadrotor system.

Fig. 4 .

 4 Fig. 4. Three MPC controlled UAVs in the V-Rep environment following waypoints provided by CACOC for three UAVs.

Fig. 5 .

 5 Fig. 5. Simulation environment of use case of CACOC exploration with 10 MPC controlled UAVs.

Fig. 6 .

 6 Fig. 6. Ten controlled UAVs trajectories in the squared area where the thickest straight lines are links between the waypoints they have to reach.

Fig. 7 .

 7 Fig. 7. Metrics measurements for one simulation. Two values are extract for the coverage: the coverage after 4000 steps and the slope of coverage for the 500 first steps. For the fairness we compute the slope of the fairness after the 500 first steps. The average value of the recent coverage is calculated after the 100 first steps.

TABLE I MAIN

 I EXPERIMENTAL PARAMETERS.

	Parameter Name	Parameter Value
	Simulation area
	Geographical Area	100 m × 100 m
	Number of cells	100 × 100
	UAV Waypoint
	UAVs speed	1 m/s
	Possible UAV actions ahead, 45 • left, 45 • right
	Initial UAVs position	middle of the bottom of the map
	Experiments
	Mobility models	Waypoint, Controlled
	Number of UAVs	10
	Simulation steps	4000
	Independent runs	30
	1) Coverage:

TABLE II p

 II -VALUES OF TESTS TO HIGHLIGHT STATISTICAL DIFFERENCE BETWEEN Waypoints AND Controlled METHODS.

	Metric	Test	Waypoints Controlled
	Coverage after 4000 steps	Shapiro-Wilks	0.06068	0.07211
		F-test	0.3859
		Student's test	9.095e-06
	Slope of Coverage	Shapiro-Wilks	0.06752	0.4839
		F-test	0.9323
		Student's test	4.367e-08
	Slope of Fairness	Shapiro-Wilks	0.9276	0.7917
		F-test	0.7721
		Student's test	0.02913
	Recent Coverage	Shapiro-Wilks	0.5773	2.539e-05
		F-test	5.696e-07
		Wilcoxon test	≤ 2.2e-16

A video illustrating the framework simulator is available at https: //martinrosalie.gforge.uni.lu/_downloads/icuas.mp4

Acknowledgments

The experiments presented in this paper were carried out using the HPC facilities of the University of Luxembourg [19] (see http://hpc.uni.lu).

* Supported by FNR "Fonds National de la Recherche" (Luxembourg) through AFR "Aides à la Formation-Recherche" Ph.D. grant scheme No. 9312118.