C Dossal

) M Kachour

The degrees of freedom of the Lasso for general design matrix

Keywords: Lasso, model selection criteria, degrees of freedom, SURE AMS classification code: Primary 62M 10, secondary 62M 20

. We also show, by providing a simple counterexample, that although the dof theorem of [32] is correct, their proof contains a flaw since their divergence formula holds on a different set of a full measure than the one that they claim. An effective estimator of the number of degrees of freedom may have several applications including an objectively guided choice of the regularization parameter in the Lasso through the SURE framework. Our theoretical findings are illustrated through several numerical simulations.

Introduction

Problem statement

We consider the following linear regression model

y = Ax 0 + ε, µ = Ax 0 , (1)
where y ∈ R n is the observed data or the response vector, A = (a 1 , • • • , a p ) is an n × p design matrix,

x 0 = x 0 1 , • • • , x 0 p
T is the vector of unknown regression coefficients and ε is a vector of i.i.d.

centered Gaussian random variables with variance σ 2 > 0. In this paper, the number of observations n can be greater than the ambient dimension p of the regression vector to be estimated. Recall that when n < p, (1) is an underdetermined linear regression model, whereas when n ≥ p and all the columns of A are linearly independent, it is overdetermined. Let x(y) be an estimator of x 0 , and µ(y) = A x(y) be the associated response or predictor. The concept of degrees of freedom plays a pivotal role in quantifying the complexity of a statistical modeling procedure. More precisely, since y ∼ N (µ = Ax 0 , σ 2 Id n×n ) (Id n×n is the identity on R n ), according to [START_REF] Efron | How biased is the apparent error rate of a prediction rule[END_REF], the degrees of freedom (dof) of the response µ(y) is defined by

df = n i=1 cov( µ i (y), y i ) σ 2 . (2)
Many model selection criteria involve df , e.g. C p (Mallows [START_REF] Mallows | Some comments on C p[END_REF]), AIC (Akaike Information Criterion, [START_REF] Akaike | Information theory and an extension of the maximum likelihood principle[END_REF]), BIC (Bayesian Information Citerion, [START_REF] Schwarz | Estimating the dimension of a model[END_REF]), GCV (Generalized Cross Validation, [START_REF] Craven | Smoothing Noisy Data with Spline Functions: estimating the correct degree of smoothing by the method of generalized cross validation[END_REF]) and SURE (Stein's unbiased risk estimation [START_REF] Stein | Estimation of the mean of a multivariate normal distribution[END_REF], see Section 2.2). Thus, the dof is a quantity of interest in model validation and selection and it can be used to get the optimal hyperparameters of the estimator. Note that the optimality here is intended in the sense of the prediction µ(y) and not the coefficients x(y).

The well-known Stein's lemma [START_REF] Stein | Estimation of the mean of a multivariate normal distribution[END_REF] states that if y → µ(y) is weakly differentiable then its divergence is an unbiased estimator of its degrees of freedom, i.e. df (y) = div( µ(y)) = n i=1 ∂ µ i (y) ∂y i , and E( df (y)) = df .

Here, in order to estimate x 0 , we consider solutions to the Lasso problem, proposed originally in [START_REF] Tibshirani | Regression shrinkage and selection via the Lasso[END_REF]. The Lasso amounts to solving the following convex optimization problem

min x∈R p 1 2 y -Ax 2 2 + λ x 1 , (P 1 (y, λ))
where λ > 0 is called the Lasso regularization parameter and • 2 (resp. • 1 ) denotes the 2 (resp. 1 ) norm. An important feature of the Lasso is that it promotes sparse solutions. In the last years, there has been a huge amount of work where efforts have focused on investigating the theoretical guarantees of the Lasso as a sparse recovery procedure from noisy measurements. See, e.g., [START_REF] Fan | Variable selection via nonconcave penalized likelihood and its oracle properties[END_REF][START_REF] Fan | Nonconcave penalized likelihood with a diverging number of parameters[END_REF][START_REF] Zhao | On model selection consistency of Lasso[END_REF][START_REF] Zou | The adaptive Lasso and its oracle properties[END_REF][START_REF] Ravikumar | Spam: Sparse additive models[END_REF][START_REF] Nardi | On the asymptotic properties of the group Lasso estimator for linear models[END_REF][START_REF] Osborne | A new approach to variable selection in least squares problems[END_REF][START_REF] Efron | Least angle regression (with discussion)[END_REF][START_REF] Fuchs | On sparse representations in arbitrary redundant bases[END_REF][START_REF] Tropp | Just relax: convex programming methods for identifying sparse signals in noise[END_REF], to name just a few.

Contributions and related work

Let µ λ (y) = A x λ (y) be the Lasso response vector, where x λ (y) is a solution of the Lasso problem (P 1 (y, λ)). Note that all minimizers of the Lasso share the same image under A, i.e. µ λ (y) is uniquely defined; see Lemma 2 in Section 5 for details. The main contribution of this paper is first to provide an unbiased estimator of the degrees of freedom of the Lasso response for any design matrix. The estimator is valid everywhere except on a set of (Lebesgue) measure zero. We reach our goal without any additional assumption to ensure uniqueness of the Lasso solution. Thus, our result covers the challenging underdetermined case where the Lasso problem does not necessarily have a unique solution. It obviously holds when the Lasso problem (P 1 (y, λ)) has a unique solution, and in particular in the overdetermined case studied in [START_REF] Zou | On the "degrees of freedom" of the Lasso[END_REF]. Using the estimator at hand, we also establish the reliability of the SURE as an unbiased estimator of the Lasso prediction risk. While this paper was submitted, we became aware of the independent work of Tibshirani and Taylor [START_REF] Tibshirani | Degrees of Freedom in Lasso Problems[END_REF], who studied the dof for general A both for the Lasso and the general (analysis) Lasso.

Section 3 is dedicated to a thorough comparison and discussion of connections and differences between our results and the one in [START_REF] Zou | On the "degrees of freedom" of the Lasso[END_REF]Theorem 1] for the overdetermined case, and that of [START_REF] Kato | On the degrees of freedom in shrinkage estimation[END_REF][START_REF] Tibshirani | Degrees of Freedom in Lasso Problems[END_REF][START_REF] Vaiter | Robust sparse analysis regularization[END_REF] for the general case.

Overview of the paper

This paper is organized as follows. Section 2 is the core contribution of this work where we state our main results. There, we provide the unbiased estimator of the dof of the Lasso, and we investigate the reliability of the SURE estimate of the Lasso prediction risk. Then, we discuss relation of our work with concurrent one in the literature in Section 3. Numerical illustrations are given in Section 4. The proofs of our results are postponed to Section 5. A final discussion and perspectives of this work are provided in Section 6.

Main results

An unbiased estimator of the dof

First, some notations and definitions are necessary. For any vector x, x i denotes its ith component. The support or the active set of x is defined by

I = supp(x) = {i : x i = 0},
and we denote its cardinality as | supp(x)| = |I|. We denote by x I ∈ R |I| the vector built by restricting x to the entries indexed by I. The active matrix A I = (a i ) i∈I associated to a vector x is obtained by selecting the columns of A indexed by the support I of x. Let • T be the transpose symbol. Suppose that A I is full column rank, then we denote the Moore-Penrose pseudo-inverse of A I ,

A + I = (A T I A I ) -1 A T I . sign(•) represents the sign function: sign(a) = 1 if a > 0; sign(a) = 0 if a = 0; sign(a) = -1 if a < 0. For any I ⊆ {1, 2, • • • , p}, let V I = span(A I ), P V I the orthogonal projector onto V I and P V ⊥ I that onto the orthogonal complement V ⊥ I .
Let S ∈ {-1, 1} |I| be a sign vector, and j ∈ {1, 2, • • • , p}. Fix λ > 0. We define the following set of hyperplanes

H I,j,S = {u ∈ R n : P V ⊥ I (a j ), u = ±λ(1 -a j , (A + I ) T S )}. (4)
Note that, if a j does not belong to V I , then H I,j,S becomes a finite union of two hyperplanes. Now, we define the following finite set of indices

Ω = {(I, j, S) : a j ∈ V I } (5)
and let G λ be the subset of R n which excludes the finite union of hyperplanes associate to Ω, that is

G λ = R n \ (I,j,S)∈Ω H I,j,S . (6)
To cut a long story short, (I,j,S)∈Ω H I,j,S is a set of (Lebesgue) measure zero (Hausdorff dimension n -1), and therefore G λ is a set of full measure.

We are now ready to introduce our main theorem.

Theorem 1. Fix λ > 0. For any y ∈ G λ , consider M y,λ the set of solutions of (P 1 (y, λ)). Let x * λ (y) ∈ M y,λ with support I * such that A I * is full rank. Then,

|I * | = min b x λ (y)∈M y,λ | supp( x λ (y))|. (7)
Furthermore, there exists ε > 0 such that for all z ∈ Ball(y, ε), the n-dimensional ball with center y and radius ε, the Lasso response mapping z → µ λ (z) satisfies

µ λ (z) = µ λ (y) + P V I * (z -y). (8)
As stated, this theorem assumes the existence of a solution whose active matrix A I * is full rank. This can be shown to be true; see e.g. [5, Proof of Theorem 1] or [20, Theorem 3, Section B.1]1 . It is worth noting that this proof is constructive, in that it yields a solution x * λ (y) of (P 1 (y, λ)) such that A I * is full column rank from any solution x λ (y) whose active matrix has a nontrivial kernel. This will be exploited in Section 4 to derive an algorithm to get x * λ (y), and hence I * .

A direct consequence of our main theorem is that outside G λ , the mapping µ λ (y) is C ∞ and the sign and support are locally constant. Applying Stein's lemma yields Corollary 1 below. The latter states that the number of nonzero coefficients of x * λ (y) is an unbiased estimator of the dof of the Lasso.

Corollary 1. Under the assumptions and with the same notations as in Theorem 1, we have the following divergence formula df λ (y

) := div( µ λ (y)) = |I * |. (9)
Therefore,

df = E( df λ (y)) = E(|I * |). (10)
Obviously, in the particular case where the Lasso problem has a unique solution, our result holds true.

Reliability of the SURE estimate of the Lasso prediction risk

In this work, we focus on the SURE as a model selection criterion. The SURE applied to the Lasso reads

SURE( µ λ (y)) = -nσ 2 + µ λ (y) -y 2 2 + 2σ 2 df λ (y), (11)
where df (y) is an unbiased estimator of the dof as given in Corollary 1. It follows that the SURE( µ λ (y)) is an unbiased estimate of the prediction risk, i.e.

Risk(µ) = E µ λ (y) -µ 2 2 = E (SURE( µ λ (y))) .
We now evaluate its reliability by computing the expected squared-error between SURE( µ λ (y)) and SE( µ λ (y)), the true squared-error, that is

SE( µ λ (y)) = µ λ (y) -µ 2 2 . ( 12
)
Theorem 2. Under the assumptions of Theorem 1, we have

E (SURE( µ λ (y)) -SE( µ λ (y))) 2 = -2σ 4 n + 4σ 2 E µ λ (y) -y 2 2 + 4σ 4 E (|I * |) . (13)
Moreover,

E SURE( µ λ (y)) -SE( µ λ (y)) nσ 2 2 = O 1 n . ( 14
)
3 Relation to prior work

Overdetermined case [START_REF] Zou | On the "degrees of freedom" of the Lasso[END_REF] The authors in [START_REF] Zou | On the "degrees of freedom" of the Lasso[END_REF] studied the dof of the Lasso in the overdetermined case. Precisely, when n ≥ p and all the columns of the design matrix A are linearly independent, i.e. rank(A) = p. In fact, in this case the Lasso problem has a unique minimizer x λ (y) = x * λ (y) (see Theorem 1). Before discussing the result of [START_REF] Zou | On the "degrees of freedom" of the Lasso[END_REF], let's point out a popular feature of x λ (y) as λ varies in ]0, +∞[:

• For λ ≥ A T y ∞ , the optimum is attained at x λ (y) = 0.
• The interval 0, A T y ∞ is divided into a finite number of subintervals characterized by the fact that within each such subinterval, the support and the sign vector of x λ (y) are constant. Explicitly, let (λ m ) 0≤m≤K be the finite sequence of λ's values corresponding to a variation of the support and the sign of x λ (y), defined by

A T y ∞ = λ 0 > λ 1 > λ 2 > • • • > λ K = 0.
Thus, in ]λ m+1 , λ m [, the support and the sign of x λ (y) are constant, see [START_REF] Efron | Least angle regression (with discussion)[END_REF][START_REF] Osborne | A new approach to variable selection in least squares problems[END_REF][START_REF] Osborne | On the LASSO and its dual[END_REF]. Hence, we call (λ m ) 0≤m≤K the transition points. Now, let λ ∈]λ m+1 , λ m [. Thus, from Lemma 1 (see Section 5), we have the following implicit form of x λ (y), ( x λ (y))

Im = A + Im y -λ(A T Im A Im ) -1 Sm, (15)
where I m and S m are respectively the (constant) support and sign vector of x λ (y) for λ ∈]λ m+1 , λ m [. Hence, based on ( 15), [START_REF] Zou | On the "degrees of freedom" of the Lasso[END_REF] showed that for all λ > 0, there exists a set of measure zero N λ , which is a finite collection of hyperplanes in R n , and they defined

K λ = R n \ N λ , (16)
so that ∀ y ∈ K λ , λ is not any of the transition points. Then, for the overdetermined case, [START_REF] Zou | On the "degrees of freedom" of the Lasso[END_REF] stated that for all y ∈ K λ , the number of nonzero coefficients of the unique solution of (P 1 (y, λ)) is an unbiased estimator of the dof. In fact, their main argument is that, by eliminating the vectors associated to the transition points, the support and the sign of the Lasso solution are locally constant with respect to y, see [32, Lemma 5].

We recall that the overdetermined case, considered in [START_REF] Zou | On the "degrees of freedom" of the Lasso[END_REF], is a particular case of our result since the minimizer is unique. Thus, according to the Corollary 1, we find the same result as [START_REF] Zou | On the "degrees of freedom" of the Lasso[END_REF] but valid on a different set y ∈ G λ = R n \ (I,j,S)∈Ω H I,j,S . A natural question arises: can we compare our assumption to that of [START_REF] Zou | On the "degrees of freedom" of the Lasso[END_REF] ? In other words, is there a link between K λ and G λ ?

The answer is that, depending on the matrix A, these two sets may be different. More importantly, it turns out that although the dof formula [32, Theorem 1] is correct, unfortunately, their proof contains a flaw since their divergence formula [32, Lemma 5] is not true on the set K λ . We prove this by providing a simple counterexample.

Example of vectors in G λ but not in K λ Let {e 1 , e 2 } be an orthonormal basis of R 2 and let's define a 1 = e 1 and a 2 = e 1 + e 2 , and A the matrix whose columns are a 1 and a 2 .

Let's define I = {1}, j = 2 and S = 1. It turns out that A + I = a 1 and (A + I ) T S, a j = 1 which implies that for all λ > 0,

H I,j,S = {u ∈ R n : P V ⊥ I (a j ), u = 0} = span(a 1
) . Let y = αa 1 with α > 0, for any λ > 0, y ∈ H I,j,S (or equivalently here y / ∈ G λ ). Using Lemma 1

(see Section 5), one gets that for any λ ∈]0, α[, the solution of (P 1 (y, λ)) is x λ (y) = (αλ, 0) and that for any λ ≥ α, x λ (y) = (0, 0). Hence the only transition point is λ 0 = α. It follows that for λ < α, y belongs to K λ defined in [START_REF] Zou | On the "degrees of freedom" of the Lasso[END_REF], but y / ∈ G λ .

We prove then that in any ball centered at y, there exists a vector z 1 such that the support of the solution of (P 1 (z 1 , λ)) is different from the support of (P 1 (y, λ)). Let's choose λ < α and ε ∈]0, αλ[ and let's define z 1 = y + εe 2 . From Lemma 1 (see Section 5), one deduces that the solution of (P

1 (z 1 , λ)) is equal to x λ (z 1 ) = (α -λ -ε, ε) whose support is different from that of x λ (y) = (α -λ, 0).
More generally, when there are sets {I, j, S} such that (A + I ) T S, a j = 1, a difference between the two sets G λ and K λ may arise. Clearly, G λ is not only the set of transition points associated to λ.

According to the previous example, in this specific situation, for any λ > 0 there may exist some vectors y that are not transition points associated to λ where the support of the solution of (P 1 (y, λ)) is not stable to infinitesimal perturbations of y. This situation may occur for under or overdetermined problems. In summary, even in the overdetermined case, excluding the set of transition points is not sufficient to guarantee stability of the support and sign of the Lasso solution.

General case [START_REF] Kato | On the degrees of freedom in shrinkage estimation[END_REF][START_REF] Tibshirani | Degrees of Freedom in Lasso Problems[END_REF][START_REF] Vaiter | Robust sparse analysis regularization[END_REF] In [START_REF] Kato | On the degrees of freedom in shrinkage estimation[END_REF], the author studies the degrees of freedom of a generalization of the Lasso where the regression coefficients are constrained to a closed convex set. When the latter is a 1 ball and p > n, he proposes the cardinality of the support as an estimate of df but under a restrictive assumption on A under which the Lasso problem has a unique solution.

In [START_REF] Tibshirani | Degrees of Freedom in Lasso Problems[END_REF]Theorem 2], the authors proved that df = E(rank(A I ))

where I = I(y) is the active set of any solution x λ (y) to (P 1 (y, λ)). This coincides with Corollary 1 when A I is full rank with rank(A I ) = rank(A I * ). Note that in general, there exist vectors y ∈ R n where the smallest cardinality among all supports of Lasso solutions is different from the rank of the active matrix associated to the largest support. But these vectors are precisely those excluded in G λ . In the case of the generalized Lasso (a.k.a. analysis sparsity prior in the signal processing

Numerical experiments

Experiments description In this section, we support the validity of our main theoretical findings with some numerical simulations, by checking the unbiasedness and the reliability of the SURE for the Lasso. Here is the outline of these experiments.

For our first study, we consider two kinds of design matrices A, a random Gaussian matrix with n = 256 and p = 1024 whose entries are ∼ iid N (0, 1/n), and a deterministic convolution design matrix A with n = p = 256 and a Gaussian blurring function. The original sparse vector x 0 was drawn randomly according to a mixed Gaussian-Bernoulli distribution, such that x 0 is 15-sparse (i.e. | supp(x 0 ) = 15|). For each design matrix A and vector x 0 , we generate K = 100 independent replications y k ∈ R n of the observation vector according to the linear regression model [START_REF] Akaike | Information theory and an extension of the maximum likelihood principle[END_REF]. Then, for each y k and a given λ, we compute the Lasso response µ λ (y k ) using the now popular iterative soft-thresholding algorithm [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF] 2 , and we compute SURE( µ λ (y k )) and SE( µ λ (y k )). We then compute the empirical mean and the standard deviation of SURE( µ λ (y k )) 1≤k≤K , the empirical mean of SE( µ λ (y k )) 1≤k≤K , which corresponds to the computed prediction risk, and we compute R T the empirical normalized reliability on the left-hand side of (13),

R T = 1 K K k=1 SURE( µ λ (y k )) -SE( µ λ (y k )) nσ 2 2 . ( 17
)
Moreover, based on the right-hand side of (13), we compute R T as

R T = - 2 n + 4 n 2 σ 2 1 K K k=1 µ λ (y k ) -y k 2 2 + 4 n 2 1 K K k=1 (|I * | k ) , (18)
where at the kth replication, |I * | k is the cardinality of the support of a Lasso solution whose active matrix is full column rank as stated in Theorem 1. Finally, we repeat all these computations for

Avg SURE Risk

± std SURE RT (eq. ( 18)) RT (eq. ( 19)) (a) Gaussian

SURE 1 realization

Avg SURE Risk ± std SURE R T (eq. ( 18)) RT (eq. ( 19)) Construction of full rank active matrix As stated in the discussion just after Theorem 1, in situations where the Lasso problem has non-unique solutions, and the minimization algorithm returns a solution whose active matrix is rank deficient, one can construct an alternative optimal solution whose active matrix is full column rank, and then get the estimator of the degrees of freedom. More precisely, let x λ (y) be a solution of the Lasso problem with support I such that its active matrix A I has a non-trivial kernel. The construction is as follows:

(b) Convolution
1. Take h ∈ ker A I such that supp h ⊂ I.

2. For t ∈ R, A x λ (y) = A ( x λ (y) + th) and the mapping t → x λ (y) + th 1 is locally affine in a neighborhood of 0, i.e. for |t| < min j∈I |( x λ (y)) j |/ h ∞ . x λ (y) being a minimizer of (P 1 (y, λ)), this mapping is constant in a neighborhood of 0. We have then constructed a whole collection of solutions to (P 1 (y, λ)) having the same image and the same 1 norm, which lives on a segment.

3.

Move along h with the largest step t 0 > 0 until an entry of x 1 λ (y) = x λ (y) + t 0 h vanishes, i.e. supp( x 1 λ (y) + t 0 h) I.

4. Repeat this process until getting a vector x * λ (y) with a full column rank active matrix A I * .

Note that this construction bears similarities with the one in [START_REF] Rosset | Boosting as a Regularized Path to a Maximum Margin Classifier[END_REF]. RT (eq. ( 18))

Upper-bound (eq. ( 49))

Avg SURE Risk 5% quantiles n n RT (eq. ( 19))

(a) λ/σ = 0.1 RT (eq. ( 18))

Upper-bound (eq. ( 49))

Avg SURE Risk 5% quantiles n n RT (eq. ( 19))

(b) λ/σ = 1 RT (eq. ( 18)) RT (eq. ( 19))

Upper-bound (eq. ( 49))

Avg SURE Risk 5% quantiles n n (c) λ/σ = 10
Figure 3: The SURE and its reliability as a function of the number of observations n.

Results discussion Figure 2 depicts the obtained results. For each design matrix, we associate a panel, each containing three plots. Hence, for each case, from left to right, the first plot represents the SURE for one realization of the noise as a function of λ. In the second graph, we plot the computed prediction risk curve and the empirical mean of the SURE as a function of the regularization parameter λ. Namely, the dashed curve represents the calculated prediction risk, the solid curve represents the empirical mean of the SURE, and the shaded area represent the empirical mean of the sure ± the empirical standard deviation of the SURE. The latter shows that the SURE is an unbiased estimator of the prediction risk with a controlled variance. This suggests that the SURE is consistent, and then so is our estimator of the degrees of freedom. In the third graph, we plot the theoretical and empirical normalized reliability, defined respectively by ( 17) and ( 18), as a function of the regularization parameter λ. More precisely, the solid and dashed blue curves represent respectively R T and R T . This confirms numerically that both sides (R T and R T ) of ( 13) indeed coincide.

As discussed in the introduction, one of the motivations of having an unbiased estimator of the degrees of freedom of the Lasso is to provide a data-driven objective way for selecting the optimal Lasso regularization parameter λ. For this, one can compute the optimal λ that minimizes the SURE, i.e.

λ optimal = argmin λ>0 SURE( µ λ (y)). (19)
In practice, this optimal value can be found either by a exhaustive search over a fine grid, or alternatively by any dicothomic search algorithm (e.g. golden section) if λ → SURE( µ λ (y)) is unimodal. Now, for our second simulation study, we consider a partial Fourier design matrix, with n < p and a constant underdeterminacy factor p/n = 4. x 0 was again simulated according to a mixed Gaussian-Bernoulli distribution with 0.1p non-zero entries. For each of three values of λ/σ ∈ {0.1, 1, 10} (small, medium and large), we compute the prediction risk curve, the empirical mean of the SURE, as well as the values of the normalized reliability R T and R T , as a function of n ∈ {8, • • • , 1024}. The obtained results are shown in Figure 3. For each value of λ, the first plot (top panel) displays the normalized empirical mean of the SURE (solid line) and its 5% quantiles (dotted) as well as the computed normalized prediction risk (dashed). Unbiasedness is again clear whatever the value of λ. The trend on the prediction risk (and average SURE) is in agreement with rates known for the Lasso, see e.g. [START_REF] Bickel | Simultaneous analysis of Lasso and Dantzig selector[END_REF]. The second plot confirms that the SURE is an asymptotically reliable estimate of the prediction risk with the rate established in Theorem 2. Moreover, as expected, the actual reliability gets closer to the upper-bound (48) as the number of samples n increases.

Proofs

First of all, we recall some classical properties of any solution of the Lasso (see, e.g., [START_REF] Osborne | A new approach to variable selection in least squares problems[END_REF][START_REF] Efron | Least angle regression (with discussion)[END_REF][START_REF] Fuchs | On sparse representations in arbitrary redundant bases[END_REF][START_REF] Tropp | Just relax: convex programming methods for identifying sparse signals in noise[END_REF]). To lighten the notation in the two following lemmas, we will drop the dependency of the minimizers of (P 1 (y, λ)) on either λ or y.

Lemma 1. x is a (global) minimizer of the Lasso problem (P 1 (y, λ)) if and only of:

1. A T I (y -A x) = λ sign( x I )
, where I = {i : x i = 0}, and

2. | a j , y -A x | ≤ λ, ∀ j ∈ I c ,
where I c = {1, . . . , p}\I. Moreover, if A I is full column rank, then x satisfies the following implicit relationship:

x I = A + I y -λ(A T I A I ) -1 sign( x I ) . (20)
Note that if the inequality in condition 2 above is strict, then x is the unique minimizer of the Lasso problem (P 1 (y, λ)) [START_REF] Fuchs | On sparse representations in arbitrary redundant bases[END_REF].

Lemma 2 below shows that all solutions of (P 1 (y, λ)) have the same image by A. In other words, the Lasso response µ λ (y), is unique, see [START_REF] Dossal | A necessary and sufficient condition for exact recovery by l1 minimization[END_REF].

Lemma 2. If x 1 and x 2 are two solutions of (P 1 (y, λ)), then

A x 1 = A x 2 = µ λ (y).
Before delving into the technical details, we recall the following trace formula of the divergence. Let J b µ(y) be the Jacobian matrix of a mapping y → µ(y), defined as follows

J b µ(y) i,j := ∂ µ(y) i ∂y j , i, j = 1, • • • , n. (21)
Then we can write div ( µ(y)) = tr

J b µ(y) . (22)
Proof of Theorem 1. Let x * λ (y) be a solution of the Lasso problem (P 1 (y, λ)) and I * its support such that A I * is full column rank. Let (x * λ (y)) I * be the restriction of x * λ (y) to its support and S * = sign ((x * λ (y)) I * ). From Lemma 2 we have,

µ λ (y) = Ax * λ (y) = A I * (x * λ (y)) I * .
According to Lemma 1, we know that

A T I * (y -µ λ (y)) = λS * ; | a k , y -µ λ (y) | ≤ λ, ∀ k ∈ (I * ) c .
Furthermore, from [START_REF] Rosset | Boosting as a Regularized Path to a Maximum Margin Classifier[END_REF], we get the following implicit form of x * λ (y)

(x * λ (y)) I * = A + I * y -λ(A T I * A I * ) -1 S * . (23)
It follows that

µ λ (y) = P V I * (y) -λd I * ,S * , (24)
and

r λ (y) = y -µ λ (y) = P V ⊥ I * (y) + λd I * ,S * , (25)
where d I * ,S * = (A + I * ) T S * . We define the following set of indices

J = {j : | a j , r λ (y) | = λ}. (26)
From Lemma 1 we deduce that I * ⊂ J.

Since the orthogonal projection is a self-adjoint operator and from [START_REF] Tibshirani | Degrees of Freedom in Lasso Problems[END_REF], for all j ∈ J, we have

| P V ⊥ I * (a j ), y + λ a j , d I * ,S * | = λ. (27)
As y ∈ G λ , we deduce that if j ∈ J ∩ (I * ) c then inevitably we have a j ∈ V I * , and therefore | a j , d I * ,S * | = 1.

In fact, if a j ∈ V I * then (I * , j, S * ) ∈ Ω and from [START_REF] Tropp | Just relax: convex programming methods for identifying sparse signals in noise[END_REF] we have that y ∈ H I * ,j,S * , which is a contradiction with y ∈ G λ . Therefore, the collection of vectors (a i ) i∈I * forms a basis of V J = span(a j ) j∈J . Now, suppose that x λ (y) is another solution of (P 1 (y, λ)), such that its support I is different from I * . If A I is full column rank, then by using the same above arguments we can deduce that (a i ) i∈I forms also a basis of V J . Therefore, we have

|I| = |I * | = dim(V J ).
On the other hand, if A I is not full rank, then there exists a subset I 0 I such that A I0 is full rank (see the discussion following Theorem 1) and (a i ) i∈I0 forms also a basis of V J , which implies that

|I| > |I 0 | = dim(V J ) = |I * |.
We conclude that for any solution x λ (y) of (P 1 (y, λ)), we have

| supp( x λ (y))| ≥ |I * |,
and then |I * | is equal to the minimum of the cardinalities of the supports of solutions of (P 1 (y, λ)). This proves the first part of the theorem.

Let's turn to the second statement. Note that G λ is an open set and all components of (x * λ (y)) I * are nonzero, so we can choose a small enough ε such that Ball(y, ε) G λ , that is, for all z ∈ Ball(y, ε), z ∈ G λ . Now, let x 1 λ (z) be the vector supported in I * and defined by

(x 1 λ (z)) I * = A + I * z -λ(A T I * A I * ) -1 S * = (x * λ (y)) I * + A + I * (z -y). ( 29
)
If ε is small enough, then for all z ∈ Ball(y, ε), we have

sign(x 1 λ (z)) I * = sign(x * λ (y)) I * = S * . (30)
In the rest of the proof, we invoke Lemma 1 to show that, for ε small enough, x 1 λ (z) is actually a solution of (P 1 (z, λ)). First we notice that z -Ax 1 λ (z) = P V ⊥ I (z) + λd I * ,S * . It follows that

A T I * (z -Ax 1 λ (z)) = λA T I * d I * ,S * = λS * = λ sign (x 1 λ (z)) I * . (31)
Moreover for all j ∈ J ∩ I * , from (28), we have that

| a j , z -Ax 1 λ (z) | = | a j , P V ⊥ I * (z) + λd I * ,S * | = | P V ⊥ I * (a j ), z + λ a j , d I * ,S * | = λ| a j , d I * ,S * | = λ. and for all j / ∈ J | a j , z -Ax 1 λ (z) | ≤ | a j , y -Ax * λ (y) | + | P V ⊥ I * (a j ), z -y |
Since for all j / ∈ J, | a j , y -Ax * λ | < λ, there exists ε such that for all z ∈ Ball(y, ε) and ∀ j / ∈ J,

we have | a j , z -Ax 1 λ (z) | < λ.
Therefore, we obtain

| a j , z -Ax 1 λ (z) | ≤ λ, ∀ j ∈ (I * ) c .
Which, by Lemma 1, means that x 1 λ (z) is a solution of (P 1 (z, λ)), and the unique Lasso response associated to (P 1 (z, λ)), denoted by µ λ (z), is defined by

µ λ (z) = P V I * (z) -λd I * ,S * . (32)
Therefore, from ( 24) and ( 32), we can deduce that for all z ∈ Ball(y, ε) we have

µ λ (z) = µ λ (y) + P V I * (z -y).
Proof of Corollary 1. We showed that there exists ε sufficiently small such that

z -y 2 ≤ ε ⇒ µ λ (z) = µ λ (y) + P V I * (z -y). (33)
Let h ∈ V I * such that h 2 ≤ ε and z = y + h. Thus, we have that zy 2 ≤ ε and then

µ λ (z) -µ λ (y) 2 = P V I * (h) 2 = h 2 ≤ ε. (34)
Therefore, the Lasso response µ λ (y) is uniformly Lipschitz on G λ . Moreover, µ λ (y) is a continuous function of y, and thus µ λ (y) is uniformly Lipschitz on R n . Hence, µ λ (y) is almost differentiable; see [START_REF] Meyer | On the degrees of freedom in shape restricted regression[END_REF] and [START_REF] Efron | Least angle regression (with discussion)[END_REF].

On the other hand, we proved that there exists a neighborhood of y, such that for all z in this neighborhood, there exists a solution of the Lasso problem (P 1 (z, λ)), which has the same support and the same sign as x * λ (y), and thus µ λ (z) belongs to the vector space V I * , whose dimension equals to |I * |, see [START_REF] Tibshirani | The Solution Path of the Generalized Lasso[END_REF] and [START_REF] Zou | On the "degrees of freedom" of the Lasso[END_REF]. Therefore, µ λ (y) is a locally affine function of y, and then

This holds almost everywhere since G λ is of full measure, and ( 10) is obtained by invoking Stein's lemma.

Proof of Theorem 2. First, consider the following random variable Q 1 ( µ λ (y)) = µ λ (y) 2 2 + µ 2 2 -2 y, µ λ (y) + 2σ 2 div( µ λ (y)).

From Stein's lemma, we have E ε, µ λ (y) = σ 2 E (div( µ λ (y))) .

Thus, we can deduce that Q 1 ( µ λ (y)) and SURE( µ λ (y)) are unbiased estimator of the prediction risk, i.e.

E (SURE( µ λ (y))) = E (Q 1 ( µ λ (y))) = E (SE( µ λ (y))) = Risk(µ).

Moreover, note that SURE( µ λ (y)) -Q 1 ( µ λ (y)) = y

Discussion

In this paper we proved that the number of nonzero coefficients of a particular solution of the Lasso problem is an unbiased estimate of the degrees of freedom of the Lasso response for linear regression models. This result covers both the over and underdetermined cases. This was achieved through a divergence formula, valid almost everywhere except on a set of measure zero. We gave a precise characterization of this set, and the latter turns out to be larger than the set of all the vectors associated to the transition points considered in [START_REF] Zou | On the "degrees of freedom" of the Lasso[END_REF] in the overdetermined case. We also highlight the fact that even in the overdetermined case, the set of transition points is not sufficient for the divergence formula to hold. We think that some techniques developed in this article can be applied to derive the degrees of freedom of other nonlinear estimating procedures. Typically, a natural extension of this work is to consider other penalties such as those promoting structured sparsity, e.g. the group Lasso.

Figure 1 :

 1 Figure 1: A counterexample for n = p = 2 of vectors in G λ but not in K λ . See text for a detailed discussion.

Figure 2 :

 2 Figure 2: The SURE and its reliability as a function of λ for two types of design matrices. (a) Gaussian; (b) Convolution. For each kind of design matrix, we associate three plots.

  y) = P V I * . (35)Then the trace formula[START_REF] Schwarz | Estimating the dimension of a model[END_REF] implies that div ( µ λ (y)) = tr (P V I * ) = |I * |.

  Then, since |I * | = O(n) and from (47), we haveE SURE( µ λ (y)) -SE( µ λ (y))

		nσ 2	2	≤	6 n	+	4 µ 2 2 n 2 σ 2 .	(48)
	Finally, since µ 2 < +∞, we can deduce that
	E	SURE( µ λ (y)) -SE( µ λ (y)) nσ 2	2		= O	1 n	.
		2 2 -E y 2 2 , where
	E y 2 2 = nσ 2 + µ 2 2 , and V y 2 2 = 2σ 4 n + 2	µ 2 2 σ 2	.	(37)

This proof is alluded to in the note at the top of[START_REF] Sardy | Block coordinate relaxation methods for nonparametric wavelet denoising[END_REF] Page 363].

Iterative soft-thresholding through block-coordinate relaxation was proposed in[START_REF] Sardy | Block coordinate relaxation methods for nonparametric wavelet denoising[END_REF] for matrices A structured as the union of a finite number of orthonormal matrices.

Acknowledgement This work was partly funded by the ANR grant NatImages, ANR-08-EMER-009.

Now, we remark also that Q 1 ( µ λ (y)) -SE( µ λ (y)) = 2 σ 2 div( µ λ (y))ε, µ λ (y) .

(38)

After an elementary calculation, we obtain

where

with

and

Hence, by using the fact that a Gaussian probability density ϕ(ε i )

and integrations by parts, we find that

Moreover, from [13, Property 1], we know that

Thus, since J b µ λ (y) = P V I * which is an orthogonal projector (hence self-adjoint and idempotent), we have tr

Furthermore, observe that

Therefore, by combining (37), (39), ( 43) and (45), we obtain

On the other hand, since x * λ (y) is a minimizer of the Lasso problem (P 1 (y, λ)), we observe that

Therefore, we have

(47)


B Fabrèges

B Maury

Approximation of single layer distributions by Dirac masses in Finite Element computations

We are interested in the finite element solution of elliptic problems with a right-hand side of the single layer distribution type. Such problems arise when one aims at accounting for a physical hypersurface (or line, for bi-dimensional problem), but also in the context of fictitious domain methods, when one aims at accounting for the presence of an inclusion in a domain (in that case the support of the distribution is the boundary of the inclusion). The most popular way to handle numerically the single layer distribution in the finite element context is to spread it out by a regularization technique. An alternative approach consists in approximating the single layer distribution by a combination of Dirac masses. As the Dirac mass in the right hand side does not make sense at the continuous level, this approach raises particular issues. The object of the present paper is to give a theoretical background to this approach. We present a rigorous numerical analysis of this approximation, and we present two examples of application of the main result of this paper. The first one is a Poisson problem with a single layer distribution as a right-hand side and the second one is another Poisson problem where the single layer distribution is the Lagrange multiplier used to enforce a Dirichlet boundary condition on the boundary of an inclusion in the domain. Theoretical analysis is supplemented by numerical experiments in the last section.

Introduction

We are interested in the numerical handling of elliptic problems with a righthand side of the single layer distribution type, i.e. problems of the form 1

-∆u = ϕδ γ in Ω, u = 0 on ∂Ω
where Ω is a bounded domain in R d , and ϕδ γ , in H -1 (Ω), can be formally written as

ϕδ γ , v H -1 , H 1 0 = γ ϕv
where γ is a d -1 dimensional smooth manifold (typically the boundary of a connected subdomain ω ⊂⊂ Ω) and ϕ is in a Sobolev space H -1/2+s (γ) for s > 0. This type of problem is commonly encountered in numerical modeling. Peskin and Mc Queen [START_REF] Peskin | A three-dimensional computational method for blood flow in the heart I. Immersed elastic fibers in a viscous incompressible fluid[END_REF] modeled the motion of the heart wall by a collection of elastic fiber immersed in a fluid. The force applied to the fluid by these fibers is a single layer distribution. In a similar context, Gerbeau et al. [START_REF] Santos | A partitioned fluidstructured algorithm for elastic thin valves with contact[END_REF] modeled more recently a valve by an elastic interface immersed in a fluid. Likewise, Cottet and Maitre in [START_REF] Cottet | A level set method for fluid-structure interactions with immersed surfaces[END_REF] track the interface of an elastic membrane immersed in a fluid by a level set method and the elastic forces on the membrane appears in the right-hand side of the Navier-Stokes equations. In [START_REF] Herrmann | A balanced force refined level set grid method for twophase flows on unstructured flow solver grids[END_REF] the interface between two different fluids is tracked with a level set method and the single layer distribution is the surface tension.

This kind of problem also appears in the context of fictitious domain methods. For examples in [START_REF] Glowinski | A fictitious domain approach to the direct numerical simulation of incompressible viscous flow past moving rigid bodies : application to particulate flow[END_REF] and [START_REF] Bertrand | A three-dimensional fictitious domain method for incompressible fluid flow problems[END_REF] the authors use Lagrange multipliers to enforce Dirichlet boundary conditions. These Lagrange multipliers act as singular distributions of forces supported by the boundary. In [START_REF] Pan | A Lagrange multipliers / fictitious domain / collocation method for solid-liquid flows[END_REF], such an approach is proposed to simulate the sedimentation of rigid particles in a fluid. The rigid body constraint is enforced with Lagrange multipliers supported on the boundary of the rigid bodies. The single layer distribution also appears in the so called Fat Boundary Method (see [START_REF] Maury | A Fat Boundary Method for the Poisson Problem in a Domain with Holes[END_REF], or [START_REF] Bertoluzza | Analysis of the fully discrete fat boundary method[END_REF] for a full analysis of the method). The resolution of a Poisson like problem in a domain with holes with this method consists of splitting the Poisson problem into two new problems: a global problem, which is solved on a regular mesh, and a local problem around the holes, which introduces the single layer distribution on the boundary of the holes.

There exist several approaches to solve these problems numerically, the choice of one of them depending on the general framework of the discrete problem. In [START_REF] Peskin | A three-dimensional computational method for blood flow in the heart I. Immersed elastic fibers in a viscous incompressible fluid[END_REF], the authors have chosen to solve the problem in a finite difference framework. They take a lattice on the whole domain to solve the fluid part of the problem and a uniform collection of points on the fibers, which are not necessarily the same as the points of the lattice. In order to take into account the force of the fibers to the fluid, they regularize the Dirac Delta function on every points on the fibers to spread the single layer distribution out on the computational lattice. This is the so called Immersed Boundary (IB) method. Applications of the IB method can be found in section 9 of [START_REF] Peskin | The immersed boundary method[END_REF] and, for instance, in [START_REF] Givelberg | A comprehensive three-dimensional model of the cochlea[END_REF][START_REF] Mittal | Immersed Boundary Methods[END_REF][START_REF] Tyson | Modelling anguilliform swimming at intermediate Reynolds number : A review and a novel extension of immersed boundary method applications[END_REF][START_REF] Kim | 3-D Parachute simulation by the immersed boundary method[END_REF] for more recent ones. Another way to deal with the Dirac Delta functions is to directly inject the jump of the normal derivatives of the corresponding solution in the finite difference scheme. This method, called the Immersed Interface Method (IIM), has been introduced by Leveque and Li in [START_REF] Leveque | The Immersed Interface Method for elliptic equations with discontinuous coefficients and singular sources[END_REF]. For applications of the IIM, one can look in [START_REF] Tan | An immersed interface method for solving incompressible viscous flows with piecewise constant viscosity across a moving elastic membrane[END_REF] and the references therein.

In a variational context, one can compute the integrals involving the single layer distribution of the variational formulation in some ways. It is the case in a finite element context as in [START_REF] Glowinski | A fictitious domain approach to the direct numerical simulation of incompressible viscous flow past moving rigid bodies : application to particulate flow[END_REF] where the authors use a constant piece-wise element to discretize the Lagrange multiplier space and then compute the integrals exactly. The IB method has been adapted to the finite element framework where one can deal with the Dirac Delta function in a variational way (see [START_REF] Boffi | A finite element approach to the immersed boundary method[END_REF][START_REF] Ilinca | A finite element immersed boundary method for fluid flows around rigid objects[END_REF]). In the same spirit, the idea is to discretize the single layer distribution by using a combination of Dirac masses at a collection of points describing the interface. That is the method used in [START_REF] Santos | A partitioned fluidstructured algorithm for elastic thin valves with contact[END_REF][START_REF] Pan | A Lagrange multipliers / fictitious domain / collocation method for solid-liquid flows[END_REF][START_REF] Bertrand | A three-dimensional fictitious domain method for incompressible fluid flow problems[END_REF] to approximate the Lagrange multiplier defined on an interface of the computational domain. This method requires no additional mesh but only the collection of points on the interface and is convenient to implement in a code. Moreover the integrals which appear in the variational formulation are fast and easy to compute.

Up to our knowledge, the latter approach has not been justified from a theoretical point of view. One of the difficulties is that a well defined linear functional, the initial right-hand side which we will consider in H -1 , is replaced by a combination of Dirac masses, which do not make sense at the continuous level (as soon as the dimension d is greater than 1).

The aim of this paper is to give a numerical analysis of this method, in the case of a scalar Poisson problem with a single layer distribution as a right hand side in the two dimensional setting. This paper is structured as follow: In section 2 we establish an error estimate when ϕ is in H -1/2+s with 0 < s ≤ 1/2 and show that there is a saturation of the order once s is greater than 0.5 (when ϕ is at least in L 2 ). Then, in section 3 we present some examples on which this result can be used to estimate the error. In section 4, we present the numerical results to validate the error analysis done in section 2.

General theorems

Let Ω be a domain in R 2 and let ω be a smooth subdomain with boundary γ. We consider a single layer distribution supported by γ:

ϕδ γ ∈ H -1 (Ω) : v ∈ V = H 1 0 (Ω) -→ ϕ, v H -1/2+s (γ), H 1/2-s (γ) ,
with ϕ ∈ H -1/2+s (γ), 0 ≤ s < 1. The expression above makes sense as v ∈ H 1 0 (Ω), so that its trace on γ is in H 1/2 (γ) ⊂ H 1/2-s . Given a conforming triangulation T h of Ω, we denote by V h the associated P 1 finite element space:

V h = v h ∈ C 0 (Ω) , v h|K is affine , ∀K ∈ T h .
We are interested in approximating ϕδ γ over V h by a combination of Dirac masses. Note that it does not make sense at the continuous level, as V ′ = H -1 (Ω) contains no Dirac mass. Approximation properties of such an appropriate combination will be expressed by Proposition 1. It is mainly based on the following two lemmas. Lemma 1. Let I be the unit interval (0, 1), h > 0, and (T h ) a family of quasi-uniform triangulations of I. We shall represent T h by its subintervals γ 1 , . . . , γ N (we drop the explicit dependence of

γ i upon h). Quasi-uniformity expresses ch ≤ |γ i | ≤ Ch , with 0 < c < C.
We now consider x 1 , . . . , x N , with x i ∈ γ i . For any v ∈ H 1 (I), we denote by v h the corresponding piecewise constant interpolant (see Fig. 1)

v h = N i=1 v(x i )1 γ i .
For any r, 0 ≤ r < 1/2, one has

v -v h H r ≤ Ch 1-r |v| 1
where |v| 1 is the H 1 seminorm, and the fractional derivative Sobolev norm is defined by

w H r = I |w| 2 1/2 + I I |w(y) -w(x)| 2 |y -x| 1+2r 1/2 , 0 < r < 1/2.
Proof. Let us start with the integral over

I × I. Setting w h = v -v h , one has I I |w h (y) -w h (x)| 2 |y -x| 1+2r = i j γ i γ j |v(y) -v(x j ) -v(x) + v(x i )| 2 |y -x| 1+2r = i γ i γ i |v(y) -v(x)| 2 |y -x| 1+2r + i j =i γ i γ j |w h (y) -w h (x)| 2 |y -x| 1+2r = A + B. x 1 x 2 x N Figure 1: P 0 interpolant
The first term (diagonal terms) writes

A = i γ i γ i |v(y) -v(x)| 2 |y -x| 1+2r ≤ i γ i γ i γ i |v ′ (t)| 2 dt |y -x| 2r ≤ 1 (1 -2r) h 2-2r I v ′ (t) 2 .
As for extradiagonal terms, let us first note that, on any γ i ,

|w h | 2 = x x i v ′ (t) dt 2 ≤ h γ i v ′ (t) 2 .
The second term B can then be estimated

B ≤ 4 i j =i γ i γ j |w h (x)| 2 |y -x| 1+2r ≤ 4h i γ i v ′ (t) 2 j =i γ i γ j 1 |y -x| 1+2r .
In the case r > 0, the last sum above (over j = i) is less than twice

h 0 dx 1 h dy |y -x| 1+2r = - 1 2r h 0 dx (y -x) -2r 1 h ≤ 1 2r h 0 (h -x) -2r dx ≤ 1 2r(1 -2r) h 1-2r , so that B ≤ 4 r(1 -2r) h 2-2r I v ′ (t) 2 .
Thus, one obtains

I I |w h (y) -w h (x)| 2 |y -x| 1+2r ≤ C r(1 -2r) h 2-2r |v| 2 1,I .
The 0-th order term is simply

I |w h | 2 = i γ i |w h | 2 ≤ h 2 |v| 2 1,I .
Hence the order 1 -r of the lemma.

Remarks

• If r = 0, the H r norm is only the L 2 norm and the result comes from the 0-th order term of the proof.

• This result may be ruled out in the case -1/2 < r < 0. Indeed, even though the semi-norm part of the definition of the H r norm can be extended for r < 0, it is not the case for the L 2 norm. One cannot get rid of this L 2 norm either, because the constant functions must be controlled in a way. Moreover, the integral over I × I is of order h 2-2r even if r < 0, so it must also be true for the expression that replaces the L 2 norm and controls the constant functions. For instance, taking the average of the function is not sufficient because it is still a term of order 1. That is why it may be that the order cannot go further than 1, even for r < 0.

The second lemma is a trace theorem for discrete functions.

Lemma 2. Let Ω be the unit square in R 2 , and ω ⊂⊂ Ω a smooth subdomain, with boundary γ. Let (T h ) h be a regular family of triangulations, and V h the associated P 1 Finite Element space. The trace operator γ 0 maps V h (equipped with

H 1 norm) onto H 1 (γ), with |γ 0 (v h )| 1,γ ≤ C h 1/2 |v h | 1,Ω .
Proof. Because v h is a P 1 function, its gradient is constant in every mesh cell of T h . Let Q γ be the set of all the mesh cells which have a non empty intersection with γ, one has

|γ 0 (v h )| 2 1,γ = γ |∇v h | 2 = Q∈Qγ |∇v h | 2 ∞,Q |Q ∩ γ| = 1 h 2 Q∈Qγ |v h | 2 1,Q |Q ∩ γ| .
Assuming that the radius of curvature of γ is greater than h, the following inequality holds,

|Q ∩ γ| ≤ Ch ,
where C is a constant independent of h. Thus, the previous equality writes,

|γ 0 (v h )| 2 1,γ ≤ C h Q∈Qγ |v h | 2 1,Q ≤ C h |v h | 2 1,Ω ,
which is the estimation of the lemma.

The next proposition is the approximation error of a single layer distribution by a combination of Dirac masses.

Proposition 1.

Let Ω be the unit square in R 2 , and ω ⊂⊂ Ω a smooth subdomain, with boundary γ. Let (T h ) h be a regular family of triangulations, and V h the associated P 1 Finite Element space. Let 0 < s ≤ 1/2 and ϕ ∈ H -1/2+s (γ).

Let (S h) be a family of quasi-uniform triangulations of γ. We shall represent S h by its subintervals γ 1 , . . . , γ N h .

We now consider x 1 , . . . , x N h , with x i ∈ γ i and the following approximation of ϕ,

ϕ h h = N h i=1 λ i δ x i ( 1
)
where the real numbers λ i are chosen as follows,

λ i = ϕ, 1 γ i . (2)
There exists a constant C > 0 such that :

ϕ, v h -ϕ h h , v h ≤ C h h hs ϕ -1/2+s,γ |v h | 1,Ω (3)
for all function v h in V h and where . , . denotes the dual pairing between H -1/2+s (γ) and H 1/2-s (γ).

Proof. Thanks to the definition of ϕ h h , one has,

ϕ, v h -ϕ h h , v h = ϕ, v h -ϕ, N h i=1 v h (x i )1 γ i ≤ ϕ -1/2+s,γ v h - N h i=1 v h (x i )1 γ i 1/2-s,γ .
The function ṽh defined by ṽh =

N h i=1 v h (x i )1 γ i
is the piecewise constant interpolant of v h on γ. Now, in order to use lemma 1 we perform a change of variable to map the curve γ onto [0, 1],

X : [0, 1] -→ γ t -→ X(t) .
The variable t here is the arc length so X(t) is the point on γ such that the length of the arc X(0)X(t) is equal to |γ| t. Thus, denoting by w h h the function v h -ṽh , one has,

w h h 2 1/2-s,γ = γ γ w h h (y) -w h h (x) 2 |y -x| 2-2s = 1 0 1 0 w h h • X(t) -w h h • X(τ ) 2 |X(t) -X(τ )| 2-2s X ′ (t)X ′ (τ )dtdτ .
For any t in [0, 1], X ′ (t) is equal to |γ| so the change of variables writes,

w h h 1/2-s,γ = |γ| w h h • X 1/2-s,[0,1]
.

The same goes for the L 2 norm. Now, the length of the subintervals γ i of γ is h, so the length of the corresponding subintervals in [0, 1] is h/ |γ|. We apply the lemma 1 with r = 1/2 -s and h = h/ |γ|,

v h -ṽh 1/2-s,γ = |γ| v h • X -ṽh • X 1/2-s,[0,1] ≤ |γ| C |γ| 1/2+s h1/2+s |v h • X| 1,[0,1] .
Again, with the change of variables, one has,

v h -ṽh 1/2-s,γ ≤ C |γ| 1/2+s h1/2+s |v h | 1,γ .
Now, the trace lemma 2 allows us to extend the H 1 semi-norm on γ of v h over the whole domain Ω,

v h -ṽh 1/2-s,γ ≤ C |γ| 1/2+s h1/2+s h 1/2 |v h | 1,Ω ,
which ends the proof.

Remark 1 : If 1/2 < s ≤ 1, one cannot expect to have a better order than the one with s = 1/2. Indeed, here is a one dimensional example where ϕ is a smooth function and the error is still of order 1.

The domain is (0, 1) and ϕ is the constant function 1. We take h = h/2 and the function v h is defined by

v h (x) =      - x -h/2 h + 1 in [0, h] 0 elsewhere.
If we take the point h/2 to approximate ϕ in [0, h], the error writes,

ϕ, v h -ϕ h h , v h = h 2 - h 0 ϕ v h h 2 = h 2 .
Thus the error is of order one and there is a saturation of the order once s is greater than 0.5.

Remark 2 : In the case s = 0 the definition of the approximation of the single layer distribution is different because the characteristic functions are not in H 1/2 . Nevertheless, these characteristic functions can be replaced by a partition of unity and the result of Proposition 1 holds. Thus, there is still convergence if h is such that h/h goes to 0 as h goes to 0.

Applications

Error estimate for the Poisson problem

We consider here a Poisson problem with homogeneous Dirichlet boundary conditions and a single layer distribution as a right-hand side. We are interested in the error estimate when the single layer distribution is approximated by a combination of Dirac masses.

Proposition 2.

Let Ω be the unit square and ω ⊂⊂ Ω a smooth subdomain, with boundary γ. We consider the following Poisson problem,

-∆u = ϕδ γ in Ω u = 0 on ∂Ω where ϕ is in H -1/2+s (γ), 0 < s ≤ 1/2.
Let (T h ) h be a regular family of triangulations and V h the associated P 1 Finite Element space. We approximate ϕ by a function ϕ h h according to Eqs. ( 1) and [START_REF] Peskin | A three-dimensional computational method for blood flow in the heart I. Immersed elastic fibers in a viscous incompressible fluid[END_REF]. Denoting by u h the finite element solution of u and assuming that h is of the order of h, the error estimate is

u -u h 1 ≤ C √ h|u| 2 + h s ϕ -1/2+s,γ ,
where |.| 2 denotes the H 2 semi-norm.

Proof. Thanks to Strang's lemma (see [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF]), the error estimate is,

u -u h ≤ C   inf v h ∈V h u -v h + sup w h ∈V h | ϕ, w h -ϕ h h , w h | w h   .
The first part is the usual error and is of order 1/2 in our case because T h is a non conformal mesh and V h is the Finite Element space of P 1 function. The second part is given by Proposition 1. Thus, the error estimate writes,

u -u h 1 ≤ C   √ h|u| 2 + h h hs ϕ -1/2+s,γ   .
In particular, if h is equal to h, the error estimate becomes,

u -u h 1 ≤ C √ h|u| 2 + h s ϕ -1/2+s,γ ,
which ends the proof.

Error estimate for the approximation of the solution of a saddle-point problem

We consider a Poisson equation in a domain with a hole and homogeneous Dirichlet boundary conditions on each boundary. The homogeneous condition on the boundary of the hole is enforced by Lagrange multipliers defined on the boundary. These Lagrange multipliers are approximated by a combination of Dirac masses. First we prove an abstract theorem on the error approximation when the discrete Lagrange multipliers space is not included in the continuous Lagrange multipliers space. Second, we use this theorem and Proposition 1 to give an error estimate of the Poisson problem.

Proposition 3. Let V be a Hilbert space, a be a bounded elliptic bilinear form in V and f be in V ′ . Let Λ be another Hilbert space, B ∈ L(V, Λ) and K the kernel of B. We consider the problem of finding u in K such that a(u, v) = f, v for all v in K, and its sadle-point formulation

a(u, v) + B * λ, v = f, v ∀v ∈ V (µ, Bu) = 0 ∀µ ∈ Λ .
Let V h be a finite dimensional subspace of V and Λ h be a finite dimensional space, not necessarily included in Λ. Let B h h be a bounded linear application from V to Λ h, we denote by K h h the approximation of K

K h h = v h ∈ V h , (B h h ) * v h , µ h = 0 ∀µ h ∈ Λ h .
The approximate saddle-point problem is

a(u h , v h ) + (B h h ) * λ h, v h = f, v h ∀v h ∈ V h µ h, B h h u h = 0 ∀µ h ∈ Λ h .
We have the following error estimate

u -u h V ≤ 1 + a α inf w h ∈K h h w h -u + 1 α inf µ h∈Λ h ξ -B h h * µ h V ′
where α is the coercivity constant of a and ξ is the linear form in V defined by a(u, v)

+ ξ, v = f, v ∀v ∈ V .
Proof. u h is the solution of the approximate saddle-point problem, so we have,

a(u h , v h ) = f, v h ∀v h ∈ K h h .
For all

w h ∈ K h h , we write v h = u h -w h ∈ K h h . Thus, a(v h , v h ) = a(u h -w h , v h ) = f, v h -a(w h , v h ) .
As u is the solution of our problem we have,

a(u, v h ) + ξ, v h = f, v h .
Thus, for any µ h in Λ h, one has,

a(v h , v h ) = a(u, v h ) + ξ, v h -a(w h , v h ) -(B h h v h , µ h) = a(u -w h , v h ) + ξ -B h h * µ h, v h .
Taking the absolute value of the previous equality gives us

α v h 2 ≤ a u -w h v h + ξ -B h h * µ h V ′ v h . But v h = u h -w h so we have α u h -w h ≤ a u -w h + ξ -B h h * µ h V ′ .
Thus,

u -u h ≤ u -w h + w h -u h ≤ 1 + a α w h -u + 1 α ξ -B h h * µ h V ′ ∀w h ∈ K h h , µ h ∈ Λ h ,
which concludes the proof.

Now we define the problem we are interested in. Let Ω be a bounded domain in R 2 and O a non empty open subset of Ω. We denote by P the following problem

(P) -∆u = f in Ω \ O u = 0 on γ u = 0 on ∂Ω ,
where f is a function of L 2 (Ω \ O) and γ denotes the boundary of O. We use a fictitious domain method and extend the function f by 0 in O. We still denote by f this extension and we denote by V the space H 1 0 (Ω) and by K the constrained space,

K = v ∈ V ; v |γ = 0 .
The Lagrange multipliers space is Λ = L 2 0 (γ), which is the space of L 2 functions on γ with zero mean value. Let B be the following application,

B : V -→ Λ v -→ v |γ .
We have K = ker B and the corresponding saddle-point problem is:

Find (u, λ) ∈ V × Λ, such that Ω ∇u • ∇v + γ λv = Ω f v ∀v ∈ V γ µu = 0 ∀µ ∈ Λ .
Let h > 0, we discretize the boundary γ by taking N h points x i, h ∈ γ so that the distance between two consecutive points is h. Let V h ⊂ V be a finite element space approximation of V . The Lagrange multipliers space and the constrained space are approximated as follows:

Λ h = R N h K h h = v h ∈ V h ; v h (x i, h) = 0 .
Here the space Λ h is not included in Λ. The application B becomes

B h h : V h -→ Λ h v h -→ v h (x i, h) i=1...N h .
We denote by (P h, h) the following approximated problem:

Find (u h , λ h) ∈ V h × Λ h, such that a(u h , v h ) + λ H , B h h v h = f, v h ∀v h ∈ V h µ h, B h h u h = 0 ∀µ h ∈ Λ h .
Where a is the bilinear form defined in V × V by

a(u, v) = Ω ∇u • ∇v .
We can apply the result of Proposition 3 so we have the following error estimate,

u -u h V ≤ 1 + a α inf w h ∈K h h w h -u + 1 α inf µ h∈Λ h ξ -B h h * µ h V ′ .
As in the previous example, the first part of this error estimate is the usual one and the second part is the error done by approximating the linear form ξ with a combination of Dirac masses. Indeed, for any µ h = (µ i ) 1,...,N h in Λ h and for any v h in V h , one has,

B h h * µ h , v h = µ h , B h h v h = N h i=1 µ i v h (x i ) = N h i=1 µ i δ x i , v h .
The error estimate of Proposition 3 becomes,

u -u h 1 ≤ 1 + a α inf w h ∈K h H w h -u + 1 α inf µ h∈Λ h ξ - N h i=1 µ i δ x i , v h v h .
Thus, by applying the result of Proposition 1,

u -u h 1 ≤ C   √ h u 1 + h h h1/2 ξ 0,γ   .
Again, if h is equal to h, the error estimate writes,

u -u h 1 ≤ C √ h ( u 1 + ξ 0,γ ) .

Numerical results

In this section we present two examples to validate the estimate (3) of Proposition 1. The first one is a one dimensional example to show that in the case h = h the error order is the one predicted by Proposition 1. The second example is a two dimensional numerical example where the solution of a Poisson equation with a right-hand side of the single layer distribution type is computed. This right-hand side is the Laplacian of a function with a jump of its normal derivatives across an interface. The exact solution is thus known and we compute the error in the case h = h for various values of s.

1d example

Let ϕ be the function defined by

ϕ = 1 x 1-s ,
where s is positive, and v h be the affine function defined in [0, h] by

v h (x) = - x -h 2 h + 1 .
The estimate of Proposition 1 can be computed explicitly. Indeed, one has,

ϕ, v h = h 0 ϕv h = hs s + 1 s(s + 1) - 1 s(s + 1)2 s - 1 2s hs+1 h .
Let x i be in [0, h], the second part is

ϕ h h , v h = h 0 ϕ v h (x i ) =   - x i -h 2 h + 1   hs s .
Let 0 ≤ α ≤ 1, we write x i = α h. Thus, the error estimate writes

ϕ, v h -ϕ h h , v h = C(s) hs+1 h ,
where C(s) is a constant depending only on s. Now, if h = C h, where C is a constant, the error is of order s which is the order predicted by Proposition 1. where f is given by the Laplacian of u outside the particle.

ρ γ Ω R R max χ = 1 χ = 0
To compute the numerical solution of this problem we truncate the series at some order N greater than 2π/h so that there is at least one period of the function sin(N θ) in a mesh cell. We compute the numerical solution for several s using Q 1 finite elements and compare it to the exact solution (see figure 3). In all the tests, the order N is equal to 2 12 and h takes its values from 2 -6 to 2 -10 . The space step h is equal to h so that the H 1 error should be of the form,

u -u h 1 ≤ C √ h|u| 2 + h s ϕ -1/2+s,γ .
We can see in figure 3 that the error order seems to stay close to 0.6 once s is greater than 0.5. This is because of the space approximation of order 1/2 since we use a cartesian mesh with finite elements of order 1, but it could also be because of the saturation of the error order when s is greater than 0.5 (see the remark right after Proposition 1).

Figure 4 shows the order of the error as s goes to 1. When s is less than 0.5, the slope is equal to 1 which is expected. The numerical order is always better than the theoretical one, but this kind of super-convergence often happens when one compute numerical solutions of PDE.

Figure 2 :

 2 Figure 2: Notations

Figure 3 :

 3 Figure 3: H 1 error for some value of s

Figure 4 :

 4 Figure 4: Order of the H 1 error with respect to s

Two-dimensional example

In this section Ω is the unit square and γ is a circle of radius R. To investigate the behavior of the approach for a right-hand side ϕδ γ , with ϕ in H -1/2+s (γ), we build the function ϕ as an infinite series of sine functions

We now build a function u such that the jump of its normal derivatives on γ is ϕ so that u is the solution of a Poisson problem of the form

where f is to be defined later on. First we need a cut function to make sure that the solution u is equal to zero on the boundary of Ω. Let us denote by χ ǫ the C 2 function defined in R by

We define ρ as the distance between γ and ∂Ω (see figure 2). We denote by R max the sum of R and ρ. Now, we can define the following function :

It is a function which has a jump of its normal derivative on γ and oscillates very fast close to γ. Its Laplacian outside the particle is given by ∆u

The jump ϕ of its normal derivatives is :

n -s sin(nθ)