Samir Adly
email: samir.adly@unilim.fr

Alberto Seeger
email: alberto.seeger@univ-avignon.fr

A nonsmooth algorithm for cone-constrained eigenvalue problems

Keywords: Complementarity problem, Cone-constrained eigenvalue problem, Semismooth Newton method, Polyhedral convex cone, Lorentz cone, Matrix

We study several variants of a nonsmooth Newton-type algorithm for solving an eigenvalue problem of the form

Such an eigenvalue problem arises in mechanics and in other areas of applied mathematics. The symbol K refers to a closed convex cone in the Euclidean space R n and (A, B) is a pair of possibly asymmetric matrices of order n. Special attention is paid to the case in which K is the nonnegative orthant of R n . The more general case of a possibly unpointed polyhedral convex cone is also discussed in detail.

1 K ∈ (R n ) ≡ set of closed convex cones in R n
is considered as general as possible, but one always assumes that x, Bx = 0 for all x ∈ K\{0}.

If one chooses B as the identity matrix, as we shall do in all our numerical experiments, then this assumption is automatically satisfied.

A fundamental problem arising in mechanics (cf. [START_REF] Figueiredo | A complementarity eigenproblem in the stability of finite-dimensional elastic systems with frictional contact[END_REF][START_REF] Figueiredo | The directional instability problem in systems with frictional contacts[END_REF][START_REF] Quittner | Spectral analysis of variational inequalities[END_REF]) and in other areas of applied mathematics is that of finding λ ∈ R and a nonzero vector x ∈ R n such that

K x ⊥ (Ax -λBx) ∈ K + (1)
with ⊥ indicating orthogonality and K + standing for the positive dual cone of K. This is what one calls a cone-constrained eigenvalue problem. One refers to x as a K-eigenvector of (A, B) and to λ as a K-eigenvalue of (A, B). The set of all K-eigenvalues, denoted by σ (A, B, K), is called the K-spectrum of (A, B). The mechanical or physical interpretation of a K-eigenvalue depends of course on the specific context. The orthogonality condition in [START_REF] Aubin | Set-Valued Analysis[END_REF] implies that λ = x, Ax x, Bx , as in a classical generalized eigenvalue problem, but Ax -λBx may be different from zero.

The mathematical theory of K-spectra is nowadays well developed [START_REF] Gowda | Some inertia theorems in Euclidean Jordan algebras[END_REF][START_REF] Lavilledieu | Existence de valeurs propres pour les systèmes multivoques: résultats anciens et nouveaux[END_REF][START_REF] Seeger | Eigenvalue analysis of equilibrium processes defined by linear complementarity conditions[END_REF][START_REF] Seeger | On eigenvalues induced by a cone constraint[END_REF][START_REF] Seeger | Local minima of quadratic forms on convex cones[END_REF], but the design of efficient algorithms is still in its infancy. In [START_REF] Humes | The symmetric eigenvalue complementarity problem[END_REF][START_REF] Júdice | On the solution of the symmetric eigenvalue complementarity problem by the spectral projected gradient algorithm[END_REF] it is assumed that both A and B are symmetric, but we are not interested in such a framework. We focus the attention to the more complicated case in which A is asymmetric. The Scaling-and-Projection Algorithm and the Power Iteration Method have been introduced and studied in [START_REF] Pinto Da Costa | Numerical resolution of cone-constrained eigenvalue problems[END_REF][START_REF] Pinto Da Costa | Cone-constrained eigenvalue problems: theory and algorithms[END_REF]. Global optimization and Branch-and-Bound techniques have been explored by Júdice et al. [START_REF] Júdice | The eigenvalue complementarity problem[END_REF]. In the present paper we discuss a solution technique that is very different in spirit from those just mentioned. The basic idea consists in reformulating (1) as a system (z) = 0 of nonsmooth equations and then applying a Newton-type algorithm. Such an approach makes sense only if the vector function : R d → R d satisfies a minimal set of assumptions. The details are explained next.

Brief description of the semismooth Newton method

To start with, one asks : R d → R d to be locally Lipschitz. Rademacher's theorem ensures then the existence of the Jacobian matrix J (z) for almost all z ∈ R d . Hence, at any reference point z, the Clarke generalized Jacobian ∂ (z) = co M ∈ M d : M = lim k→∞ J (z k ) for some {z k } k∈N → z with z k ∈ D is a nonempty compact convex set. Here, D stands for the set of differentiability points of and "co" refers to the convex hull operation.

The standard formulation of the Semismooth Newton Method (SNM) reads as follows:

• Initialization. Choose an initial point z 0 and set t = 0.

• Iteration. One has a current point z t . Choose M t ∈ ∂ (z t ) and compute h t by solving the linear system

M t h t = -(z t ). ( 2
)
Set z t+1 = z t + h t and increment t by one.

This algorithm has been studied in depth by many authors. Some additional assumptions on are needed for ensuring that (2) admits a unique solution and that {z t } t∈N converges. The theorem below is taken from Qi and Sun [START_REF] Qi | A nonsmooth version of Newton's method[END_REF]Sect. 3], but other results in the same vein can be found in [START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF][START_REF] Klatte | Nonsmooth Equations in Optimization[END_REF][START_REF] Pang | Newton's method for B-differentiable equations[END_REF].

Theorem 1 Let z be a zero of the locally Lipschitz function : R d → R d . Suppose that is semismooth at z, and

(3)

all matrices in ∂ (z) are nonsingular. ( 4
)
Then there exists a neighborhood V of z such that the SNM initialized at any z 0 ∈ V generates a sequence {z t } t∈N that converges superlinearly to z.

We opted for a formulation of Theorem 1 that is as simple as possible, but the specialized literature in the area provides also conditions guaranteeing a quadratic rate of convergence of the sequence {z t } t∈N . The technical assumption (3) refers to the existence of the limit

(z; h) := lim M∈∂ (z+tw) w→h, t→0 + Mw ( 5
)
for all direction h ∈ R d . This is equivalent to saying that sup

M∈∂ (z+v) (z + v) -(z) -Mv = o( v ) as v → 0.
Despite the unpleasant appearance of the expression [START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF], there are plenty of helpful rules for checking semismoothness. While applying the SNM to cone-constrained eigenvalue problems, the main source of annoyances is the nonsingularity hypothesis (4). We shall comment on this issue in due course.

Eigenvalue problems under nonnegativity constraints

Theoretical background

The most interesting example of a cone-constrained eigenvalue problem is the one in which the convex cone K is the nonnegative orthant of R n . In such a case, (1) takes the form

x ≥ 0, Ax -λBx ≥ 0, x, Ax -λBx = 0, (6)
where the notation x ≥ 0 indicates that each component of x ∈ R n is nonnegative. This is the so-called Pareto eigenvalue problem. One refers to x as a Pareto eigenvector of (A, B) and to λ as a Pareto eigenvalue of (A, B). The model ( 6) is not as sui generis as it might appear. In fact, a cone-constrained eigenvalue problem associated to a pointed polyhedral convex cone can be brought to the form [START_REF] Fan | A minimax inequality and applications[END_REF]. The details are explained in the next proposition borrowed from [START_REF] Seeger | Eigenvalue analysis of equilibrium processes defined by linear complementarity conditions[END_REF]. Recall that a convex cone is pointed if it contains no line.

Proposition 2 Let K = {Gu : u ∈ R p + }
for some real matrix G whose columns are positively linearly independent vectors of R n . Then

σ (A, B, K) = σ G T AG, G T BG, R p + , ( 7
)
where the superscript "T" indicates transposition.

The integer p corresponds to the number of generators of K. A pointed polyhedral convex cone in R n can always be represented as in Proposition 2, but sometimes the columns of G are not easily available. Needless to say, the evaluation of the righthand side of (7) could be very expensive if p is large.

As a first step towards a reformulation of the Pareto eigenvalue problem (6) as a system of equations, we write x ≥ 0 primal or eigenvector feasibility, y ≥ 0 dual feasibility,

x, y = 0 complementarity slackness, [START_REF] Figueiredo | The directional instability problem in systems with frictional contacts[END_REF] Ax -λBx = y stationarity or equilibrium law,

1 n , x = 1 normalization
with 1 n ∈ R n standing for a vector of one's. The last condition in (8) makes sure that x is a nonzero vector. Normalizing x with the help of the linear function 1 n , • is just one option among many. One could use instead the normalization condition x 2 = 1 and this would not change substantially the overall presentation of our work. In the sequel, we refer to y ∈ R n as a "dual vector" (or as a vector of dual variables). We are borrowing the terminology of optimization theory, but one should not view (8) as a Karush-Kuhn-Tucker system for an optimization problem. Recall that A and B are not necessarily symmetric.

Next, we observe that the first three conditions in (8) can be written in the succinct form U ϕ (x, y) = 0, where U ϕ : R 2n → R n is the vector function given by

U ϕ (x, y) = ⎡ ⎢ ⎣ ϕ(x 1 , y 1 ) . . . ϕ(x n , y n ) ⎤ ⎥ ⎦
and ϕ : R 2 → R is a complementarity function, i.e.,

ϕ(a, b) = 0 ⇐⇒ a ≥ 0, b ≥ 0, ab = 0. (9)
In short, we are led to solve a system

U ϕ (x, y) = 0, ( 10
)
Ax -λBx -y = 0, ( 11
)
1 n , x -1 = 0 ( 12
)
of 2n + 1 equations involving the same number of variables. The troublesome part of the above nonlinear system is the block [START_REF] Galántai | ABS methods for nonlinear systems of algebraic equations[END_REF], and this is because a complementarity function is usually nonsmooth. The next lemma explains why it is reasonable to apply the SNM to the system (10)- [START_REF] Gowda | Some inertia theorems in Euclidean Jordan algebras[END_REF]. We mention in passing that concept of semismoothness applies to any locally Lipschitz function between finite dimensional vector spaces, even if the output space has different dimension than the input space. The same remark applies to the definition of a Clarke generalized Jacobian. Lemma 3 Consider a vector function : R 2n+1 → R 2n+1 of the form

(x, y, λ) = ⎡ ⎣ U ϕ (x, y) Ax -λBx -y 1 n , x -1 ⎤ ⎦ ( 13
)
with ϕ being any of the following complementarity functions:

ϕ min (a, b) = min{a, b}, ( 14
)
ϕ FB (a, b) = a + b -a 2 + b 2 , ( 15
)
ϕ EP (a, b) = 2ab -(min{0, a + b}) 2 . ( 16
)
Then is locally Lipschitz and semismooth. Furthermore, its Clarke generalized Jacobian at z = (x, y, λ) admits the block representation

∂ (z) = ⎧ ⎨ ⎩ ⎡ ⎣ E F 0 A -λB -I n -Bx 1 T n 0 0 ⎤ ⎦ : [E, F ] ∈ ∂U ϕ (x, y) ⎫ ⎬ ⎭ ( 17
)
with I n standing for the identity matrix of order n.

Proof The last n + 1 components of are continuously differentiable, so they are locally Lipschitz and semismooth. Concerning the first n components of , we recall the following known facts:

-The "lattice" complementarity function ( 14) is a minimum of two linear functions, therefore it is globally Lipschitz and semismooth. -The Fischer-Burmeister complementarity function [START_REF] Júdice | On the solution of the symmetric eigenvalue complementarity problem by the spectral projected gradient algorithm[END_REF] is globally Lipschitz. Furthermore, it is continuously differentiable everywhere except at the origin, where it is semismooth (cf. [START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF]Sect. 7.4]). -The Evtushenko-Purtov complementarity function ( 16) is continuously differentiable because the nonsmooth term min{0, a + b} has been squared. That [START_REF] Kanzow | A new class of semismooth Newton-type methods for nonlinear complementarity problems[END_REF] satisfies the "only if" part of ( 9) is mentioned for the first time in [START_REF] Evtushenko | Sufficient conditions for a minimum for nonlinear programming problems[END_REF]. The "if" part can be seen by direct inspection. Concerning the name attributed to ( 16), we follow references [START_REF] Galántai | ABS methods for nonlinear systems of algebraic equations[END_REF][START_REF] Tseng | Growth behavior of a class of merit functions for the nonlinear complementarity problem[END_REF].

This takes care of the announced behavior of and the representation formula [START_REF] Klatte | Nonsmooth Equations in Optimization[END_REF].

For any of the complementarity functions listed in Lemma 3, the computation of the Clarke generalized Jacobian ∂U ϕ (x, y) offers no difficulty. For saving space, we omit writing down the explicit formulas. Many other complementarity functions could also be considered (cf. [START_REF] Chen | A family of NCP functions and a descent method for the nonlinear complementarity problem[END_REF][START_REF] Kanzow | A new class of semismooth Newton-type methods for nonlinear complementarity problems[END_REF]), but there is no need of multiplying the examples. For notational convenience, we write

U min = U ϕ min , U FB = U ϕ FB , U EP = U ϕ EP ,
and define the functions min , FB , EP according to [START_REF] Humes | The symmetric eigenvalue complementarity problem[END_REF].

Of course, the SNM specialized to the system (10)-( 12) takes a particular form. It reads as follows:

• Initialization. Choose an initial point z 0 = (x 0 , y 0 , λ 0 ) and set t = 0. • Iteration. 1 One has a current point z t = (x t , y t , λ t ). Choose (E t , F t ) ∈ ∂U ϕ (x t , y t ) and compute z t+1 = (x t+1 , y t+1 , λ t+1 ) by solving the linear system

E t x t+1 + F t y t+1 = E t x t + F t y t -U ϕ (x t , y t ), (A -λ t B)x t+1 -y t+1 -λ t+1 Bx t = -λ t Bx t , 1 n , x t+1 = 1.
There are no theoretical results on the best way of selecting the initial point. Guided just by common sense, we apply a rule of the sort:

generate a random vector ξ ∈ R n with uniform distribution on [-1, 1] n , ( 18
)
and set

x 0 = ξ/ 1 n , ξ , λ 0 = x 0 , Ax 0 / x 0 , Bx 0 , y 0 = Ax 0 -λ 0 Bx 0 . ( 19
)
Generating an n-dimensional Gaussian vector ξ is another option, but we stick to the uniform distribution rule. A question raised by one of the referees is why we do not generate ξ uniformly on [0, 1] n . After all, the sequence {x t } t∈N is expected to converge to a vector with nonnegative components. This observation is very pertinent, but it happens in practice that an overwhelming majority of Pareto eigenvectors lie on the boundary of the cone R n + . Now, a boundary type Pareto eigenvector can equally well be approximated from the interior or from the exterior of the cone. Anyhow, we carried out some numerical experiments with the referee's choice and we did not observe a substantial difference with respect to [START_REF] Kong | Fischer-Burmeister complementarity function on Euclidean Jordan algebras[END_REF].

An important question now is this: if z = ( x, ȳ, λ) is a solution to ( 10)-( 12), are we sure that every matrix in ∂ (z) is nonsingular? The answer to this question is not always yes. It depends not only on the matrices A and B, but also on the solution z under consideration.

Example 4 Consider the Pareto eigenvalue problem for the pair

A = 3 -1 4 -1 , B= 1 0 0 1
and the lattice complementarity function [START_REF] Júdice | The eigenvalue complementarity problem[END_REF]. Although min is differentiable at both solutions

z bad = ((1/3, 2/3), (0, 0), 1)
,

z good = ((1, 0), (0, 4), 3), the corresponding Jacobian matrices ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 0 0 1 0 0 0 0 0 1 0 2 -1 -1 0 -1/3 4 -2 0 -1 -2/3 1 1 0 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ , ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 0 0 1 0 0 0 1 0 0 0 0 -1 -1 0 -1 4 -4 0 -1 0 1 1 0 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦
are somewhat different in nature. The former is singular and the latter is nonsingular.

The situation described by Example 4 warns us from being overly optimistic, but it does not impede us from using the SNM. What happens in practice is that some solutions could remain undetected. From a theoretical point of view, it is interesting to distinguish between different types of solutions to the Pareto eigenvalue problem [START_REF] Fan | A minimax inequality and applications[END_REF].

Definition 5 Let x be a Pareto eigenvector of (A, B) and write

λ = x, A x / x, B x , ȳ = A x -λ x, z = ( x, ȳ, λ).
The Pareto eigenvector x is called:

(i) Strict if xi + ȳi > 0 for all i ∈ {1, . . . , n}.
(ii) Nondegenerate relative to ϕ if the function ( 13) is such that every matrix in ∂ (z) is nonsingular.

Strictness simply means that the variables xi and ȳi cannot vanish at the same time. Such a requirement has a clear interpretation in mechanics and, more generally, in any mathematical model involving complementarity slackness. The analytical benefit derived from strictness is differentiability: Proposition 6 Under the notation of Definition 5, the following statements are equivalent:

(a) The Pareto eigenvector x is strict. (b) U min is continuously differentiable on a neighborhood of ( x, ȳ). (c) U FB is continuously differentiable on a neighborhood of ( x, ȳ).

Proof For a Pareto eigenvector x, strictness can be expressed in the equivalent form xi = ȳi for all i ∈ {1, . . . , n}.

The proposition is a consequence of this observation and the following two facts: firstly, ϕ min is continuously differentiable in a neighborhood of ( ā, b) if and only if ā = b, and, secondly, ϕ FB is continuously differentiable everywhere except at the origin.

Unfortunately, strictness does not prevent degeneracy. The risk of degeneracy persists even in a smooth context and simply cannot be avoided. Anyway, it is good to have the following result in mind.

Proposition 7 Let x be a strict Pareto eigenvector of (A, B). Then the following statements are equivalent:

(a) x is nondegenerate with respect to ϕ min . (b) x is nondegenerate with respect to ϕ FB . (c) x is nondegenerate with respect to ϕ EP .

Proof The complementarity functions ϕ min and ϕ FB are calibrated so that

∂ϕ ∂a ( ā, 0) = 0, ∂ϕ ∂b ( ā, 0) = 1, for all ā > 0, ( 20
)
∂ϕ ∂a (0, b) = 1, ∂ϕ ∂b (0, b) = 0, for all b > 0. ( 21
)
Hence, U min and U FB have the same Jacobian matrix at the reference point ( x, ȳ). Consequently, min and FB have the same Jacobian matrix at z. This is a lot more than merely saying that (a) and (b) are equivalent. The partial derivatives of ϕ EP are not exactly as in ( 20)-( 21), but they satisfy a calibration condition that is somewhat similar. Indeed, one has

∂ϕ EP ∂a ( ā, 0) = 0, ∂ϕ EP ∂b ( ā, 0) = 0, for all ā > 0, ∂ϕ EP ∂a (0, b) = 0, ∂ϕ EP ∂b (0, b) = 0, for all b > 0.
The Jacobian matrices J min (z) and J EP (z) are different in general. However, we claim that det J EP (z) = μ det J min (z) [START_REF] Pinto Da Costa | Numerical resolution of cone-constrained eigenvalue problems[END_REF] for some scalar μ = 0. To see this, suppose for instance that xj > 0 for j = 1, . . . , r, xj = 0 for j = r + 1, . . . , n with r ∈ {1, . . . , n}. Any other sign configuration of x can be treated in the same way.

The strictness assumption implies that ȳj = 0 for j = 1, . . . , r, ȳj > 0 for j = r + 1, . . . , n.

Hence,

J min (z) = ⎡ ⎣ E F 0 A -λB -I n -B x 1 T n 0 0 ⎤ ⎦ (23)
with

E = 0 0 0 I n-r , F = I r 0 0 0 .
The Jacobian matrix J EB (z) has the same structure as [START_REF] Pinto Da Costa | Cone-constrained eigenvalue problems: theory and algorithms[END_REF], except that E and F are slightly different: I r must be changed by a diagonal matrix with nonzero entries

∂ϕ EP ∂b ( x1 , 0), . . . , ∂ϕ EP ∂b ( xr , 0), ( 24
)
and I n-r must be changed by a diagonal matrix with nonzero entries

∂ϕ EP ∂a (0, ȳr+1 ), . . . , ∂ϕ EP ∂a (0, ȳn ). ( 25
)
Hence, the equality [START_REF] Pinto Da Costa | Numerical resolution of cone-constrained eigenvalue problems[END_REF] holds if one defines μ as the product of all the terms in ( 24)-( 25). This proves that (a) and (c) are equivalent.

Numerical experiments

In all our tests2 , B = I n and A is an arbitrary matrix, not necessarily symmetric. By a Pareto eigenvalue of A one understands a scalar λ ∈ R such that the complementarity system

x ≥ 0, Ax -λx ≥ 0, x, Ax -λx = 0 ( 26
)
admits a nonzero solution x ∈ R n . Such an x is called a Pareto eigenvector of A.

It is worthwhile mentioning that the integer

π n = max A∈M n card[σ (A, I n , R n + )]
is finite, but grows exponentially with n. By definition, π n is the largest number of Pareto eigenvalues that can be found in a matrix of order n. As shown in [START_REF] Pinto Da Costa | Cone-constrained eigenvalue problems: theory and algorithms[END_REF], one has

2 n -1 ≤ π n ≤ n2 n-1 -(n -1).
Note, for instance, that a matrix of order 40 could have more than 10 12 Pareto eigenvalues! Finding all of them is simply not realistic, so we aim at detecting just a bunch of them.

Testing on a matrix of order 3 with 9 Pareto eigenvalues

Our first numerical experiment concerns a matrix of order 3 that is known to have 9 Pareto eigenvalues:

A = ⎡ ⎣ 8 -1 4 3 4 1/2 2 -1/2 6 ⎤ ⎦ . ( 27
)
As seen in the Table 1, the SNM was able to detect all of them. 3 We considered 10 4 random initial points, but we could have used a smaller sample. What Table 1 essentially says is that some Pareto eigenvalues are more likely to be found than others. For easy of visualization, proportions presented as percentages are given without decimals. As one can see, the percentage of convergence towards each particular solution depends on the choice of complementarity function.

⎤ ⎥ ⎥ ⎦ . ( 28
)
We tested the SNM on (28) by using the complementarity function ϕ FB and a sample of 10 3 random initial points. Table 2 displays what we got. The matrix (28) possesses exactly 23 Pareto eigenvalues, so the SNM found all of them. The Pareto eigenvectors are all strict in the sense of Definition 5. Some faces of the cone R 4 + contain several Pareto eigenvectors, but other faces are free of them.

Testing on random matrices of various dimensions

Table 3 has been produced as follows. For each n ∈ {10, 20, . . . , 50} ∪ {100, 200}, one generates a sample of 10 4 random pairs (A, ξ ) with ξ as in [START_REF] Kong | Fischer-Burmeister complementarity function on Euclidean Jordan algebras[END_REF] and A following a uniform distribution on the unit hypercube [0, 1] n×n . One applies the SNM to (A, I n ) with initial point z 0 constructed from ξ as explained in [START_REF] Lavilledieu | Existence de valeurs propres pour les systèmes multivoques: résultats anciens et nouveaux[END_REF]. One counts the number of times in which convergence to a solution occurs within 100 iterations. More than 100 iterations, or the occurrence of a badly conditioned Jacobian matrix at current point, is declared a failure. Some of the conclusions that can be drawn from Table 3 are:

-The overall performance of the SNM diminishes as the dimension n increases. This is something to be expected from any numerical method. We saw such a deterioration phenomenon already with the Scaling-and-Projection Algorithm [START_REF] Pinto Da Costa | Cone-constrained eigenvalue problems: theory and algorithms[END_REF]. -The highest percentages of convergence are obtained with the complementarity function ϕ FB . The performance of ϕ min is almost as good as ϕ FB . Despite being continuously differentiable, ϕ EB yields very disappointing percentages of convergence. Said in a provocative way, a nonsmooth complementarity function may perfectly well beat a smooth one. -Concerning the average number of iterations, the quality of ϕ FB deteriorates as n increases, but the quality of ϕ min remains rather stable.

Remark 8 Similar numerical tests with A following a uniform distribution on [-1, 1] n×n lead essentially to the same conclusions. However, the percentages of convergence are weaker, regardless of the choice of complementarity function.

3 Partially constrained eigenvalue problems

Coping with unpointedness

In some applications only a group of components of x ∈ R n is cone-constrained. To fix the ideas, suppose that the first m components of x are required to be nonnegative, but the remaining ones are free. For notational convenience we write x in the form

x = x c x f ,
where the constrained block x c is an m-dimensional vector and the free block x f is an (nm)-dimensional vector. Consistently, we write

y = y c y f , A= A c,c A c,f A f,c A f,f , B= B c,c B c,f B f,c B f,f .
We are dealing thus with a generalized eigenvalue problem constrained by a convex cone

K m,n-m = R m + × R n-m ( 29
)
that is no longer pointed. Hence, we go beyond the framework of Sect. 2 and of reference [START_REF] Pinto Da Costa | Numerical resolution of cone-constrained eigenvalue problems[END_REF]. Recall that the lineality of a convex cone K in R n is defined as the integer

lin(K) = dim[K ∩ -K].
The lineality of ( 29) is clearly nm. The dual cone K + m,n-m = R m + × {0} n-m has empty interior and the spectral problem (1) becomes

x c ≥ 0,
x f free,

y c ≥ 0, y f = 0, x c , y c + x f , y f = 0, A c,c A c,f A f,c A f,f x c x f -λ B c,c B c,f B f,c B f,f x c x f = y c y f .
After adding the usual normalization condition and simplifying, one ends up with

U ϕ (x c , y c ) = 0, (30)
A c,c x c + A c,f x f -λ(B c,c x c + B c,f x f ) -y c = 0, (31)
A f,c x c + A f,f x f -λ(B f,c x c + B f,f x f ) = 0, (32)
1 m , x c + 1 n-m , x f -1 = 0. (33)
Note that ( 30)-( 33) is a nonsmooth system of n + m + 1 equations involving the same number of variables. Table 4 has been produced as follows. For each dimension n ∈ {20, 40, 60, 80} one generates a sample of 10 4 random pairs (A, ξ ) with ξ as in [START_REF] Kong | Fischer-Burmeister complementarity function on Euclidean Jordan algebras[END_REF] and A following a uniform distribution on the hypercube [0, 1] n×n . One applies the SNM to (A, I n ) with an initial point z 0 constructed from ξ as explained in [START_REF] Lavilledieu | Existence de valeurs propres pour les systèmes multivoques: résultats anciens et nouveaux[END_REF]. One counts the number of times in which convergence to a solution occurs within 100 iterations. The same procedure is carried out for each lineality nm ranging from 1 to 5. We use ϕ FB as complementarity function. The main conclusion of Table 4 is that, for each dimension n, the SNM performs slightly better as the number m of constrained variables increases.

The general polyhedral case

Dealing with a general polyhedral convex cone K in R n is not more difficult than dealing with the particular case [START_REF] Stern | Exponential nonnegativity on the ice cream cone[END_REF]. The first thing one has to do is writing

K = {Gu + F v : u ∈ R p + , v ∈ R q }, (34)
where G and F are real matrices whose columns are vectors in R n satisfying a suitable blend of linear independence and positive linear independence:

Gu + F v = 0, u∈ R p + , v ∈ R q =⇒ u = 0, v = 0. ( 35
)
Finding the matrices G and F is often a cumbersome task, but there are a number of interesting examples for which these matrices are readily available. The result stated below extends Proposition 2.

Proposition 9 Let K be as in (34) with G, F satisfying the hybrid linear independence condition (35). Then σ (A, B, K) = σ ( Ã, B, K p,q ), where

à = G T AG G T AF F T AG F T AF , B = G T BG G T BF F T BG F T BF
are matrices of order p + q.

Proof The dual cone of (34) is clearly

K + = {y ∈ R n : G T y ≥ 0, F T y = 0}. Plugging x = Gu + F v into (1) leads to the system u ≥ 0, v free, G T [A(Gu + F v) -λB(Gu + F v)] ≥ 0, F T [A(Gu + F v) -λB(Gu + F v)] = 0, Gu + F v, A(Gu + F v) -λB(Gu + F v) = 0,
which can be written in the more compact form

K p,q u v ⊥ G T AG G T AF F T AG F T AF -λ G T BG G T BF F T BG F T BF u v ∈ K + p,q . ( 36
)
Given the hybrid linear independence condition (35), x is a nonzero solution to (1) if and only if (u, v) is a nonzero solution to (36).

Example 10 By way of illustration, consider a generalized eigenvalue problem constrained by the downward monotonic cone

K down = {x ∈ R n : x 1 ≥ • • • ≥ x n }.
For this particular case, one can take

G = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 1 1 . . . 1 0 1 . . . 1 . . . . . . . . . . . . 0 0 . . . 1 0 0 . . . 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , F = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 1 . . . . . . . . . 1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ .
Note that G is formed by p = n -1 columns and F is formed by q = 1 column. This is consistent with the fact that K down ∩ -K down = R1 n is a one-dimensional linear subspace. The case of an upward monotonic cone can be handled in a similar way.

Extensions

Beyond a polyhedral setting

The Lorentz eigenvalue problem is the prototype example of a cone-constrained eigenvalue problem involving a nonpolyhedral convex cone. It can be formulated in the succinct form

0 x ⊥ (Ax -λBx) 0, ( 37
)
where stands for the Lorentz ordering relation on R n . By definition, u v means that uv belongs to the Lorentz cone

L n = x ∈ R n : [x 2 1 + • • • + x 2 n-1 ] 1/2 ≤ x n .
Some authors refer to L n as the n-dimensional ice-cream cone [START_REF] Loewy | Positive operators on the n-dimensional ice cream cone[END_REF][START_REF] Stern | Exponential nonnegativity on the ice cream cone[END_REF] or as the second order cone [START_REF] Boyd | Applications of second-order cone programming[END_REF]. A nonzero vector x satisfying (37) is called a Lorentz eigenvector of the pair (A, B). The reader should be aware that the Lorentz spectrum σ (A, B, L n ) may contain infinitely many elements, even if one assumes that

x, Bx = 0 for all x ∈ L n \{0}.

The lack of polyhedrality in L n is at the origin of many complications. Anyhow, the SNM can be applied also to (37). Indeed, the Lorentzian eigenvalue problem can be reformulated as a nonsmooth system of equations, namely

U (x, y) = 0, (38)
Ax -λBxy = 0, (39)

x n -1 = 0. ( 40
)
The equilibrium law (39) is as in Sect. 2.1. The normalization condition (40) guarantees that x is a nonzero vector. This is the simplest and most natural way of normalizing a Lorentz eigenvector. The function U : R n × R n → R n appearing in (38) is a vector complementarity function for the Lorentz ordering. This simply means that U (x, y) = 0 ⇐⇒ x 0, y 0, x ⊥ y.

For making sure that the SNM works properly, one asks U to be locally Lipschitz and semismooth. A vector complementarity function U with such properties does exist, the two most popular examples being

U min (x, y) = x -L n (x -y), ( 41
)
U FB (x, y) = x + y - √ x • x + y • y . ( 42
)
The symbol L n in (41) stands for the metric projection onto L n . The product • and the square root in (42) are to be understood in the Jordan algebra framework associated to L n . The details concerning the definition and properties of U min and U FB can be found in [START_REF] Fukushima | Smoothing functions for second-order-cone complementarity problems[END_REF][START_REF] Gowda | Some P-properties for linear transformations on Euclidean Jordan algebras[END_REF][START_REF] Kong | Fischer-Burmeister complementarity function on Euclidean Jordan algebras[END_REF].

Beyond affine matrix pencils

The expression A -λB is a particular instance of a matrix pencil. So, the coneconstrained eigenvalue problem (1) is a particular case of the model

K x ⊥ M(λ)x ∈ K + (43)
with M standing for the pencil associated to a finite collection {A 0 , A 1 , . . . , A r } of real matrices of order n, that is,

M(λ) = r k=0 λ k A k . ( 44
)
The way we handle an eigenvalue problem like (43) is as follows. To fix the ideas, consider for instance the Paretian version

x ≥ 0, M(λ)x ≥ 0, x, M(λ)x = 0. (45)
It should be clear that solving (45) is a matter of finding a zero of

(x, y, λ) = U ϕ (x, y) M(λ)x -y 1 n , x -1 (46)
with ϕ being any of the complementarity functions listed in Lemma 3. The function (46) is locally Lipschitz and semismooth. Furthermore, its Clarke generalized Jacobian at z = (x, y, λ) is given by

∂ (z) = ⎧ ⎨ ⎩ ⎡ ⎣ E F 0 M(λ) -I n M (λ)x 1 T n 0 0 ⎤ ⎦ : [E, F ] ∈ ∂U ϕ (x, y) ⎫ ⎬ ⎭ with M (λ) = r k=1 kλ k-1 A k .
In principle, everything is more or less the same as in Sect. 2.1. Of course, the expression "dual vector" refers now to y = M(λ)x. There is however an important theoretical difference between an affine pencil like A -λB and a polynomial pencil like (44). In the latter case, the existence of Pareto eigenvalues is not automatically guaranteed. We end this work with an existential result for the case of a quadratic pencil

M(λ) = A 0 + λA 1 + λ 2 A 2 .
(47)

Cone-constrained eigenvalue problems involving quadratic pencils arise, for instance, in the analysis of mechanical systems with frictional contacts (cf. [START_REF] Figueiredo | A complementarity eigenproblem in the stability of finite-dimensional elastic systems with frictional contact[END_REF]). An application to electronic circuits with ideal diodes will be discussed in a forthcoming paper of ours.

Theorem 11 Let M be the quadratic pencil associated to a triplet (A 0 , A 1 , A 2 ) of real matrices of order n such that

x, A 2 x = 0 for all x ∈ R n + \{0}, x, A 1 x 2 -4 x, A 2 x x, A 0 x ≥ 0 for all x ∈ R n + .
Then there exist a scalar λ ∈ R and a nonzero vector x ∈ R n solving (45).

Proof We rely on a theorem by Ky Fan on the existence of solutions to variational inequalities (cf. [1, Theorem 3.1.1] or the original source [START_REF] Fan | A minimax inequality and applications[END_REF]). Let = {x ∈ R n + : 1 n , x = 1}. Note that is a convex compact set not containing the origin. Both roots

λ ± (x) = -x, A 1 x ± [ x, A 1 x 2 -4 x, A 2 x x, A 0 x ] 1/2 2 x, A 2 x
of the quadratic equation x, M(λ)x = 0 are real and well defined for every x ∈ . Consider for instance the function : × → R given by (x, w) = M(λ + (x))x, w .

One readily sees that ⎧ ⎨ ⎩ for all w ∈ , ( • , w) is continuous, for all x ∈ , (x, • ) is linear, for all x ∈ , (x, x) = 0. Ky Fan's theorem ensures the existence of x ∈ such that ( x, w) ≥ 0 for all w ∈ . If one sets λ = λ + ( x), then it is clear that ( λ, x) solves (45) with M(λ) as in (47).

Example 12

The quadratic pencil satisfies the hypotheses of Theorem 11. In fact, (48) satisfies even stronger assumptions, namely, A 2 is positive definite and

M(λ) = λ 2 ⎡ ⎣ 2
x, A 1 x 2 -4 x, A 2 x x, A 0 x > 0 for all x ∈ R 3 \{0}.

The existence of solutions to (45) is then guaranteed. Table 5 displays the solutions we got by implementing the SNM with the complementarity function ϕ FB and a sample of 10 4 random initial points. We used the initialization rule

x 0 = ξ/ 1 n , ξ , λ 0 = λ ± (x 0 ), y 0 = M(λ 0 )x 0 with ξ being a random vector as in [START_REF] Kong | Fischer-Burmeister complementarity function on Euclidean Jordan algebras[END_REF].

Note that λ 6 and λ 7 are eigenvalues of (48) in the classical sense. Indeed, the corresponding eigenvectors are in the interior of R 3 + and therefore the corresponding dual vectors are equal to zero. All the other Pareto eigenvalues are not classical, but created by the cone constraint imposed on the pencil.

By way of conclusion

The SNM is well suited for solving the cone-constrained eigenvalue problem (1), regardless of whether the matrices A and B are symmetric or not. The SNM is also successful when it comes to solve the more general problem (43) involving a polynomial matrix pencil.

More important than the structure of A and B, what really makes a practical difference is the complexity of the cone K. If this cone is polyhedral, then (1) can be converted into a system of equations like in ( 30)-(33), and the SNM becomes directly applicable. As reasonable to expect, the nonpolyhedral case is more involved. One has to worry about finding a suitable vector complementarity function for the cone. Such function must be semismooth and, at the same time, it must have an easily computable Clarke generalized Jacobian.

It is hard to compare the SNM with other methods for solving [START_REF] Aubin | Set-Valued Analysis[END_REF]. For instance, the numerical results reported in [START_REF] Júdice | The eigenvalue complementarity problem[END_REF] are different in spirit from ours and focus on different issues. What is safe to say, at least, is this:

-While working with some particular examples, like the one of Sect. 2.2.1, the SNM was able to find the solutions that were left undetected by the Scaling-and-Projection Algorithm [START_REF] Pinto Da Costa | Cone-constrained eigenvalue problems: theory and algorithms[END_REF] and by the Power Iteration Method [START_REF] Pinto Da Costa | Numerical resolution of cone-constrained eigenvalue problems[END_REF][START_REF] Pinto Da Costa | Cone-constrained eigenvalue problems: theory and algorithms[END_REF].

This does not mean that the SMN is systematically better than the methods suggested in [START_REF] Pinto Da Costa | Numerical resolution of cone-constrained eigenvalue problems[END_REF][START_REF] Pinto Da Costa | Cone-constrained eigenvalue problems: theory and algorithms[END_REF]. As said before, implementing the SNM could be a headache if the cone K is not polyhedral. By contrast, implementing the Scaling-and-Projection Algorithm offers no difficulty if one has to deal with a revolution cone of arbitrary half-aperture angle [START_REF] Pinto Da Costa | Numerical resolution of cone-constrained eigenvalue problems[END_REF].

Table 1

 1 Testing on a matrix of order 4 with 23 Pareto eigenvaluesAn example of a matrix of order 4 with 23 Pareto eigenvalues was kindly communicated to us by Victor Qi, a student at Tongji University, Shanghai. Qi's example is this:

	SNM applied to the matrix (27) with a sample of 10 4	Solution	% of convergence to each solution
	random initial points. The		ϕ min	ϕ FB	ϕ EP
	figures on the last three columns
	refer to percentages of convergence to each particular solution. The cases of	λ 1 = 4.1340 λ 2 = 4.6021	24 13	31 7	22 14
	divergence are not counted, so	λ 3 = 5.0000	8	6	7
	the percentages in each column add up to 100%	λ 4 = 5.8660 λ 5 = 6.0000	13 13	10 13	8 10
		λ 6 = 7.0000	9	5	6
		λ 7 = 8.0000	15	15	13
		λ 8 = 9.3979	2	6	10
		λ 9 = 10.0000	3	7	10

Table 2

 2 SNM applied to the matrix[START_REF] Seeger | Local minima of quadratic forms on convex cones[END_REF] with a sample of 10 3 random initial points. Pareto eigenvectors are found in various faces of the Pareto cone R 4

				+
	Pareto	Pareto eigenvector			Dual vector
	eigenvalue	x 1	x 2	x 3	x 4	y 1	y 2	y 3	y 4
	λ 1 = 26.2823	0.4314 0.0762 0	0.4924	0	0	0.7693	0
	λ 2 = 26.4149	0.4558 0.0368 0.0581 0.4493	0	0	0	0
	λ 3 = 28.7114	0.4527 0	0.1913 0.3559	0	1.1100	0	0
	λ 4 = 29.1341	0.2266 0.2491 0	0.5243	0	0	7.7457	0
	λ 5 = 32.6080	0	0.4461 0	0.5539	2.4261	0	15.7526	0
	λ 6 = 32.8635	0.4258 0	0.2844 0.2897	0	3.0863	0	0
	λ 7 = 37.5767	0.2238 0	0.7762 0	0	1.9626	0	4.7001
	λ 8 = 41.0162	0.1241 0.0681 0.8078 0	0	0	0	5.1944
	λ 9 = 46.4681	0	0.1771 0.8229 0	3.9579	0	0	6.7290
	λ 10 = 49.1435	0.1561 0.1589 0.4874 0.1976	0	0	0	0
	λ 11 = 66.9700	0	0.3429 0.4566 0.2005	11.8334	0	0	0
	λ 12 = 77.4251	0.7814 0	0.0010 0.2176	0	49.8944	0	0
	λ 13 = 77.4575	0.7823 0	0	0.2177	0	49.9815	0.0406	0
	λ 14 = 99.4233	0.9690 0	0.0310 0	0	88.3988	0	20.3481
	λ 15 = 100.0000 1	0	0	0	0	9 2 .0000	2.0000 21.0000
	λ 16 = 107.5010 0	0.5132 0.3019 0.1849	33.9922	0	0	0
	λ 17 = 127.3920 0	0.7674 0	0.2326	62.5095	0	32.1390	0
	λ 18 = 148.5319 0	0.7171 0.2829 0	70.9204	0	0	27.2498
	λ 19 = 158.0000 0	1	0	0	1 0 6 .0000	0	44.0000 38.0000
	λ 20 = 197.1730 0.3415 0.4238 0.1155 0.1193	0	0	0	0
	λ 21 = 204.5836 0.3874 0.4820 0	0.1306	0	0	21.0694	0
	λ 22 = 226.2813 0.3935 0.4888 0.1178 0	0	0	0	26.8356
	λ 23 = 231.9223 0.4455 0.5545 0	0	0	0	25.2880 30.4261

Table 4

 4

	SNM applied to (A, I n ) with A following a	Dimension	Number m of constrained variables
	uniform distribution on	n	m = n/4	m = n/2	m = 3n/4
	[0, 1] n×n . The percentages of
	convergence are given without decimals and are estimated by using a sample of 10 4 randomly	20 40	91% 93%	95% 96%	98% 98%
	generated pairs (A, ξ )	60	92%	95%	97%
		80	91%	94%	96%

Table 5

 5 SNM applied to the quadratic pencil (48) with a sample of 10 4 random initial points. Pareto eigenvectors were found in various faces of R 3

			+
	Pareto	Pareto eigenvector		Dual vector
	eigenvalue	x 1	x 2	x 3	y 1	y 2	y 3
	λ 1 = -4.3930	0	1	0	6.0000	0	5.0000
	λ 2 = -3.7656	1	0	0	0	2.0000	0
	λ 3 = -3.6524	0.8712	0.1288	0	0	0	0.6438
	λ 4 = -2.0000	0	0	1	0	3.0000	0
	λ 5 = -1.9613	0	0.1318	0.8682	0.7909	0	0
	λ 6 = -1.9580	0.0954	0.1277	0.7769	0	0	0
	λ 7 = -0.7689	0.3877	0.4006	0.2116	0	0	0
	λ 8 = -0.6986	0.5036	0.4964	0	0	0	2.4820
	λ 9 = -0.6820	0	0.6426	0.3574	3.8554	0	0
	λ 10 = -0.6070	0	1	0	6.0000	0	5.0000
	λ 11 = 0.0000	0	0	1	0	3.0000	0
	λ 12 = 0.2656	1	0	0	0	2.0000	0

As stopping criterion we use (z t ) < 10 -8 .

The numerical experiments are carried out in a Powerbook Mac OS 10.5.8 with a processor 2.33 GHz Intel Core 2 Duo and 2Go memory. The software used is Matlab 7.7.0.471 (R2008b).

In[START_REF] Pinto Da Costa | Numerical resolution of cone-constrained eigenvalue problems[END_REF], the Power Iteration Method implemented with 10 8 random initial points detected only four Pareto eigenvalues of[START_REF] Seeger | On eigenvalues induced by a cone constraint[END_REF].

Acknowledgements

The authors are grateful to Professor C. Kanzow (Würzburg) for useful hints concerning the semismooth Newton method.


G Chardon

Optimal design of sampling sets for least-squares signal recovery with the Frank-Wolfe algorithm

Keywords: Least-squares, sensor selection, convex optimization, optimal design

We consider the sensors selection problem in a leastsquares setting. The sensors selection is replaced by the relaxed problem of designing a sampling density minimizing the number of samples needed to ensure stability of the recovery, shown to be equivalent to the D-optimal design problem. We propose to use the Frank-Wolfe algorithm to solve this optimization problem, with low space and time computational complexity, linear with respect to the number of possible sensors positions. As the optimal densities are usually sparse, sampling points are drawn from the optimized density using resampling methods. The optimization problem and procedure can be easily modified to account for additional design constraints.

I. INTRODUCTION

The problem of sampling signals is of particular importance in digital signal processing, as one necessary step linking the real, analog world, to the discrete computations of digital devices.

In the classical sampling theorem, the compact support of the Fourier transform of a signal is leveraged to allow its perfect reconstruction from its values on a discrete set of finite density. Starting from low-pass signal models and uniform sampling, this line of work continues through sampling of signal with bandpass spectrum [START_REF] Vaughan | The theory of bandpass sampling[END_REF] or arbitrary known Fourier support in higher dimensions [START_REF] Landau | Necessary density conditions for sampling and interpolation of certain entire functions[END_REF], to sampling of signal with unknown Fourier support of finite measure. In this case, vertices of quasi-crystals are easily built sets that are known to be stable sampling sets [START_REF] Matei | Quasicrystals are sets of stable sampling[END_REF]. Compressed sensing and its applications put once again the question of sampling set design on the foreground of signal processing and harmonic analysis, by considering the recovery of vectors that are sparse in a given basis with a number of samples lower than the dimension of the basis. [START_REF] Foucart | A Mathematical Introduction to Compressive Sensing, ser. Applied and Numerical Harmonic Analysis[END_REF]. The design of compressed sensing measurement schemes is difficult, and usually relies on randomness, either by using random matrices [START_REF] Candes | Near-Optimal Signal Recovery From Random Projections: Universal Encoding Strategies?[END_REF], or by choosing random lines from a deterministic matrix [START_REF] Candes | Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information[END_REF].

In this article, we consider the least-squares recovery of a signal f defined on a set X. We use a simple linear model where f is approximated by a signal f m in a space V m of finite dimension m spanned by m basis functions (φ i ) 1≤i≤m . An estimate fm of f is found by matching fm to noisy samples of f in a least-squares sense. A suboptimal sampling set can lead to instability of the recovery, i.e. an estimation error fm -f 2 2 Gilles Chardon is with the Laboratoire des Signaux et Systèmes (L2S), CentraleSupélec-CNRS-Université Paris-Sud, Université Paris-Saclay, 3 rue Joliot-Curie, 91192 Gif-sur-Yvette, France. email: gilles.chardon@centralesupelec.fr. larger than the approximation error f m -f 2 2 or than the noise level by several orders of magnitude.

These instabilities can be reduced by regularizing the estimation of fm . In this approach, priors on the noise and on the signal f are needed, which are not necessarily available. Moreover, regularization reduces the variance of the estimation, but also introduces bias. Our approach will be to design a sampling set such that no regularization is required, i.e. that no prior is necessary apart the knowledge that f can be approximated by f m ∈ V m with sufficient precision.

The problem of sensors placement is an important problem in signal processing [START_REF] Joshi | Sensor Selection via Convex Optimization[END_REF], [START_REF] Ranieri | Near-Optimal Sensor Placement for Linear Inverse Problems[END_REF], [START_REF] Chepuri | Sparsity-Promoting Sensor Selection for Non-Linear Measurement Models[END_REF], but such considerations are no limited to this field. Similar problems are found in numerical analysis, in particular when solving PDEs using the reduced basis method [START_REF] Maday | A general multipurpose interpolation procedure: the magic points[END_REF], or in optimal design of experiments in statistics [START_REF] Pukelsheim | Optimal Design of Experiments[END_REF] (in particular, it will be shown that the problem of D-optimal design is equivalent to the problem of stable least-squares recovery).

A particular difficulty of the design of sampling sets is the non-convexity of the optimization problem. Methods such as convex relaxation of the problem as the optimization of a probability density [START_REF] Joshi | Sensor Selection via Convex Optimization[END_REF], [START_REF] Castro | Approximate Optimal Designs for Multivariate Polynomial Regression[END_REF], [START_REF] Chepuri | Sparsity-Promoting Sensor Selection for Non-Linear Measurement Models[END_REF], as well as greedy algorithms [START_REF] Ranieri | Near-Optimal Sensor Placement for Linear Inverse Problems[END_REF], [START_REF] Wynn | The Sequential Generation of D-Optimum Experimental Designs[END_REF], [START_REF] Shamaiah | Greedy sensor selection: Leveraging submodularity[END_REF], can be used to construct approximate solutions to the optimal sampling set problem. We will here use the former approach, also known as variable density compressive sampling in sparse recovery [START_REF] Puy | On Variable Density Compressive Sampling[END_REF], [START_REF] Chauffert | Variable Density Sampling with Continuous Trajectories[END_REF].

Our main goal is the minimization of the number of samples necessary to ensure stability of the reconstruction. To this end, we use a lower bound on the necessary number of samples needed for stability in function of the probability density used to draw the sampling points, introduced in [START_REF] Cohen | On the Stability and Accuracy of Least Squares Approximations[END_REF]. We propose a fast algorithm for the design of the measurement density to minimize this bound, as well as a sampling scheme to construct the actual measurement set from the optimized density.

A. Contributions

The problem of designing a sampling set for stable leastsquares estimation is shown to be equivalent to the D-optimal design problem, in the sense that optimizing a sampling density to minimize the necessary number a samples for stable recovery is equivalent to maximizing the determinant of the information matrix of a linear regression problem. Furthermore, in generic settings, the optimal densities are discrete.

The application of the Frank-Wolfe algorithm to this problem was introduced in [START_REF] Chardon | Design of variable densities for least-squares approximations[END_REF]. It is here shown that the complexity in time and space of this algorithm is lower than existing methods, with similar or better results. In particular, the complexity is linear with respect to the number of possible sensor positions. This algorithm, yielding sparse densities, is well adapted to the properties of optimal sampling densities.

The algorithm can take in account additional constraints or modified objective functions. As an example, the joint optimization of a common sampling density for two different spaces V m is demonstrated. The number of samples necessary for stable approximation in both spaces is also given by the optimization algorithm.

Finally, as the optimal sampling densities are sparse, the use of resampling methods is discussed. These methods yield sampling sets that ensure stability with a number of measurements close to the dimension of the space V m , with a minimal amount of randomness. The code necessary to reproduce the figures is available online 1 .

B. Structure of the article

Results on the stability of least-squares recovery are recalled in section II, as well as connections with D-optimal design and other methods. The optimization problem and algorithm are formulated in section III, with numerical comparisons with existing algorithms. The problem of sampling the optimized density and numerical results on estimation errors are discussed in section IV. Extensions of the method are introduced in section V, and concluding remarks are given in section VI.

II. STABLE LEAST-SQUARES RECOVERY AND RELATED

PROBLEMS

The standard setting for the least-squares recovery of a signal is as follows. We assume that the signal f to be recovered can be approximated by f m in a space of finite dimension V m , and that we observe samples y i = f (x i ) + b i , where x i ∈ X are n fixed sampling points, and the measurement errors b i are i.i.d. centered random variables of variance σ 2 . The space V m depends on the application. Examples include standard or trigonometric polynomials [START_REF] Cohen | On the Stability and Accuracy of Least Squares Approximations[END_REF], Fourier-Bessel functions or spherical harmonics [START_REF] Chardon | Sampling and Reconstruction of Solutions to the Helmholtz Equation[END_REF], [START_REF] Chardon | Design of Spatial Microphone Arrays for Sound Field Interpolation[END_REF], cosine basis [START_REF] Ranieri | Near-Optimal Sensor Placement for Linear Inverse Problems[END_REF], etc.

The least-squares estimate fm of f is given by

fm = argmin fm∈Vm n i=1 |y i -fm (x i )| 2 . ( 1
)
The Gauss-Markov theorem shows that in the case where f = f m , fm is the best linear unbiased estimate of f . Our goal is to ensure the stability of the estimation with respect to the measurement noise, as well as the approximation error, i.e. that the estimation error can be bounded:

fm -f 2 2 ≤ C 1 f m -f 2 2 + C 2 σ 2 (2)
with C 1 as close as possible to 1, and C 2 as close as possible to 0 for a given number of samples. Depending on the choice of the sampling points x i , a large number of samples can be necessary to ensure stability of the estimation, even in noiseless cases. A famous example is 1 https://gilleschardon.fr/fwkm known as the Runge phenomenon, where polynomial interpolation of a function with regularly spaced samples fails for a number of samples linear with respect to the dimension m.

Optimizing the sampling points x i leads to an untractable non-convex optimization problem. We will use a relaxed approach, introduced in [START_REF] Cohen | On the Stability and Accuracy of Least Squares Approximations[END_REF], where the sampling points are drawn from a probability measure µ. The problem of choosing points x i is replaced by choosing a density µ. A similar approach is adopted for sparse recovery, where results similar to (2) are obtained through the use of random matrices.

A. Stability of least-squares estimation

In this setting, Cohen et al. [START_REF] Cohen | On the Stability and Accuracy of Least Squares Approximations[END_REF] give a criterion for the stability of least-squares recovery. With (L j ) 1≤j≤m an orthogonal basis of V m with respect to the measure µ, they define

K(x, µ) = m j=1 |L j (x)| 2 (3)
K(µ) = sup x∈X K(x, µ). (4)
K(µ) is bounded from below by m. The following theorem shows that the expectation of the estimation error is bounded in a way similar to eq. ( 2).

Theorem 1. For a given r > 0, κ = (1 -log 2)/(2 + 2r), and

ε(n) = 4κ/ log n, if the number of measurements n is such that K(µ) ≤ κ n log n (5)
then the expectation of the estimation error is bounded by

E fm -f µ ≤ (1+2ε(n)) f m -f 2 µ +8M 2 n -r +8σ 2 m n ( 6
)
where M is an upper bound on |f | and

f 2 µ = X |f | 2 dµ (7)
K(µ) essentially measures the number of samples necessary to ensure stability with respect to deterministic approximation errors and random measurement noise. Before introducing the optimization algorithm to design a sampling density µ such that K(µ) is close to m, we recall alternative state of the art methods that will be tested against our proposed method.

B. D-Optimal design

The criterion K(µ) also appears in the analysis of Doptimal design [START_REF] Pukelsheim | Optimal Design of Experiments[END_REF]. In this setting, the quantity of interest are the parameters α i of the decomposition of f in the basis

(φ i ) 1≤i≤m of V m . f (x) = m i=1 α i φ i (x). (8)
Here, deterministic approximation errors are not considered (i.e., f = f m ).

The parameters α i are estimated using measurements at sampling points x i ∈ X, that can be repeated w i times with independent noise realizations. The problem of D-optimal design is the maximization of the determinant

D = det M (9)
where M is the information matrix with coefficients

M ij = 1 n n l=1 w n φ i (x l )φ j (x l
). This determinant is proportional the inverse of the volume of the uncertainty ellipsoid for the estimation of the parameters α i . This problem can also be relaxed, by defining M (µ) with coefficients M (µ) ij = X φ i (x)φ j (x)dµ. In this approximate design problem, it has been shown that maximizing D(µ) is equivalent to minimizing K(µ) [START_REF] Kiefer | The equivalence of two extremum problems[END_REF], in the sense that

argmin µ K(µ) = argmax µ D(µ). (10)
Furthermore, the following bound [START_REF] Kiefer | Optimum Experimental Designs V, with Applications to Systematic and Rotatable Designs[END_REF]

D(µ) ≥ exp(m -K(µ)) sup µ D(µ) (11)
shows that the maximization of the determinant can be replaced by the minimization of K(µ), of which the convergence is easier to control. Indeed, it will be shown that the lower bound m can be reached. Moreover, the variance of the estimation of f (x) is σ 2 K(x, µ). Different algorithms have been proposed for the design of sampling sets in the D-optimal sense:

• greedy algorithms [START_REF] Wynn | The Sequential Generation of D-Optimum Experimental Designs[END_REF], [START_REF] Shamaiah | Greedy sensor selection: Leveraging submodularity[END_REF],

• multiplicative algorithms [START_REF] Silvey | An algorithm for optimal designs on a design space[END_REF], and combination with greedy algorithms [START_REF] Yu | D-optimal designs via a cocktail algorithm[END_REF],

• convex optimization methods, either on a relaxed problem [START_REF] Joshi | Sensor Selection via Convex Optimization[END_REF], or on a moment problem in the particular case of polynomial approximation [START_REF] Castro | Approximate Optimal Designs for Multivariate Polynomial Regression[END_REF],... In [START_REF] Joshi | Sensor Selection via Convex Optimization[END_REF], instead of optimizing a probability density, Joshi and Boyd constrain the sampling density to sum to the number of desired sensors n. The optimization problem is solved by Newton's method, and sampling points are selected as the n largest values of the density. A local optimization can then refine the choice of the sampling set.

C. FrameSense

The FrameSense algorithm [START_REF] Ranieri | Near-Optimal Sensor Placement for Linear Inverse Problems[END_REF], proposed by Ranieri et al., is based on a proxy for the mean squared estimation error when no deterministic measurement errors are present. For a given set of points x i , the rows of the matrix with coefficients φ ij = φ i (x j ) define a frame of R m . In the case where this frame has elements with unit norm, it is shown that minimizing the estimation error is equivalent to minimizing the frame potential

i,j | l φ il φ jl |
2 . An argument of sub-modularity is invoked to justify the use of a greedy algorithm to minimize the frame potential.

D. Weighted least-squares

An alternative to the standard least-squares is to use weights in eq. ( 1) :

fm = argmin fm∈Vm n i=1 w(x i )|y i -fm (x i )| 2 . ( 12
)
A joint optimization procedure for the sampling density µ and weights w(x) based on a modified version of Theorem 1 is proposed in [START_REF] Cohen | Optimal weighted least-squares methods[END_REF]. This procedure is not well adapted to the case of measurement design, as imposing weights implies that not only the location of the measurements, but also the variance of the noise applied to each measurement can be independently controlled.

III. THE MINIMIZATION PROBLEM

We aim at designing sampling sets ensuring stability of the estimation with the smallest number of samples. To this end, supported by Theorem 1, we consider the following optimization problem:

µ = argmin µ K(µ) (13)
under the constraint that µ is a probability measure. As our objective is the numerical optimization of a sampling set, we will limit ourselves to a discretized version of the problem. The space V m is a subspace of a larger space E of dimension L. V m is described by a matrix P containing the basis vectors. The sampling density µ is replaced by a positive vector w of weights w i , and the basis Q is an orthogonal basis of V m with respect to the weights w (Q diag(w)Q = I). Its row vectors are denoted q i (i.e. each q i contains the i-th coefficient of the basis vectors of V m ).

In this discrete setting, the criterion K is defined by

K(i, w) = q i 2 2 (14) K(w) = max i K(i, w) (15)
and the optimization problem becomes:

w = argmin w K(w) subject to w i ≥ 0, L i=1 w i = 1 (16)
Theorem 2. The optimization problem ( 16) is convex and at the optimal density w

K(w ) = m, (17)
and

w i = 0 if K(i, w ) = m. ( 18
)
A proof is given in [START_REF] Kiefer | The equivalence of two extremum problems[END_REF] in the continuous case, where the support of µ is included in the set {x|K(x, µ ) = m}. An alternative proof is given in appendix A using standard convex analysis tools in the discrete case. An important consequence of the support condition is that in most cases, the optimal density is sparse. Indeed, the functions used to build the space V m are usually analytic (exponentials, polynomials, cosines, etc.), and so is K(x, µ ). Except in the unlikely cases where K(x, µ ) is constant (e.g. for a Fourier basis), K(x, µ ) = m almost everywhere, implying that the support of the measurement density has Lebesgue measure 0.

The optimization problem [START_REF] Chauffert | Variable Density Sampling with Continuous Trajectories[END_REF], despite its apparent complexity (cost of evaluation of the objective function which involves the orthogonalization of a large matrix, large number of constraints), can be solved in an efficient way using the Frank-Wolfe algorithm. As will be shown, this algorithm allows to leverage the particular form of the objective function and of the constraints to efficiently compute a minimizer of [START_REF] Chauffert | Variable Density Sampling with Continuous Trajectories[END_REF].

A. Frank-Wolfe algorithm

Input: function f , feasible set D, initialization x 0 , number of iterations N Output: last iterate x N n = 0 while n < N do γ = 2 n+2 s = argmin s∈D s, grad f (x n ) x n+1 = (1 -γ)x n + γs n = n + 1 end Algorithm 1: Frank-Wolfe algorithm
The Frank-Wolfe algorithm, also known as the conditional gradient method, is a convex optimization algorithm introduced in 1956 by Frank and Wolfe [START_REF] Frank | An algorithm for quadratic programming[END_REF], [START_REF] Jaggi | Revisiting Frank-Wolfe: Projection-Free Sparse Convex Optimization[END_REF].

It is similar to a gradient descent, with the difference that the next iterate is not chosen in the direction of the gradient, but in the direction of the minimizer in the feasible set of the linear approximation of the objective function around the current iterate. The domain being convex, the iterates are guaranteed to be feasible, and no projection in the feasible set is necessary. The algorithm is outlined in Alg. 1.

In our case K(w) is the maximum of the smooth functions K(i, w). The Franke-Wolfe algorithm can be extended to such non-smooth functions [START_REF] White | Extension of the Frank-Wolfe algorithm to concave nondifferentiable objective functions[END_REF]. The gradient of K(w) is the convex hull of the gradients of the maximal K(i, w), and at most points, the gradient of K(w) is simply the gradient of the maximum K(i, w), with index i .

The linear problem to solve at each iteration is:

min s s, gradK(i , w) subject to L i=1 s i = 1, s i ≥ 0. ( 19
)
Its solution is a vector with zero coefficients, except at the minimal value of the gradient, where it is one. A straightforward application of the Cauchy-Schwartz inequality to eq. ( 36) shows that the minimal coordinate of the gradient of K(i , w) is the i -th coordinate.

An iteration of the Frank-Wolfe algorithm applied to our optimization problem has thus the simple implementation:

• find the index i of the maximum value of K(i, w), • add weight on the i -th coefficient of the sampling density, • rescale the sampling density such that L i=1 w i = 1. As at most one point is added to the support at each iterations, the Frank-Wolfe algorithm yields sparse measurement densities.

The results of the Frank-Wolfe algorithm for K(w) minimization, Newton's method for the equivalent problem of determinant maximization [START_REF] Joshi | Sensor Selection via Convex Optimization[END_REF], and the measurement density for weighted least-squares [START_REF] Cohen | Optimal weighted least-squares methods[END_REF] are plotted on figure 1, for a polynomial basis of size m = 20. For Newton's method, the solution of the relaxed problem for n = 2m sensors is plotted. As expected, the obtained density is sparse: K(i, w ) being a polynomial of order 2m -2, the support of w , where K(i, w ) = m, contains at most 2m -2 points. The Frank-Wolfe optimization procedure for K(w) minimization can also be used to solve the approximate Doptimality problem. As shown above, the determinant maximization problem and the K(µ) minimization problem have the same solutions. Moreover, applying the Frank-Wolfe algorithm to the determinant maximization problem gives the same optimization procedure. Indeed, the Fisher information matrix of the approximate D-optimal design problem (which is also the Gram matrix of P with respect to the weights w) is given by

M = P diag(w)P (20)
and the determinant lemma shows that the derivative with respect to the weight w i of the determinant D = det M is

∂D ∂w i = p i M -1 p i D = q i q i D = K(i, w)D (21)
where the orthogonal basis Q is defined as Q = PM -1/2 and p i are the row vectors of P. The maximal value of the gradient is at index i , identical to the index selected by eq. ( 19).

B. Time and space complexity

The computational cost of an iteration of the Frank-Wolfe algorithm is dominated by the orthogonalization of the design matrix P of size L × m. While a standard orthogonalization algorithm would have a complexity in O(Lm 2 ), the orthogonalization can be accelerated by using the particular way the sampling density is updated. Indeed, an iteration of the algorithm modifies the density at a unique point, and renormalizes the density. As this modification amounts to a rank-1 perturbation of the Gramian matrix, its squared root can be computed using the matrix inversion lemma with complexity O(Lm), which governs the cost of an iteration.

While no theoretical bound on the number of iterations necessary to achieve a value of K(w) close to the optimum

Time (s)

Frank-Wolfe Weighted LS Newton FrameSense Fig. 3. Computational time of the proposed Frank-Wolfe algorithm, Newton's method [START_REF] Joshi | Sensor Selection via Convex Optimization[END_REF], FrameSense [START_REF] Ranieri | Near-Optimal Sensor Placement for Linear Inverse Problems[END_REF] and weighted least-squares [START_REF] Cohen | Optimal weighted least-squares methods[END_REF], in function of the size of the discretization L, with m = 20, for gaussian matrices. is yet available, numerical experiments show that in various settings, a number of iterations proportional to m is sufficient to reach values of K(w) close to the optimum m. Fig. 2 shows on the same plot the value of K(w) along the iterations, for various dimensions L and m, and various bases (polynomials and random matrices with Rademacher, gaussian, and uniformly distributed independent coefficients). In all settings, the number of iterations necessary to reach a density w such that (K(w) -m) ≤ Cm is approximately m/C. For a given C, the computation complexity of the algorithm is thus O(Lm 2 ).

The time and space complexity of the proposed algorithm, Newton's method [START_REF] Joshi | Sensor Selection via Convex Optimization[END_REF], the FrameSense algorithm [START_REF] Ranieri | Near-Optimal Sensor Placement for Linear Inverse Problems[END_REF] and the joint design of a sampling density and weights for weighted least-squares [START_REF] Cohen | Optimal weighted least-squares methods[END_REF] are given in table I. The proposed algorithm, having time and space complexities linear in the number of possible sensor locations, is much more efficient for increasing L and fixed m, and can attain higher dimensions because of its low memory requirements, as can be seen on figure 3 (the computation times are given for a laptop with an Intel Core i7 quad core at 2.10 GHz and 16 GB of RAM,

) O(Lm) O(m) O(Lm 2 ) O(Lm) Weigted LS [25] O(Lm 2 ) - - O(Lm 2 ) O(Lm) FrameSense [8] O(L 2 m) O(L) L -n O(L 2 m) O(L 2 ) Newton [7] O(1) O(L 3 ) ≤ 30 O(L 3 ) O(L 2 )
running MATLAB R2017b). The complexity of the joint weights and density design for weighted least-squares is of the same order as the Frank-Wolfe algorithm, but is faster as it is equivalent to its initialization. While the method implies control of the measurement noise variance and does not yield sparse densities, its reduced numerical cost makes it a valuable alternative.

IV. SAMPLING THE OPTIMAL DENSITY

Once a density is designed, the measurement points are obtained by sampling this density. After introducing various sampling methods, numerical results on estimation errors are given.

A. Sampling methods

The sampling methods we consider are the following:

a) Maximal values of the density : Joshi and Boyd [START_REF] Joshi | Sensor Selection via Convex Optimization[END_REF] proposed to select the locations of the n largest values of the weights w i . While this method works in random settings, it does not produce efficient sampling sets in more typical cases. Indeed, as is visible on figure 1, the n largest values of the optimized density are not necessarily located at all sensor positions necessary for a stable recovery. They also propose a local optimization to improve the sampling set, by testing swaps of a sampling point with a non-selected point. b) i.i.d. samples: Theorem 1, based on a i.i.d. sampling of the density, is valid for a number of measurements such that

K(µ) ≤ κ n log n (22)
This quantity of samples ensures that the empirical Gram matrix of the basis (L j ), with coefficients

G e ij = 1 n n l=1 L i (x l )L j (x l ) (23)
is close enough to the coefficients of the Gram matrix with respect to the sampling measure µ

G ij = X L i (x)L j (x)dµ (24) = δ ij . ( 25
)
In addition to the large number of samples compared to K(m) and the dimension of the space m, using i.i.d. sampling in the cases where the sampling density is concentrated on few points leads to multiple selections of a given point, which is usually to be avoided in practice.

c) Resampling: The coefficients of the empirical Gram matrix in eq. ( 23) are a Monte-Carlo evaluation of the integrals [START_REF] Yu | D-optimal designs via a cocktail algorithm[END_REF]. Moreover, in most cases, the optimal density is discrete. Sampling of discrete densities is involved in the resampling stage of particle filters, or sequential Monte-Carlo methods. In this setting, the performances of an i.i.d. sampling of a discrete density (also called multinomial sampling in this case) can be improved by considering other resampling methods.

Stratified sampling and systematic sampling are well known simple resampling methods in dimension 1 [START_REF] Douc | Comparison of resampling schemes for particle filtering[END_REF]. With F (x) the cumulative distribution function of the sampling density, the samples are chosen to be

x i = F -1 i -1 n + u i ( 26
)
for stratified sampling, and

x i = F -1 i -1 n + u ( 27
)
for systematic sampling, where u, u 1 , . . . , u n are i.i.d uniform random variables in the interval [0, 1/n].

To avoid multiple selections of a sampling point, duplicates can be replaced by samples drawn from the uniform density. In higher dimensions, methods based on parametrizing the domain using an Hilbert curve and stratified or systematic sampling along this curve can be used [START_REF] Schretter | A direct inversion method for non-uniform quasi-random point sequences[END_REF], [START_REF] Gerber | Negative association, ordering and convergence of resampling methods[END_REF].

d) FrameSense as a sampling method: The above sampling methods are based on the sampling density only, and do not take into account the row matrix q i associated to a point x i . In the case when two rows are correlated, selecting both the rows does not add much information. Another less probable but uncorrelated row could yield more information. Removing similar rows is the mode of operation of FrameSense, which can be used as a sampling method after optimization of the density using the proposed method. Here only the support of the optimized density is used.

Conversely, the Frank-Wolfe algorithm can be considered as a pre-processing for FrameSense, selecting important points and lowering the computational complexity. The total cost of the proposed method followed by FrameSense is O(Lm 2 ), to be compared with the cost O(L 2 m) for FrameSense alone.

B. Numerical results

1) Polynomials:

(t) = (1 + 25t 2 ) -1 in the interval [-1, 1].
The mean squared error is estimated by repeating the experiment 10000 times, for 20 to 40 measurements, and the interval is discretized using 2000 points. The FrameSense algorithm, as well as systematic sampling combined with Newton's method algorithm, the proposed Frank-Wolfe algorithm and weighted least-squares, have performances similar to the optimal density (Dirac masses of weight 1/m at the extremities of the interval and at zeros of the Legendre polynomial of order m -2 [START_REF] Pukelsheim | Optimal Design of Experiments[END_REF]). Designing the sampling set using stratified sampling of the optimized density reaches the same performances only when n ≈ 2m (the increase of the error between n = 31 and n = 32 is a consequence of the non identical distributions for each sampling point). The performances of a set of regularly spaced points in the interval is also plotted for completeness.

The mean squared errors of sets of independent samples of the density optimized by the Frank-Wolfe algorithm, as well as the n largest values of the result of Newton's method and the locally optimized sampling set, are too large to be included in the figure .  2) Fourier-Bessel functions: Acoustical fields can be estimated using the approximation of solutions to the Helmholtz equation ∆p + k 2 p = 0 by so-called generalized harmonic polynomials (in R 2 ) or spherical harmonics (in R 3 ) [START_REF] Chardon | Design of Spatial Microphone Arrays for Sound Field Interpolation[END_REF], [START_REF] Moiola | Plane wave approximation of homogeneous Helmholtz solutions[END_REF]. In polar coordinates:

p(r, θ) ≈ m n=-m p n J n (kr)e inθ (28)
where J n is the Bessel function of order n, and m = 2m + 1.

In the case of a circular domain, it was proven that values of K(µ) close to the lower bound are obtained with densities that are concentrated on the boundary of the domain [START_REF] Chardon | Sampling and Reconstruction of Solutions to the Helmholtz Equation[END_REF]. Numerical evidence showed that this was also the case for square domains. Figure 5 shows the optimized density w in the case of a polygonal domain with a reentrant corner. Here m = 10, m = 2m + 1 = 21, and the square is discretized using L = 100 points on each axis, i.e. L = 7500 points in the domain. As expected, most of the optimized density (66%) is concentrated on the boundary of the domain.

The mean square error for the approximation of a plane wave on this domain is plotted on figure 6 for sampling sets obtained by the FrameSense algorithm, as well as set sampled from densities optimized with the Newton and the Frank-Wolfe algorithms. The convex relaxation approach (Frank-Wolfe or Newton) or weighted least-squares with Hilbert parametrization and systematic sampling yields accurate reconstructions. Using FrameSense as a postprocessing of the Frank-Wolfe algorithm improves the results. FrameSense alone is less accurate than the relaxed approached, and selecting the maximal values of the density optimized with Newton's method does not yield accurate estimations, but is somewhat improved by the local optimization. The computation times are 0.12s for the proposed Frank-Wolfe method, 1.8s for the FrameSense algorithm, and one minute for Newton's method (with L = 250, L = 46875, 2 seconds are necessary for the Frank-Wolfe algorithm, while the FrameSense and Newton's method fails because of too large memory requirements).

V. EXTENSIONS

Formulating the sensors selection problem as a convex optimization problem allows to modify the constraints and/or the objective functions to take into account more complex settings. We demonstrate here two simple modifications.

A. Bounding the error

A consequence of the concentration of the optimal sampling density µ on a set of measure 0 in important settings (e.g. polynomial bases, cosine bases, etc.) is that the bound ( 6) is unable to control the L 2 -norm of the error with respect to the Lebesgue measure λ.

To ensure stability in the L 2 -norm, a bound can be applied to the measure, e.g. µ(I) ≥ αλ(I) for any subset I of the domain and α > 0. The L 2 -norm of the error with respect to µ can now be used to bound the error with respect to the Lebesgue measure : A bound on the optimal measure with this additional constraint can be found. The measure µ α = (1 -α)µ + αλ being feasible, a simple computation (see ( 35)) shows that

f λ ≤ 1 α f µ . ( 29
)
K(µ α ) ≤ K((1 -α)µ + αλ) ≤ 1 1 -α K(µ ) (30)
In the discrete formulation, the optimization problem becomes ŵ = argmin w K(w) subject to

L i=1 w i = 1, w i ≥ α/L. ( 31
)
Algorithm 1 is easily modified to solve this problem. The value of K(w) for a polynomial basis and a gaussian matrix with L = 500, m = 20 is plotted for α between 0 and 1. The bound [START_REF] Schretter | A direct inversion method for non-uniform quasi-random point sequences[END_REF] appears to be pessimistic in these cases. In the polynomial case in particular, a bound of the error over the entire interval can be obtained at a minimal cost as K(w) remains almost equal to m for α < 0.6.

B. Joint optimization for multiple bases

In some applications, a unique sensor array has to be used to sample several signals that are not necessarily described by a unique linear model. An example is soundfield measurement (or more generally, wavefield estimation), where the same microphone array is used at different frequencies. The approched used in [START_REF] Chepuri | Sparsity-Promoting Sensor Selection for Non-Linear Measurement Models[END_REF] for sensor placement in non-linear models involves a similar problem.

In the case of S different subspaces, not necessarily of the same dimension, a quantity K s (w) can be computed for each subspace. We define

K(w) = max s K s (w), (32)
and the problem [START_REF] Chauffert | Variable Density Sampling with Continuous Trajectories[END_REF] with this modified objective function can be solved with a simple application of Alg. 1. This method is tested for the approximation of a function using polynomials (space V p m , with criterion K p ) and a cosine basis (space V d m , with criterion K d ). The following densities are computed:

• w p , optimized for the polynomial basis, • w d , optimized for the cosine basis,

• the mixture w m = (w p + w d )/2 • and w j , jointly optimized for the two bases. In this case, a simple analysis of the optimization problem shows that K p (w j ) = K d (w j ). The values of K p and K d for the mixed and joint densities are plotted on figure 8 in function of m. While the mixed density yields values of K p (w m ) close to the upper bound 2m, for the jointly optimized density, K d and K p can be bounded by 1.3m, making it possible to construct a sampling set capable of stable recovery in two different spaces with a reasonable number of measurements.

Estimation errors for the different combinations of spaces and sampling densities are plotted on figure 9, using systematic sampling. Here m = 20 and K p (w j ) = K d (w j ) = 24.1. As expected, 24 samples are sufficient to ensure the stability of the estimation simultaneously in the polynomial and cosine bases.

VI. CONCLUSION

The problem of sampling set design for least-squares estimation is here considered in a relaxed setting. The Frank-Wolfe algorithm offers an efficient optimization procedure in terms of computational complexity in time and space. Variants of this algorithm, such as the away step variant [START_REF] Guélat | Some comments on Wolfe's 'away step[END_REF] could improve the speed of convergence of the optimization.

The optimal densities, as well as the Frank-Wolfe iterates, being sparse in generic cases, the sampling points are obtained by resampling techniques. This approach is shown to be more efficient than i.i.d. sampling or selecting the maximal values of the density with local optimization. While efficient in dimension 1 (they allow stable estimation with a number of measurements close to the number of parameters to estimate), the resampling methods considered here do not seem optimal in higher dimensions. The development of sampling methods aware of the rows q i associated to a point w i along to its probability should improve the performances of the sampling sets.

Finally, the combination of the convex optimization formulation with resampling methods offers a wide arrays of possibilities, such as the optimization of sampling densities for non-linear problems, or the optimization of E-criterion in a way similar to [START_REF] Chepuri | Sparsity-Promoting Sensor Selection for Non-Linear Measurement Models[END_REF].

APPENDIX A MINIMUM OF K(w)

In this appendix, we prove the convexity of K(w) and that the optimal value of K(w) is the dimension m of the approximation space V m .

We consider a discrete setting. A basis of the space V m is given by the matrix P of dimension L × m, with row vectors p i . We call G w = P diag(w)P the Gram matrix of the basis P with respect to the set of weights w.

A. Convexity of K(w)

As K(w) = max i K(i, w), it is sufficient to prove the convexity of K(i, w). We consider two sets of weights u and v, and their convex combination w = αu + (1 -α)v. Then K(i, w) = p i G -1 w p i . By using G w = αG u + (1 -α)G v and the derivative of the inverse of a matrix, we find

d 2 K(i, w) dα 2 = 2 p i G -1 w (G u -G v ) G -1 w (G u -G v )G -1 w p i (33) = 2r G -1 w r ≥ 0 (34)
where r = (G u -G v )G -1 w p i and G -1 w is a positive matrix.

The inequality

K(w) ≤ K(u) α (35)
can also be proved by using the matrix inequality G -1 w ≤ G -1 u /α.

1) Subgradient of K(w): The subgradient of K(w) is the convex combination of the gradients of the K(i, w) such that K(i, w) = K(w). With similar computations as the previous results and Q the matrix of a basis orthogonal with respect to w, we find ∂K(i, w) ∂w j = -| q i , q j | 2 = -J ij (36)

The matrix J, as the Hadamard product of QQ and its conjugate, is semi-definite by the Schur product theorem. The subgradient of K(w) is the set of vectors -Jα, for α a positive vector of sum 1, and α i (K(i, w) -K(w)) = 0.

B. Karush-Kuhn-Tucker conditions

The Karush-Kuhn-Tucker conditions for the optimization problem [START_REF] Chauffert | Variable Density Sampling with Continuous Trajectories[END_REF] implies that there exists a scalar λ and a positive vector γ such that the subgradient of K(w) verifies

Jα = λ1 -γ (37)
with w i γ i = 0, appropriate conditions on α, and 1 a vector of L coefficients equal to 1. Combining the definitions of J, K(w), Q, and the properties of α, w and γ, we find:

• α Jw = K(w) • w Jα = λ • α Jα = λ -α γ = K(w) -α γ ≤ K(w)
• w Jw = m Using the symmetry and the positivity of J, we have :

|α Jw| 2 ≤ |w Jw||α Jα| (38) K(w) 2 ≤ mK(w) (39)
which implies K(w) ≤ m. Combined with the lower bound K(w) ≥ m, we have K(w) = m. By remarking that i w i K(i, w) = i w i q i q i = m from the orthogonality of the columns of Q, we also have that w i = 0 if K(i, w) < K(w).

Fig. 1 .

 1 Fig. 1. Measurement densities for a polynomial basis of size m = 20. The 2m largest values of the output of Newton's method are indicated by circles.

Fig. 2 .

 2 Fig. 2. Convergence of the proposed algorithm in various settings.

  Fig 4 shows the mean squared estimation error for different design methods, in the case of polynomial approximation in function of the number of samples. Polynomials of order 19 (i.e. m = 20) are used to approximate the function f

Fig. 4 .

 4 Fig. 4. Mean square error for polynomials approximation, m = 20.

Fig. 5 .

 5 Fig. 5. Optimized density for recovery of a wavefield using Fourier-Bessel functions, m = 21. (a) Frank-Wolfe algorithm. The density is here multiplied by m for an easier comparison with the maximal weight 1/m of a Dirac mass. (b) Weighted least-squares.

Fig. 6 .

 6 Fig. 6. Mean squared error for recovery using Fourier-Bessel functions, m = 21.

Fig. 7 .

 7 Fig. 7. Optimization of a sampling density for a polynomial basis and a gaussian matrix ensuring an upper bound on the error in the interval, m = 20.

Fig. 8 .

 8 Fig.8. K(w) for a polynomial basis and a cosine basis, with the jointly optimized density, with comparison to the mixed density and the lower and upper bounds.

Fig. 9 .

 9 Fig. 9. Mean square error for estimation using polynomial (top) and cosine (bottom) bases.

ACKNOWLEDGMENT

The author thanks Piotr Bojanowski for introducing him to the Frank-Wolfe algorithm.