L Pronzato
email: pronzato@i3s.unice.fr

A Zhigljavsky
email: zhigljavskyaa@cf.ac.uk

Algorithmic construction of optimal designs on compact sets for concave and differentiable criteria *

Keywords: Approximate design, optimum design, construction of optimal designs, global optimization MSC 62K05, 90C30, 90C26

We consider the problem of construction of optimal experimental designs (approximate theory) on a compact subset X of R d with nonempty interior, for a concave and Lipschitz differentiable design criterion φ(•) based on the information matrix. The proposed algorithm combines (a) convex optimization for the determination of optimal weights on a support set, (b) sequential updating of this support using local optimization, and (c) finding new support candidates using properties of the directional derivative of φ(•). The algorithm makes use of the compactness of X and relies on a finite grid X ℓ ⊂ X for checking optimality. By exploiting the Lipschitz continuity of the directional derivatives of φ(•), efficiency bounds on X are obtained and ǫ-optimality on X is guaranteed. The effectiveness of the method is illustrated on a series of examples.

Introduction and motivation

A design measure ξ on a finite set X ⊂ R d with ℓ elements is characterized by the ℓ-dimensional vector of weights W (nonnegative and summing to one) allocated to the ℓ elements of X . The determination of an optimal measure ξ * which maximizes a concave differentiable criterion φ(•) then forms a finite-dimensional convex problem for which many optimization algorithms are available, see, e.g., Hiriart-Urruty and [START_REF] Hiriart-Urruty | Convex Analysis and Minimization Algorithms, part 1 and 2[END_REF]; den [START_REF] Den Hertog | Interior Point Approach to Linear, Quadratic and Convex Programming[END_REF]; [START_REF] Nesterov | Interior-Point Polynomial Algorithms in Convex Programming[END_REF]; [START_REF] Ben-Tal | Lectures on Modern Convex Optimization: Analysis, Algorithms, and Engineering Applications[END_REF]; [START_REF] Boyd | Convex Optimization[END_REF]; [START_REF] Nesterov | Introductory Lectures to Convex Optimization: A Basic Course[END_REF] for recent developments on convex programming. In particular, the cutting plane method of [START_REF] Kelley | The cutting plane method for solving convex programs[END_REF] is considered by [START_REF] Sibson | Coefficients in D-optimal experimental design[END_REF], and a variant of it (with a modified version of the Equivalence Theorem) by [START_REF] Gribik | Equivalence theorems and cutting plane algorithms for a class of experimental design problems[END_REF]; see also Pronzato and Pázman (2013, Chap. 9). However, design problems are usually such that (a) the cardinality ℓ of the design space, which determines the dimension of the optimization problem to be solved, is large, and (b) there always exists an optimal measure ξ * with a few support points only, i.e., such that only a few components of W are positive (say m, with m ≪ ℓ). These particularities have motivated the development of specific methods which happen to be competitive compared to general-purpose convex-programming algorithms. One may distinguish three main families, see Pronzato and Pázman (2013, Chap. 9) for a recent survey.

(i) Vertex-direction methods only increase one component of W at each iteration, all other components being multiplied by the same factor smaller than one; see [START_REF] Wynn | The sequential generation of D-optimum experimental designs[END_REF]; [START_REF] Fedorov | Theory of Optimal Experiments[END_REF] and also [START_REF] Frank | An algorithm for quadratic programming[END_REF] for a method originally proposed in a more general context. A significant improvement can be achieved by decreasing some individual components of W only, which allows us to set some components of W to zero, i.e., to remove support points from a poorly chosen initial design and sometimes to exchange two components of W (vertex-exchange methods); see [START_REF] Atwood | Sequences converging to D-optimal designs of experiments[END_REF]; [START_REF] St | D-optimality for regression designs: a review[END_REF]; [START_REF] Böhning | Numerical estimation of a probability measure[END_REF][START_REF] Böhning | A vertex-exchange-method in D-optimal design theory[END_REF]; Molchanov andZuyev (2001, 2002).

(ii) In gradient methods the gradient direction projected on the set of design measures is used as the direction of movement at each iteration; see Wu (1978a,b), and [START_REF] Atwood | Convergent design sequences for sufficiently regular optimality criteria[END_REF] for an extension to a Newton-type method.

(iii) At each iteration of a multiplicative method, each component of W is multiplied by a suitably chosen positive function; see, e.g., [START_REF] Titterington | Algorithms for computing D-optimal designs on a finite design space[END_REF]; [START_REF] Silvey | An algorithm for optimal designs on a finite design space[END_REF]; [START_REF] Torsney | A moment inequality and monotonicity of an algorithm[END_REF]; [START_REF] Fellman | An empirical study of a class of iterative searches for optimal designs[END_REF]; [START_REF] Dette | Improving updating rules in multiplicative algorithms for computing D-optimal designs[END_REF]; Yu (2010a,b), and [START_REF] Torsney | W-iterations and ripples therefrom[END_REF] for a historical review.

In multiplicative methods all initial weights must be positive, those which should be zero at the optimum decrease continuously along iterations but typically stay strictly positive, i.e., never achieve exactly zero. The convergence is inevitably slow close to the optimum. The same phenomenon occurs for vertex-direction methods, unless the decrease of some particular individual components of W is allowed at some iterations, so that poor initial support points can be removed. Even in this case, if ℓ is large many iterations are required to identify the components of W which should be positive at the optimum. Gradient-type methods are also not efficient when ℓ is large.

Several authors tried combinations of different methods to make use of their respective advantages. A very sensible algorithm is proposed in (Wu, 1978a,b); it combines a gradient method (always working in a small dimensional subspace) with a vertex-direction method (which allows a suitable updating of this subspace). A mixture of multiplicative and vertex-exchange algorithms is proposed in [START_REF] Yu | D-optimal designs via a cocktail algorithm[END_REF] for D-optimum design; it includes a nearest-neighbor exchange strategy which helps apportioning weights between adjacent points in X and has the property that poor support points are quickly removed from the support of the initial measure.

Attractive performance is reported. All the methods above, however, are restricted to the case where the design space X is finite with l being not too large and, for some of them, to particular design criteria (D-optimality, for instance).

When one is interested in the determination of an optimal design on a compact subset X of R d with nonempty interior, the usual practice consists in discretizing X into a finite set X ℓ with ℓ elements and applying one of the methods above. When a precise solution is required, then ℓ is necessarily very large (in some cases, ℓ = 10 6 should be considered as a small number) and none of the methods above is efficient. Refining iteratively a finite grid contained in X is a possible option, see Wu (1978a), but the search for the optimal design is still performed in a discrete set.

In contrast, the algorithm we propose makes use of the compactness of X and looks for the optimal support in the whole set X . A finite grid X ℓ ⊂ X is only used to check optimality on X ℓ and, using the Lipschitz continuity of directional derivatives of φ(•), to construct an efficiency bound on X . A key ingredient in the algorithm is the separation between the determination of the support points of ξ * (knowing that there are at most m of them) from the determination of the associated weights (an m-dimensional convex problem).

The determination of the support of ξ * is a non-convex problem, usually multimodal, for which the straightforward application of a global search method (such as the simulated annealing or one of the genetic algorithms) cannot be recommended. Indeed, these heuristic methods do not use the crucial information about the objective function provided by the Lipschitz constants, derivatives and convexity of the weight-optimization problem, and most of them do not provide any indication of the closeness of the returned solution to an optimum, global or local. On the other hand, by using properties of the directional derivative of φ(•), we can easily locate good candidates for the support of ξ * and thereby we do not need to perform a global search in the m × d dimensional space X d . The aim of this paper is to show how these properties can be combined with algorithms for convex optimization which are classically used in design for the determination of optimal weights to yield an efficient method of construction of optimal designs on compact sets.

The paper is organized as follows. Section 2 defines the problem and introduces the notation.

Notation and problem statement

Let the design space X be a compact subset of R d with nonempty interior (typically X = [-1, 1] d ) and denote by Ξ(X ) the set of probability measures on X . Any element ξ ∈ Ξ(X ) will be called design measure, or shortly design.

Consider a model (the main example being a regression model) defined on X with p unknown parameters θ ∈ R p , p ≥ 2. Denote by M ≥ and M > the sets of p × p symmetric non-negative and positive definite matrices respectively. Let

M (ξ) = X µ(x) ξ(dx)
(1) be the information matrix related to the estimation of θ, with µ(x) ∈ M ≥ . The matrix µ(x) is the elementary information matrix for the estimation of θ associated with one single observation at x and may have rank larger than one (think of Bayesian optimal design or multivariate regression, see, e.g., Fedorov (1972, Chap. 5), [START_REF] Harman | Approximate D-optimal designs of experiments on the convex hull of a finite set of information matrices[END_REF]). In the case of simple linear regression model with observations y j = θ T f (x j ) + ε j , where the ε j denote i.i.d. random errors with zero mean, we have µ(x) = f (x)f T (x). We suppose that µ(•) is continuous on X (this assumption could be weakened, but it would make the presentation unduly complicated).

Although µ(•), and thus M (ξ), may depend on the model parameters θ in a nonlinear situation, this dependence will be omitted and, when it occurs, we only consider locally optimum designs and evaluate M (•) at a nominal value for θ. The set M X = {M (ξ) : ξ ∈ Ξ(X )} is a convex and compact subset of a linear space of dimension p(p + 1)/2 + 1; from Caratheodory's Theorem, see, e.g., Silvey (1980, p. 72), any element of M X can be written as M (ξ discr ), where ξ discr is a discrete measure supported on a finite set with m ≤ p(p + 1)/2 + 1 elements.

We are interested in constructing a φ-optimal design on X ; that is, a design ξ * ∈ Ξ(X ) which maximizes the functional

φ(ξ) = Φ[M (ξ)]
with respect to ξ ∈ Ξ(X ). The function Φ : M ≥ -→ R is normally assumed to be (see, e.g., Pukelsheim (1993, Chap. 5)):

(i) concave, i.e., Φ[(1 -α)M 1 + αM 2 ] ≥ (1 -α)Φ(M 1 ) + αΦ(M 2 ) for all M 1 , M 2 in M ≥ and all α in [0, 1]; (ii) proper, i.e., the set {M ∈ M ≥ : Φ(M ) > -∞} is nonempty and Φ(M ) < ∞ for all M ∈ M ≥ ;
(iii) isotonic (monotonic for Loewner ordering), i.e., Φ(M

1 ) ≥ Φ(M 2 ) for all M 1 , M 2 in M ≥ such that M 1 -M 2 is non-negative definite.
We suppose that a nonsingular design measure ξ 0 exists in Ξ(X ) (i.e., such that M (ξ 0 ) ∈ M > ). We shall denote by Φ + (•) the positively homogeneous version of Φ(•), which satisfies Φ + (I p ) = 1, with I p the p-dimensional identity matrix, and Φ + (aM ) = aΦ + (M ) for all M in M ≥ and all a positive. We also suppose that the design criterion Φ(•) is global, that is, its with Φ + (M ) = [(1/p)trace(M -1 )] -1 and ψ(t) = -p/t; see, e.g., Pukelsheim (1993, Chap. 5).

positively homogeneous version Φ + (•) satisfies Φ + (M ) > 0 if and only if M ∈ M > . We can write Φ + (M ) = ψ[Φ(M )],
The φ-efficiency of a design measure ξ ∈ Ξ(X ) is defined as

E φ (ξ) = φ + (ξ) max ξ∈Ξ(X) φ + (ξ) , (2)
where φ

+ (ξ) = Φ + [M (ξ)].
We assume that Φ(•) :

M ≥ -→ R is differentiable on M > with Lipschitz continuous gradient; that is, ∇Φ(M 1 ) -∇Φ(M 2 ) ≤ L(c) M 1 -M 2 (3) for all M 1 , M 2 in the set M(c) = {M ∈ M > : Φ(M ) ≥ c} ,
for some constant L(c), where ∇Φ(M ) denotes the gradient of Φ(•) at M , the norm • is the usual L 2 -norm (for any real matrix M , M = trace 1/2 (M ⊤ M )) and c is any number between inf ξ∈Ξ(X ) φ(ξ) and max ξ∈Ξ(X ) φ(ξ). We also suppose that

∇Φ(M ) ≤ B(c) for all M ∈ M(c) (4)
for some constant B(c). Note that L(c) and B(c) will only be used in the proofs of Theorems 2 and 3 and that their knowledge is not required for the construction of optimal designs with the algorithms of Sect. 3.

Isotonicity implies that any optimal information matrix M (ξ * ) lies on the boundary of the set M X , and therefore there always exists an optimal design measure with no more than m * (p) = p(p + 1)/2 support points. We can thus restrict our attention to such designs, i.e., to designs ξ defined by their support points {x (1) , . . . , x (m) }, m ≤ m * (p), and associated weights

w i = ξ({x (i) }) ≥ 0 which satisfy m i=1 w i = 1. The notation ξ = x (1) x (2) • • • x (m) w 1 w 2 • • • w m
is standard. One should notice that the bound m * (p) is often pessimistic since many situations exist where optimal designs are concentrated on a much smaller number of support points or are even saturated, i.e., are supported on p points only; see [START_REF] Yang | Support points of locally optimal designs for nonlinear models with two parameters[END_REF]; [START_REF] Yang | On de la Garza phenomenon[END_REF] and the paper [START_REF] Dette | A note on de la Garza phenomenon for locally optimal designs[END_REF] which makes use of results in (Karlin and Studden, 1966b) on Chebyshev systems and moment spaces.

For any n ≥ 1 and any design measure ξ on X = {x (1) , . . . , x (n) } ⊂ X n , we shall denote by W = W (ξ) = (w 1 , . . . , w n ) ⊤ ∈ P n the vector formed by the weights w i = {W } i = ξ(x (i) ), with P n the probability simplex

P n = {W = (w 1 , . . . , w n ) ⊤ ∈ R n : w i ≥ 0 , n i=1 w i = 1} . (5)
Some components w i of W may equal zero, and we shall denote by W + (ξ) the vector formed by strictly positive weights and by S = S(ξ) the support of ξ, defined by the corresponding design points. With a slight abuse of notation, we shall consider X and S sometimes as sets with respectively n and m ≤ n elements, and sometimes as vectors in R n×d and R m×d . Also, we shall respectively denote by φ(S|W + ) and φ(W + |S) = φ(W |X) the value of the design criterion φ(ξ) when considered as a function of S with W + fixed and as a function of W + with S fixed, and by ∇φ(W + |S) the usual gradient ∂φ(W + |S)/∂W + (a m-dimensional vector). For any ν and ξ in Ξ(X ) such that φ(ξ) > -∞, we can compute the directional derivative of φ(•) at ξ in the direction ν,

F (ξ; ν) = lim α→0 + φ[(1 -α)ξ + αν] -φ(ξ) α = trace{[M (ν) -M (ξ)]∇Φ[M (ξ)]} (with F (ξ; ν) ∈ R ∪ {+∞}).
When ν is the delta-measure δ x at x, we obtain

F (ξ, x) = F (ξ; δ x ) = trace{[µ(x) -M (ξ)]∇Φ[M (ξ)]} . (6)
Suppose that ξ and ν are finitely supported and let X denote the union of their supports.

Denoting by W and W ′ the vectors of weights which ξ and ν respectively allocate to points in X, with some components of W and W ′ possibly equal to zero, we obtain

F (ξ; ν) = (W ′ -W ) ⊤ ∇φ(W |X) .
The concavity and differentiability of Φ(•) yield a necessary and sufficient condition for the optimality of a design measure ξ * , see, e.g., [START_REF] Silvey | Optimal Design[END_REF], usually named the Equivalence Theorem as a tribute to the work of [START_REF] Kiefer | The equivalence of two extremum problems[END_REF]; see also Karlin and Studden (1966a).

Theorem 1 The following properties are equivalent (i) ξ * is φ-optimal on X , i.e., ξ * maximizes φ(ξ) with respect to ξ ∈ Ξ(X );

(ii) ξ * minimizes max x∈X F (ξ, x) with respect to ξ ∈ Ξ(X );

(iii) max x∈X F (ξ * , x) = 0.
Note that F (ξ; ξ) = 0 for all ξ ∈ Ξ(X ) and that F (ξ; ν) = X F (ξ, x) ν(dx) for all ξ and ν ∈ Ξ(X ), so that (iii) implies that ξ * -almost everywhere we have F (ξ * , x) = 0; that is,

F (ξ * , x) = 0 on the support of ξ * .
We shall say that ξ ǫ is ǫ-optimal for φ(•) on X if and only if

max x∈X F (ξ ǫ , x) < ǫ .
The concavity of φ(•) then implies that φ(ξ ǫ ) > max ξ∈Ξ(X) φ(ξ) -ǫ.

The directional derivative of the positively homogeneous criterion φ

+ (•) = ψ[φ(•)] at ξ in the direction ν is F φ + (ξ; ν) = ψ ′ [φ(ξ)]F (ξ; ν) . We denote F φ + (ξ, x) = F φ + (ξ; δ x ) for all x ∈ X
and all ξ ∈ Ξ(X ). Due to the concavity of φ + (•), the ǫ-optimality of ξ ǫ for φ(•) implies that

φ + (ξ ǫ ) > max ξ∈Ξ(X) φ + (ξ) -ǫ ψ ′ [φ(ξ ǫ )]; this yields the efficiency bound E φ (ξ ǫ ) > 1 - ǫ ψ ′ [φ(ξ ǫ )] max ξ∈Ξ(X) φ + (ξ) ≥ E φ (ξ ǫ ) = 1 - ǫ ψ ′ [φ(ξ ǫ )] φ + (ξ ǫ ) , (7)
where the φ-efficiency E φ (•) is defined in (2).

Algorithms

We denote by A 0 a convex optimization algorithm for the determination of optimal weights.

When initialized at an arbitrary nonsingular design on X, for any ǫ > 0, the algorithm A 0 is assumed to return an ǫ-optimal design in a finite number of iterations; we shall denote this ǫ-optimal design by ξ

= A 0 [X, ǫ].
In the examples of Sect. 4, we use a combination of the projected-gradient and vertex-exchange methods, following ideas similar to those in (Wu, 1978a), but other methods could be used as well, see a discussion in Sect. 1.

Algorithm A 0

Input: Discrete set X = {x (1) , . . . , x (m) } ∈ X m (such that M (ν 0 ) has full rank, with ν 0 the uniform measure on X) and ǫ > 0.

Output:

Design ξ = A 0 [X, ǫ] such that max x∈X F (ξ, x) < ǫ.

Construction of the main algorithm

We shall construct an algorithm for the maximization of φ(•) on the set of design measures on the compact set X . We denote by ξ k the design measure obtained at iteration k (always with finite support), by S k = S(ξ k ) its support and by W + k the vector of associated weights, and write

φ k = φ(ξ k ).
The construction is decomposed into three steps. First, we present a prototype algorithm A 1 ; at each iteration k of this algorithm we determine an (ǫ/2)-optimal design ξ k = A 0 [X k , ǫ/2] for an increasing sequence of sets X k , where X k has one element more than X k-1 . The main difference between algorithm A 1 and the vertex-direction methods discussed in Sect. 1 is step 1, where the weights at the current support are optimized. Although closely related, the algorithm also differs from a direct application of the cutting-plane method as considered by [START_REF] Sibson | Coefficients in D-optimal experimental design[END_REF]; [START_REF] Gribik | Equivalence theorems and cutting plane algorithms for a class of experimental design problems[END_REF].

Algorithm A 2 is very similar to A 1 but in A 2 the number of points in the sets X k is bounded.

In both A 1 and A 2 we assume that max x∈X F (ξ k , x) can be determined with arbitrary precision ǫ. In practice, we usually perform the search for the maximum of F (ξ k , x) over a finite test-set X ℓ . Using Lipschitz continuity arguments and (7), we can then guarantee some efficiency bound over X when X ℓ ⊂ X is a fine enough grid, see Sect. 3.3. However, the computational effort required becomes too high when ℓ gets large if we search for the maximum of F (ξ k , x) over the whole grid X ℓ at each iteration. In algorithm A 3 , we use a two-level strategy where the maximization of F (ξ k , •) is performed over a coarse test-set T k which is progressively enriched by points taken from X ℓ . Moreover, the maximization over T k is complemented by a local search over X .

The addition of a local maximization of φ(S|W + k ) with respect to the support S of ξ k at each iteration preserves the convergence of φ k to the optimum φ * = φ(ξ * ) and yields our final algorithm A 4 , which also incorporates a merging strategy which avoids having clusters of points (that is, groups of points positioned closely) in T k and X k . This is important as the presence of such clusters slows down the convergence of algorithms.

Algorithm A 1

Step 0) Choose X 0 = {x (1) , . . . , x (m) } ∈ X m such that M (ν 0 ) has full rank, with ν 0 the uniform measure on X 0 ; choose ǫ > 0, set k = 0.

Step

1) Compute ξ k = A 0 [X k , ǫ/2].
Step 2) Find x * k = arg max x∈X F (ξ k , x).

Step 3

) If F (ξ k , x * k ) < ǫ, stop; otherwise, set X k+1 = {X k , x * k }, k ← k + 1, go to step 1.
Note that we can always use m ≤ p(p + 1)/2 at step 0 above. The next theorem establishes an important property of algorithm A 1 .

Theorem 2 Algorithm A 1 stops after a finite number of iterations. The design measure ξ k obtained at step 3 is then ǫ-optimal, in the sense

φ k = φ(ξ k ) > φ * -ǫ.
Proof. Denote by φ * k the maximum value of φ(•) on the set of probability measures supported on X k . From the (ǫ/2)-optimality of the design ξ k returned by algorithm A 0 , we have

φ 0 = φ(ξ 0 ) ≥ φ * 0 -ǫ/2. Since X 0 ⊂ X k for any k ≥ 1, φ * k ≥ φ * 0 and φ k ≥ φ * k -ǫ/2 ≥ φ * 0 -ǫ/2 for all k. Now, for all k and all x, x ′ ∈ X we have |F (ξ k , x ′ ) -F (ξ k , x)| = trace{[µ(x ′ ) -µ(x)]∇Φ[M (ξ k )] ≤ µ(x ′ ) -µ(x) ∇Φ[M (ξ k )] ≤ µ(x ′ ) -µ(x) B(φ * 0 -ǫ/2) ,
see ( 4) and ( 6). Since X is compact and µ(x) is continuous, there exists ρ > 0 such that for

any x and x ′ in X , x ′ -x < ρ implies that |F (ξ k , x ′ ) -F (ξ k , x)| < ǫ/2 for all k. Since ξ k = A 0 [X k , ǫ/2], F (ξ k , x (i) ) < ǫ/2 for all x (i) in X k , so that any x in a ball B(x (i) , ρ) with centre x (i)
and radius ρ is such that F (ξ k , x) < ǫ. This implies that x * k constructed at step 2 is at distance larger than ρ from all points in X k when algorithm A 1 does not stop at step 3. In view of the fact that X is compact, A 1 necessarily stops after a finite number of iterations.

The concavity of φ(•) implies that φ * ≤ φ k + max x∈X F (ξ k , x) and thus φ * < φ k + ǫ when the algorithm stops.

Remark 1 Note that the algorithm stops in less than k * iterations, with

k * = max{k : there exists (x 1 , . . . , x k ) ∈ X k such that min i =j x i -x j ≥ ρ} .
An upper bound on k * (very pessimistic) can easily be constructed when X is an hypercube in R d , say [-1, 1] d , using packing radius. Indeed, the existence of (x 1 , . . . , x k ) in the definition of k * implies that k non-intersecting balls of radius R = ρ/2 can be packed in the hypercube

[-(1 + ρ/2), (1 + ρ/2)] d , so that k * ≤ (2 + 4/ρ) d /V d , with V d = π d/2 /Γ(d/2 + 1) the volume of the d-dimensional unit ball.
In algorithm A 1 , the size of X k increases by one at each iteration (see step 3), so that the dimension of the optimization problem to be solved by A 0 at step 1 is also increasing. In order to facilitate the work of A 0 we shall consider another construction for X k .

The support S k = S(ξ k ) of the design ξ k determined at step 1 may be strictly included in

X k . It is then tempting to replace step 3 of A 1 by Step 3 ′ ) If F (ξ k , x * k ) < ǫ stop; otherwise set X k+1 = {S k , x * k }, k ← k + 1, go to step 1.
However, the arguments used in the proof of Theorem 2 for the non-existence of cluster points for the sequence {x * k } are no longer valid since X k does not contain all previous points x * k ′ , k ′ < k. We show in Theorem 3 that convergence is ensured by choosing a suitably decreasing sequence of constants

γ k in ξ k = A 0 [X k , γ k ǫ] at step 1, e.g., γ k = 1/(k + 1) or γ k = γ k for some γ ∈ (0, 1
). The algorithm is then as follows.

Algorithm A 2 Let {γ k }, k ≥ 0, denote a decreasing sequence of numbers in (0, 1), with lim k→∞ γ k = 0.

Step 0) Same as in A 1 .

Step 1) Compute

ξ k = A 0 [X k , γ k ǫ].
Step 2) Same as in A 1 .

Step 3) Use step 3 ′ above.

Theorem 3 Algorithm A 2 stops after a finite number of iterations and the design measure ξ k obtained at step 3 is ǫ-optimal.

Proof. We use the same notation as in the proof of Theorem 2. The maximum φ * k of φ(•) in the space of probability measures supported on X k is also the maximum of φ(•) for a measure

supported on S k . Since S k ⊂ X k+1 , φ * k+1 ≥ φ * k and therefore φ k > φ * k -γ k ǫ ≥ φ * 0 -γ k ǫ ≥ φ * 0 -ǫ for all k ≥ 0.
We suppose that A 2 never stops and show that we reach a contradiction. Denote

ξ ′ k (α) = (1-α)ξ k +αδ x * k , with δ x * k the delta measure allocating mass 1 at x * k . The function α -→ φ[ξ ′ k (α)
] is concave, reaches its maximum for some α k ∈ [0, 1) and equals φ k for α = 0 and α = α ′ k , with α ′ k some number larger than α k . For any α ∈ [0, α ′ k ] we may write

φ[ξ ′ k (α)] = φ k + α trace{[µ(x * k ) -M (ξ k )]∇Φ[M (ξ ′ k (α ′ ))]} for some α ′ ∈ [0, α]. Using (3) and the fact that φ k > φ * 0 -ǫ, we obtain φ[ξ ′ k (α)] = φ k + αtrace{[µ(x * k ) -M (ξ k )]∇Φ[M (ξ k )]} +α trace{[µ(x * k ) -M (ξ k )] [∇Φ[M (ξ ′ k (α ′ ))] -∇Φ[M (ξ k )]]} ≥ φ k + αF (ξ k , x * k ) -α µ(x * k ) -M (ξ k ) ∇Φ[M (ξ ′ k (α ′ ))] -∇Φ[M (ξ k )] ≥ φ k + αF (ξ k , x * k ) -α α ′ µ(x * k ) -M (ξ k ) 2 L(φ * 0 -ǫ) ≥ φ k + αǫ -α 2 C L(φ * 0 -ǫ) with C = max x,x ′ ∈X µ(x ′ ) -µ(x) 2
, where we used the property

F (ξ k , x * k ) ≥ ǫ. Therefore, φ * k+1 ≥ max α φ[ξ ′ k (α)] ≥ φ k + ǫ 2 4C L(φ * 0 -ǫ)
and

φ k+1 ≥ φ * k+1 -γ k+1 ǫ ≥ φ k + ǫ 2 4C L(φ * 0 -ǫ) -γ k+1 ǫ . (8)
This implies that there exists some k 0 such that, for all k > k 0 ,

φ k+1 > φ k + ǫ 2 /[8C L(φ * 0 -ǫ)].
The sequence {φ k } is thus unbounded, which contradicts the assumptions on X , µ(•) and φ(•).

Algorithm A 2 thus stops in a finite number of iterations. As in the proof of Theorem 2, the concavity of φ(•) implies that ξ k is then ǫ-optimal.

Remark 2 An upper bound k * on the number of iterations required by A 2 can easily be derived from (8). Indeed, k * satisfies k * ǫ 2 /A -ǫ k * i=1 γ i ≤ φ * -φ 0 , where we have denoted A = 4C L(φ * 0 -ǫ) and φ * = φ(ξ * ). This gives k * ≤ (A/ǫ 2 )[φ * -φ 0 + γǫ/(1 -γ)] for γ i = γ i with γ ∈ (0, 1), and k * ≤ A 2 /(4ǫ 2 )[1 + 1 + 4(φ * -φ 0 )/A] 2 for γ i = 1/(i + 1) (where we have used the property

(1/k) k i=1 γ i < 1/ √ k).
Step 2 of algorithms A 1 and A 2 may be difficult to implement when X has nonempty interior since the function x ∈ X -→ F (ξ k , x) is severely multimodal. We thus discretize X into a finite grid X ℓ and use the maximizer over this grid as initialization for a local maximization of

F (ξ k , •) over X .
The grid must be fine enough (i) to ensure that the maximum of F (ξ k , •) is not missed and (ii) to guarantee a reasonable level of optimality over X , see Sect. 3.3. However, the computational cost corresponding to the evaluation of F (ξ k , x) for all points of X ℓ at each iteration would be too high when ℓ is very large. For that reason, algorithm A 3 presented below uses a two-level strategy, with a test-set T k remaining much smaller than X ℓ .

We denote by x * = LM [F (ξ, •); x 0 ] the result of a local maximization of F (ξ, x) with respect to x ∈ X initialized at x 0 , see Sect. 3.2 for possible algorithms when X is the hypercube [-1, 1] d , or the probability simplex P d given by (5) (mixture experiments). As in A 2 , {γ k } denotes a decreasing sequence of numbers in (0, 1), with lim k→∞ γ k = 0.

Algorithm A 3

Step 0) Choose X 0 = {x (1) , . . . , x (m) } ∈ X m such that M (ν 0 ) has full rank, with ν 0 the uniform measure on X 0 ; choose some arbitrary test-set T 0 of n points in X ; choose ǫ > 0, set k = 0.

Step 1) Compute

ξ k = A 0 [X k , γ k ǫ/2].
Step 2)

Find x * = arg max x∈T k F (ξ k , x) and compute x * k = LM [F (ξ k , •); x * ].
Step 3) If F (ξ k , x * k ) ≥ ǫ/2, go to step 7; otherwise go to step 4.

Step 4) Find x * * = arg max x∈X ℓ F (ξ k , x).

Step 5) If F (ξ k , x * * ) < ǫ/2, stop.

Step

6) Compute x * k = LM [F (ξ k , •); x * * ].
Step 7)

Set X k+1 = {S k , x * k }, T k+1 = (T k , x * k ), k ← k + 1; go to step 1.
The test-set T 0 of step 0 is arbitrary, one may choose for instance T 0 = X 0 or the 2 d vectors with components ±1 when X = [-1, 1] d . The algorithm repeatedly goes through steps 1 to 3 and 7 until F (ξ k , x) becomes smaller than ǫ/2 for all x ∈ T k ; then the algorithm goes through step 5 which tests the (ǫ/2)-optimality of ξ k over X ℓ . The proof of convergence uses the same arguments as for A 2 : the algorithm stops after a finite number of iterations and returns a design measure ξ k that is (ǫ/2)-optimal on X ℓ . The set X ℓ can be chosen such that (ǫ/2)-optimality on X ℓ implies ǫ-optimality on the whole compact set X , see Sect. 3.3.

The following two easy improvements can be made to A 3 .

(i) When min x∈S k x -x * k < δ k at step 7, with δ k some small number, one may consider replacing the point in S k closest to x * k by x * k and accept the substitution if φ(ξ k ) < φ(ξ ′ k ), with ξ ′ k the design measure after the substitution. We denote the result of this operation by

X k+1 = EX[{S k , x * k }; δ k ].
(ii)

Step 1 can be complemented by a local maximization of φ(S|W + k ) with respect to S ∈ X m k , initialized at S k , with m k the number of support points of ξ k , S k their collection and W + k the associated weights. We shall denote the result of this local maximization by

S ′ k = LM [φ(•|W + k ); S k ]; see Sect. 3.2 for algorithms.
Since (i) and (ii) above can only increase the value of φ k , the arguments used in the proof of Theorem 3 remains valid and the convergence to an (ǫ/2)-optimal measure on X ℓ is preserved.

We also add modifications which may decrease the value of φ(•) and for which some caution is thus required.

(iii) If the substitution is rejected in

X k+1 = EX[{S k , x * k }; δ k ], see (i),
then the set X k+1 may contain points which are very close, which may complicate the task of A 0 at step 1. We thus consider aggregation of points in X k+1 and replace (sequentially) all pairs of points x i = x j such that x i -x j < δ k by (x i +x j )/2, until the set X ′ k+1 which is obtained is such that min x i =x j ∈X ′ k+1 x i -x j ≥ δ k . We denote this operation by X ′ k+1 = AG(X k+1 , δ k ).

(iv) When the set T k+1 contains many points, say more than some specified number N , one may consider an aggregation step similar to (iii) above and use AG(T k+1 , δ k ).

To preserve convergence, we take {δ k } as a positive decreasing sequence tending to zero as k → ∞ (e.g., δ k = δ/(k + 1) for some δ > 0) and, to avoid possible oscillations which may be caused by alternating local optimization in (ii) and aggregation of points in (iii), the local optimization is only used when there is an increase of φ(•) in two successive passages through step 1. The final algorithm is as follows. Remember that W + k denotes the vector formed by the strictly positive weights of ξ k (in practice we only keep weights larger than some small constant, typically 10 -6 , all other points are removed from ξ k and their total weight is reallocated to remaining points proportionally to the gradient components).

Algorithm A 4 Let {γ k } and {δ k }, k ≥ 0, denote two decreasing sequences of numbers in (0, 1), with lim k→∞ γ k = lim k→∞ δ k = 0.

Step 0) Choose X 0 = {x (1) , . . . , x (m) } ∈ X m such that M (ν 0 ) has full rank, with ν 0 the uniform measure on X 0 ; choose some arbitrary test-set T 0 of n points in X ; choose N > n and ǫ > 0, set k = 0, φ old = φ(ν 0 ).

Step

1) Compute ξ k = A 0 [X k , γ k ǫ/2]. Step 1a) Compute S ′ k = LM [φ(•|W + k ); S k ]. Step 1b) If φ(S ′ k |W + k ) > φ old , substitute S ′ k for S k as support of ξ k . Step 1c) Compute φ k = φ(ξ k ), set φ old = φ k .
Step 2)

Find x * = arg max x∈T k F (ξ k , x) and compute x * k = LM [F (ξ k , •); x * ].
Step 3) If F (ξ k , x * k ) ≥ ǫ/2, go to step 7; otherwise go to step 4.

Step 4) Find x * * = arg max x∈X ℓ F (ξ k , x).

Step 5) If F (ξ k , x * * ) < ǫ/2, stop.

Step

6) Compute x * k = LM [F (ξ k , •); x * * ]. Step 7) Set X * = EX[{S k , x * k }; δ k ] and T * = (T k , x * k ).
Step 7a) If #T * > N , aggregate points in T * : T k+1 = AG(T * , δ k ); otherwise set

T k+1 = T * .
Step 7b) Aggregate points in X * : X k+1 = AG(X * , δ k ), k ← k + 1, go to step 1.

In practice the finite set X l can be taken as an adaptive grid such that, as the number ℓ of points tends to infinity, the minimax distance ρ = max y∈X min x∈X ℓ y -x q tends to zero, for some L q -norm • q . The search for x * * = arg max x∈X ℓ F (ξ k , x) of step 4 is then replaced by a sequential inspection of grid points: the screening through the grid is stopped when either (i) a grid point x i is found such that F (ξ k , x) ≥ ǫ/2, this point being then taken as x * * , or (ii) the value of ρ is small enough to ensure that max x∈X F (ξ k , x) < ǫ, that is, ξ k is ǫ-optimal on X . This is detailed in Sect. 3.4.

Local maximizations

Steps 1a and 2 involve local maximizations, respectively with respect to the m k points forming the support S k of ξ k , S k ⊂ X m k , and with respect to x ∈ X . There exist many sophisticated constrained optimization methods (such as sequential quadratic programming for instance) which can solve this problem efficiently. Note that very high accuracy is not required, any improvement compared with the initialization is enough to ensure convergence of algorithm A4.

Simple unconstrained optimization methods can therefore also be used when X has a simple form. Some methods are suggested below for situations where X is the hypercube [-1, 1] d and for the case of mixture experiments with X given by the simplex P d defined by (5).

Remark 3 Note that a local maximization of F (ξ k , x) initialized at the support points of ξ k is generally not adequate at step 2. Consider for instance D-optimal design for the linear regression

model η(x, θ) = θ 1 + θ 2 x sin(x) on X = [0, x], with x ≥ 2. For ξ = (1/2)δ 0 + (1/2)δ 2 , F (ξ, 0) = F (ξ, 2) = 2 and F (ξ, x
) is locally maximum at x = 0 and x = 2 (with ξ being Doptimal on X when 2 ≤ x ≤ π), but max x∈X F (ξ, x) > 2, and ξ is not D-optimal, when x > π.

Reparameterization When X = [-1, 1] d , one may consider the reparameterization x = cos(y) ∈ X (to be understood componentwise) with a local maximization with respect to y ∈ R d at step 2, or with respect to the m k points y (1) , . . . , y (m k ) , each one in R d , at step 1a; see Appendix B for further details. Any classical method for unconstrained optimization can then be used, including derivative-free methods, see for instance [START_REF] Powell | An efficient method for finding the minimum of a function of several variables without calculating derivatives[END_REF].

Projected gradient Let h(•) denote a differentiable function on some convex and compact set A , with gradient ∇h(•), and consider the local maximization of h(x) with respect to x ∈ A .

The projected-gradient algorithm uses iterations of the form

x j+1 = P A x + j (α * )
with x 0 given in A , where x + j (α) = x j + α ∇h(x j ),

α * = arg max α h(P A [x + j (α)]) and P A [•] denotes the orthogonal projection onto A . When X = [-1, 1] d , P X [x]
simply amounts to a truncation of all components of x to [-1, 1].

The value α * can be obtained by a classical line-search algorithm, restricted to α in the interval [0, α max ], with

α max = min min i=1,...,d:{∇h(x j )} i ≥0 (1 -{x j } i )/{∇h(x j )} i , min i=1,...,d:{∇h(x j )} i ≤0 (1 + {x j } i )/|{∇h(x j )} i | ,
which guarantees that x + j (α) ∈ X . The same method can be applied when the optimization is with respect to S k ∈ X m k ; projections on X m k then correspond to m k projections over X .

The iterations can be stopped when the norm of the projected gradient is smaller than some specified threshold.

The projected-gradient method can also be used in the case of mixture experiments where X is the simplex P d given by ( 5). The projection of x + j (α) on X can then be obtained as the solution of a Quadratic Programming (QP) problem: P P d [x + j (α)] minimizes x -x + j (α) 2 with respect to x satisfying the linear constraints x ∈ P d . Alternatively, one may also use the property that, for any Z ∈ R d , the projection P P d (Z) is given by x(Z, t * ), where x(Z, t) = max{Z -t1 d , 0 d } (componentwise, with 1 d and 0 d respectively the d-dimensional vectors of ones and zeros) and t * maximizes L[x(Z, t), t] with respect to t, with L(x, t) the partial Lagrangian

L(x, t) = (1/2) x -Z 2 + t [1 ⊤ d x -1]
. Indeed, L(x, t) can be written as

L(x, t) = (1/2) x -(Z -t1 d ) 2 + t [1 ⊤ d Z -1] -d t 2 /2 ,
which reaches its minimum with respect to x ∈ R d+ for x = x(Z, t). One may notice that

max{max i {Z} i -t * , 0} ≤ 1 ⊤ d x(Z, t * ) = 1 ≤ d max{max i {Z} i -t * , 0}
, so that the search for t * can be restricted to the interval [max i {Z} i -1, max i {Z} i -1/d], see Pronzato and Pázman (2013, Chap. 9).

From (ǫ/2)-optimality on a grid to ǫ-optimality on a compact set

For any ξ ∈ Ξ(X ) and any x, y ∈ X we have

|F (ξ, x) -F (ξ, y)| = |trace{[µ(x) -µ(y)]∇Φ[M (ξ)]}| ≤ µ(x) -µ(y) ∇Φ[M (ξ)] .
Only the second term on the right-hand side depends on the criterion φ(•) and design ξ, the first term being related to the problem under consideration. For instance (assuming that M (ξ) has full rank), D-optimality with φ

(ξ) = log det[M (ξ)] gives ∇Φ[M (ξ)] = M -1 (ξ), A-optimality with φ(ξ) = -trace[M -1 (ξ)] gives ∇Φ[M (ξ)] = M -2 (ξ)
and, more generally, when t ≥ 1 and A is non- (Harville, 1997, Chap. 15).

negative definite φ A,t (ξ) = -trace[AM -t (ξ)] gives ∇Φ[M (ξ)] = t i=1 M -i (ξ)A ⊤ M -(t+1-i) (ξ), see
Suppose that µ(•) satisfies the Lipschitz condition ∀y ∈ X ∩ B q (x, ρ) , µ(y) -µ(x) ≤ L µ (x, ρ) ρ , with B q (x, ρ) the ball {y : y -x q ≤ ρ} for some norm • q , 1 ≤ q ≤ ∞. We then obtain

∀y ∈ X ∩ B q (x, ρ) , F (ξ, y) ≤ F (ξ, x) + ρ L F (ξ, x, ρ) , (9)
where

L F (ξ, x, ρ) = L µ (x, ρ) ∇Φ[M (ξ)] . Suppose in particular that µ(•) is differentiable on
X which is a convex set. Then, for any x, y in X we can write

F (ξ, y) = F (ξ, x) + ∂F (ξ, z) ∂z ⊤ z=x+α(y-x) (y -x)
for some α ∈ [0, 1], so that ( 9) is satisfied for q = 2 with

L F (ξ, x, ρ) = max max z∈B 2 (x,ρ) ∂F (ξ, z) ∂z .
More generally, we can use regularity properties of µ(•) to derive inequalities of the form

∀y ∈ X ∩ B q (x, ρ) , F (ξ, y) ≤ F (ξ, x) + h(ξ, x, ρ) , (10)
where, for any ξ such that M (ξ) has full rank, sup x∈X h(ξ, x, ρ) → 0 as ρ → 0.

Consider then the final design ξ k returned by algorithm A 4 and take ρ = max y∈X min x∈X ℓ y-

x q with X ℓ the set used at step 4 of the algorithm. Then (10) implies that

max y∈X F (ξ, y) ≤ max x∈X ℓ {F (ξ k , x) + h(ξ, x, ρ)}
and, if ρ is small enough to ensure that max x∈X ℓ h(ξ k , x, ρ) < ǫ/2, ξ k is ǫ-optimal on X and the efficiency bound (7) applies.

The case where µ(x) has rank one, i.e., µ(x) = f (x)f ⊤ (x) for some p-dimensional regressor vector f (x), deserves special attention. Consider φ 0 (ξ) = log det[M (ξ)] (D-optimality) and

φ t (ξ) = -trace[M -t (ξ)] for t ≥ 1.
For each of these criteria we get

|F (ξ, y) -F (ξ, x)| = c t f ⊤ (y)M -(t+1) (ξ)f (y) -f ⊤ (x)M -(t+1) (ξ)f (x) = c t ∆f ⊤ M -(t+1) (ξ)∆f + 2∆f ⊤ M -(t+1) (ξ)f (x) ≤ c t ∆f 2 λ t+1 min [M (ξ)] + 2 c t f (x) ∆f λ t+1 min [M (ξ)] , (11)
where we have denoted c t = max(t, 1) and ∆f = f (y) -f (x). This gives (10) with

h(ξ, x, ρ) = c t λ t+1 min [M (ξ)] D 2 ρ (x) + 2D ρ (x) f (x) ( 12
)
where

D ρ (x) = max y∈X ∩Bq(x,ρ) f (y) -f (x) . If f (•) is differentiable, then we also have ∂F (ξ, x)/∂x = 2 c t f ⊤ (x)M -(t+1) (ξ) ∂f (x)/∂x so that ∂F (ξ, x) ∂x ≤ 2 c t f (x) ∂f (x) ∂x 1 λ t+1 min [M (ξ)]
.

This gives ( 9) with

L F (ξ, x, ρ) = {2c t /λ t+1 min [M (ξ)]} max y∈X ∩B 2 (x,ρ) { f (y) ∂f (y)/∂y }.
Any particular problem can thus be handled by a case-by-case analysis in order to obtain a Lipschitz inequality similar to (9) or (10). Examples with polynomial regression models are considered in Appendix A.

Adaptive grids

An adaptive construction of the grid X l used at step 4 of algorithm A 4 allows us to only refine the grid at locations where the upper bound (10) needs to be improved. Meanwhile, it also sometimes speeds up the search for x * * . Here we indicate a possible construction when

X = [-1, 1] d .
For any n ≥ 0, let X 1,mm (n) = {x (1) , . . . , x (n) } denote the design (sometimes called minimaxdistance optimal) given by x (j) = (2j -1)/n -1, j = 1, . . . , n, and

X d,mm (n) = X d 1,mm (n) denote the n d -point d-dimensional grid with all coordinates in X 1,mm (n).
Step 4 of algorithm A 4 is then replaced by the following (step 5 is removed).

Step 4-0) Choose n (say, n = 100), N max (say, N max = 10 7 ), set

X (0) ℓ = X d,mm (n), G 0 = n d , ρ 0 = 1/m, j = 0. Step 4-1) Find x * * = arg max x∈X (j) ℓ F (ξ k , x). Step 4-2) If F (ξ k , x * * ) ≥ ǫ/2, go to step 6 of A 4 . p support points 2 -1 1 3 -1 0 1 4 -1 -1 √ 5 1 √ 5 ≃ 0.4472 1 5 -1 - √ 3 √ 7 0 √ 3 √ 7 ≃ 0.6547 1 6 -1 - √ 7+2 √ 7 √ 21 - √ 7-2 √ 7 √ 21 √ 7-2 √ 7 √ 21 ≃ 0.2852 √ 7+2 √ 7 √ 21 ≃ 0.7651 1
Table 1: D-optimal designs for polynomial regression on [-1, 1] (all support points are equally weighted).

Step 4-3)

If ǫ = max x∈X (j) ℓ F (ξ k , x) + h(ξ k , x, ρ j ) < ǫ or if G j > N max , stop A 4 : ξ k is ǫ-optimal on X .
Step 4-4) Set X

(j>) ℓ = {x (j) ∈ X (j) ℓ : F (ξ k , x (j) ) + h(ξ k , x (j) , ρ j ) ≥ ǫ}; for each x (j) ∈ X (j>) ℓ , divide the hypercube B ∞ (x (j) , ρ j ) into 2 d hypercubes B ∞ (x (j i ) , ρ j /2), where x (j i ) = x (j) + ρ j x (i) d,mm (2) with x (i) d,mm (2) the i-th element of the 2 d -point grid X d,mm (2), i = 1, . . . , 2 d . Collect all x (j i ) to form X (j+1) ℓ , set G j+1 = #X (j+1) ℓ , ρ j+1 = ρ j /2, j ← j + 1, go to step 4-1.
Steps 4-3 and 4-4 rely on (10); although the construction of the grids X (j) ℓ use the norm • ∞ , another norm can be used in the definition of h(ξ, x, ρ). Due to the bound N max set on the number G j of elements of X (j) ℓ , the algorithm may stop before the precision on max x∈X F (ξ k , x) reaches ǫ, i.e., one may have ǫ > ǫ at step 4-3. The use of a precise bound in (10) is then crucial to avoid a fast increase of G j and a premature stopping of the algorithm.

Examples

Description

We have tested algorithm A 4 on a series of examples with known optimal design ξ * . The design space is X = [-1, 1] d with d = 1 or 2 and the number of parameters p varies between 3 and 9. We use φ(ξ) = log det M (ξ) for D-optimal design and φ(ξ) = -trace[M -1 (ξ)] for A-optimal design.

Problems 1 to 4 correspond to D-optimal design for linear regression models with f (x) given by ( 13) and p = s varying from 3 to 6. For each p, the D-optimal design measure is uniquely defined and gives weight 1/p at the roots of the polynomial t -→ (1 -t 2 )P ′ p-1 (t) with P ′ k (•) the derivative of the k-th Legendre polynomial, see Table 1.

Problems 5 and 6 correspond to D-optimal design for additive polynomial models: f (x) is given by ( 14) with d = 2, p = 2s -1 and s = 3, 4 respectively. Problem 7 correspond to Doptimal design for a complete product-type interaction model with d = 2: f (x) = g(x 1 ) ⊗ g(x 2 ) with ⊗ denoting tensor product and g(x) given by ( 13) with s = 3; this gives p = 9. For those three problems the D-optimal design is the cross-product of the D-optimal design measures for the individual models; see Schwabe (1996, Chap. 4 and 5).

In Problems 8 and 9, the models are the same as in Problems 1 and 7 but the criterion is φ(ξ) = -trace[M -1 (ξ)] (A-optimality). The A-optimal design for Problem 8 is

ξ * A = -1 0 1 1/4 1/2 1/4 ,
the A-optimal design for Problem 9 is the product measure ξ * A ⊗ξ * A ; see Schwabe (1996, Chap. 4). Problem 10 corresponds to (local) D-optimal design for the nonlinear regression model

η(x, θ) = θ 0 + θ 1 exp(-θ 2 x 1 ) + θ 3 θ 3 -θ 4 [exp(-θ 4 x 2 ) -exp(-θ 3 x 2 )]
with five parameters θ = (θ 0 , θ 1 , θ 2 , θ 3 , θ 4 ) and two design variables

x = (x 1 , x 2 ) ∈ X = [0, 2] × [0, 10].
A numerical value must be given to the parameters θ 2 , θ 3 and θ 4 that intervene nonlinearly in η(x, θ), and we use θ 2 = 2, θ 3 = 0.7, θ 4 = 0.2. The model considered is additive with a constant term, so that the tensor product of the D-optimal designs for the two models

η 1 (x 1 , β (1) ) = β (1) 0 + β (1) 1 exp(-β (1) 2 x 1 ) and η 2 (x 2 , β (2) ) = β (2) 0 + β (2) 1 [exp(-β (2) 2 x 2 ) - exp(-β (2) 1 x 2 )]/(β (2) 1 -β (2)
2 ), respectively at β (1) = (θ 0 , θ 1 , θ 2 ) and β (2) = (θ 0 , θ 3 , θ 4 ), is D-optimal for η(x, θ) at θ, see [START_REF] Schwabe | Designing experiments for additive nonlinear models[END_REF]. These two D-optimal designs are supported on three points only, each one receiving mass 1/3, which can be computed with arbitrary precision using Theorem 1-(iii): the support points are approximately 0, 0.46268527927 and 2 for η 1 (x 1 , β (1) ) and 0, 1.22947139883, 6.85768905493 for η 2 (x 2 , β (2) ). The design space X is renormalized to [-1, 1] 2 in order to use the adaptive-grid construction of Sect. 3.4; the tests for optimality are based on (10) and (12).

Results and discussion

The results obtained are presented in Table 2. When d = 1, the algorithm is initialized with X 0 having m equidistant points, X 0 = {(2i -1)/m -1, i = 1 . . . , m}, with m = p(p + 1)/2; when d = 2, X 0 is the tensor product of two such designs, with m 2 the smallest integer larger than or equal to p(p + 1)/2. The test-set T 0 is the union of X 0 and the 2 d vertices of X . The number of points in T k always remains reasonably small (see Table 2) so that step 7a is not used.

The decreasing sequences {γ k } and {δ k } are respectively given by 1/(k + 1) and 0.1/(k + 1), ǫ is set to 10 -6 in all examples. The algorithm A 0 used at step 1 is based on a combination of projected-gradient and vertex-exchange methods, see Wu (1978a). Local maximizations at steps 1a and 2 of A 4 use reparameterization and the derivative-free algorithm of [START_REF] Powell | An efficient method for finding the minimum of a function of several variables without calculating derivatives[END_REF].

The set X ℓ of step 4 corresponds to an adaptive grid, constructed according the the algorithm of Sect. 3.4 with n = 100 for d = 1 and n = 10 for d = 2 (so that G 0 = 100).

In all problems considered the optimal design ξ * is unique and one can thus compute its distance to ξ kmax returned by A 4 . Various metrics can be considered and we use here the Wasserstein and Lévy-Prokhorov metrics, see Appendix C. Table 2 indicates that ξ kmax is very close to ξ * for all problems considered. We also report in Table 2 the distance to 1 of the efficiencies E φ (ξ kmax ) defined by (2) together with the values of 1 -E φ (ξ kmax ), with E φ (ξ kmax ) the efficiency bound given by ( 7). Note that 1 -E φ (ξ) < ǫ/p for D-optimality when max x∈X F (ξ, x) < ǫ and φ(ξ) = log det M (ξ) and 1-E φ (ξ) < ǫ/trace[M -1 (ξ)] for A-optimality when max x∈X F (ξ, x) < ǫ and φ(ξ) = -trace[M -1 (ξ)]. Due to the fact that the efficiency bound, used to terminate the algorithm, is very pessimistic in most problems, the true efficiency of ξ kmax is often much closer to 1 than the bound indicated. We thus observe in all cases considered that a very good precision is achieved although the number of iterations k max is small. Also note that the construction of the adaptive grid is only required a small number of times (only once for more than half of the problems considered); due to the local maximizations performed at steps 1a, 2 and 6, the points in the test-set T k , although there are only few of them, are generally enough to construct an optimal design on the whole set X .

Most of the computational cost of A 4 corresponds to step 1 and especially step 4.

Concerning step 4 (which we try to avoid as much as possible), the use of an adaptive grid, as proposed in Sect. 3.4, aims at minimizing its cost. Using precise bounds in the construction of h(ξ, x, δ) used in ( 10) is essential to maintain the cardinality G j of X (j) ℓ reasonably small. One may notice that the construction of Sect. 3.4 does not make use of all the information contained in the evaluations of F (ξ k , x i ) for the various x i which are used. Indeed, some x (j) are removed from X (j) ℓ when constructing the set X (j>) ℓ

. With some computational effort, other constructions making a more efficient use of the information collected on F (ξ k , •) could be considered, and the results in [START_REF] Harman | Improvements on removing non-optimal support points in D-optimum design algorithms[END_REF][START_REF] Pronzato | A delimitation of the support of optimal designs for Kiefer's φ p -class of criteria[END_REF] could be used to remove parts of X with low values of F (ξ k , x) that cannot contain support points of an optimal design. Also, Lipschitz bounds obtained from diagonal partitions, see [START_REF] Sergeyev | Global search based on efficient diagonal partitions and a set of Lipschitz constants[END_REF], might be used to get a significant reduction of the number of evaluations of F (ξ k , x i ) required when d > 1. Of course, for a given problem this number depends very much on the value of ǫ, see the bottom part of Table 2.

Concerning step 1, the algorithm is constructed in such a way that X k has a small number of elements. The choice of the optimization method used for A 0 at step 1 is then not crucial. We have used a combination of the projected-gradient and vertex-exchange methods for the results in Table 2. Numerical experimentations with a pure projected-gradient algorithm [START_REF] Mccormick | The gradient projection method under mild differentiability conditions[END_REF], the methods of cutting planes [START_REF] Kelley | The cutting plane method for solving convex programs[END_REF] and the level method (Nesterov, 2004, Sect. 3.3.3), originally developed for non-differentiable optimization (see also Pronzato and Pázman (2013, Chap. 9)), yield similar results.

Conclusion

We have proposed an algorithm for the construction of optimal design measures over a compact set X for a concave and differentiable criterion φ(•). The method exploits the property that optimal designs are often supported on a small number of points and combines updating of the support with convex optimization of the weights it receives. Efficiency bounds and guaranteed ǫ-optimality over X are obtained using Lipschitz continuity of the directional derivative of φ(•). As illustrated by some examples, the algorithm seems to be quite effective when accurate Lipschitz constants can be determined for the particular problem considered. General techniques for deriving those constants are given which could be applied to any model satisfying usual regularity conditions. Also, some general ideas have been given for the local optimization of the support of a design measure, which should allow the construction of specific methods especially adapted to particular design regions X .

On the other hand, some issues remain that call for further developments. We mention two of them.

(i) There exist problems for which the optimal design is not unique; in such situations it would be of interest to have an algorithm capable of ensuring convergence to an ǫ-optimal design with minimal support. This is not the case of our algorithm A 4 and improvements in that direction are under current investigation. Note that once an optimal design is found, the determination of another optimal design (having necessarily the same information matrix when Φ(•) is strictly concave) with minimum support corresponds to a fixed-charge problem, see, e.g., [START_REF] Sadagopan | A vertex ranking algorithm for the fixed-charge transportation problem[END_REF]. A simple reduction of the support of an optimal design can be obtained by iteratively removing elementary matrices µ(x (i) ) that are in the convex hull of others µ(x (j) ), j = i.

(ii) Although the convex optimization of weights for non-differentiable criteria (for instance, E-optimality) is not a big challenge, see, e.g., Nesterov (2004, Sect. 3.3.3), Pronzato and Pázman (2013, Chap. 9), the determination of new points being candidates for inclusion in the support and the construction of efficiency bounds over X raise specific issues due to the fact that the maximum of the directional derivative of φ(•) is not always attained at a one-point (delta) measure. The geometrical characterization of optimal designs, see in particular [START_REF] Dette | Geometry of E-optimality[END_REF], may then prove salutary for the construction of efficient algorithms. Pb. 8 10 -6 1 3 10 1 9.9 10 -9 8.9 10 -5 5.2 10 -9 < 10 -16 8.4 10 -8 8 196 9 10 -6 2 9 13 1 1.0 10 -8 5.7 10 -5 2.2 10 -8 1.3 10 -15 1.2 10 -8 53 ≃ 4.9 10 6 10 10 -6 2 5 9 1 2.8 10 -7 2.1 10 -4 3.6 10 -8 4.8 10 -14 1.8 10 -7 21 ≃ 6.4 10 6 10 10 -5 2 5 7 1 1.1 10 -6 7.7 10 -4 1.5 10 -6 1.9 10 -12 1.9 10 -6 21 ≃ 1.7 10 6 10 10 -4 2 5 7 1 5.0 10 -7 1.9 10 -6 1.6 10 -6 1.6 10 -12 1.3 10 -5 21 ≃ 5.9 10 5 10 10 -3 2 5 6 1 4.9 10 -6 7.6 10 -6 1.6 10 -5 1.7 10 -10 1.0 10 -4 21 ≃ 2.2 10 5 10 10 -2 2 5 4 1 1.6 10 -4 2.4 10 -4 6.1 10 -4 1.8 10 -7 1.5 10 -3 21 ≃ 8.7 10 4

ǫ d p k max # step 4 W 1 W 2 L 1 -E φ (ξ kmax ) 1 -E φ #T kmax G max 1 10 -6 1
Table 2: Results obtained by using algorithm A 4 on Problems 1-10 with ǫ = 10 -6 and on Problem 10 with ǫ between 10 -2 and 10 -6 : number k max of iterations before stopping and number of passages through step 4; Wasserstein distances W 1 (ξ kmax , ξ * ) and W 2 (ξ kmax , ξ * ) and Lévy-Prokhorov distance L(ξ kmax , ξ * ) between ξ kmax and the optimal design ξ * ; (1-) φ-efficiency E φ (ξ kmax ) and (1-) efficiency bound E φ (ξ kmax ) given by ( 7); number of elements of T kmax and maximum number of points in the grids X (j) ℓ used to guarantee the ǫ-optimality of ξ kmax .

Appendix A: Bounds for directional derivatives in polynomial regression

In this appendix, we construct bound for f (x) and f (y) -f (x) to be used in (12) for polynomial regression models. Here X = [-1, 1] and µ

(x) = f (x)f ⊤ (x) with f (x) = (1, x, x 2 , . . . , x s-1 ) ⊤ , s ≥ 2 . ( 13
)
We then have the global Lipschitz bound µ(x) -µ(y) ≤ A s |x -y|, with A s the s × s matrix

A s =           0 1 2 • • • s -1 1 2 2 . . . . . . . . . s -1 • • • 2(s -1)          
so that A s 2 = trace(A ⊤ s A s ) = s 2 (7s -5)(s -1)/6 . Also, for all x, y ∈ X , ∆f 2 = f (y)f (x) 2 ≤ (x -y) 2 [1 + 2 2 + • • • + (s -1) 2 ] = (x -y) 2 s(s -1)(2s -1)/6 , and max x∈X f (x) = √ s, which can be used as global bounds in (11). More precise local bounds can also be obtained: where ω s (c) = s-1 j=1 j 2 c 2(j-1) . Therefore,

D δ (x) = max y∈X ∩B 1 (x,δ)
f (y) -f (x) ≤ δ ω s (max{|x + δ|, |x -δ|}) .

to be used in (10,12).

Multi-factor polynomial regression models can be handled similarly. For additive models with f (x) = (1, x 1 , x 2 1 , . . . , x s-1 1 , x 2 , x 2 2 , . . . , x d , x 2 d , . . . , x s-1 d ) ⊤ , (14)

x = (x 1 , . . . , x d ) ∈ X = [-1, 1] d , then ∆f 2 ≤ x-y 2 s(s-1)(2s-1)/6 and max x∈X f (x) =

[1 + d(s -1)] 1/2 . We also obtain, similarly to the case d = 1 considered above, .

Consider now models with complete product-type interactions and take f (x) = g(x 1 ) ⊗ g(x 2 ) with g(x) given by ( 13), x = (x 1 , x 2 ) ∈ X = [-1, 1] 2 . Then ∆f = f (y)-f (x) = [g(x 1 )⊗g(x 2 )]-[g(y 1 ) ⊗ g(y 2 )] = [g(x 1 ) -g(y 1 )] ⊗ g(x 2 ) + g(y 1 ) ⊗ [g(x 2 ) -g(y 2 )] and ∆f ≤ [g(x 1 ) -g(y 1 )] ⊗ g(x 2 ) + g(y 1 ) ⊗ [g(x 2 ) -g(y 2 )] = [g(x 1 ) -g(y 1 )]g ⊤ (x 2 ) + g(y 1 )[g(x 2 ) -g(y 2 )] ⊤ = g(x 1 ) -g(y 1 ) g(x 2 ) + g(y 1 ) g(x 2 ) -g(y 2 ) ≤ s(s -1)(2s -1)

6 1/2 √ s [|x 1 -y 1 | + |x 2 -y 2 |] ≤ s (s -1)(2s -1) √ 3 x -y , f (x) = g(x 1 ) g(x 2 ) ≤ s .
We also get the local bounds where D(y) is the diagonal matrix diag{sin(y i ), i = 1, . . . , d}. Any vector y such that, for all i = 1, . . . d, y i = k i π with k i ∈ Z is thus a stationary solution, i.e., the reparameterization renders any point on the boundary of X stationary. However, a local optimization algorithm applied to h(•) provides a local maximizer for h(•). Indeed, suppose that the maximization of h(•) with respect to y ∈ R d yields a vector y * such that sin(y * i ) = 0 for i = 1, . . . , q and sin(y * i ) = 0 for i = q + 1, . . . , d .

Since h(•) is maximized, one may assume that ∂ 2 h(y) ∂y∂y ⊤ y=y * is semi-negative definite and (ŷ(i) -y * ) ⊤ ∂ h(y) ∂y y=ŷ(i) < 0 , ŷ(i) = y * + α i e i , i = 1, . . . , d,

with e i the i-th basis vector, for |α i | small enough. Denote x * = cos(y * ) (notice that x * ∈ X ) and take any i ∈ {q + 1, . . . , d}; (15) implies that ∂h(x)/∂x i x=x * > 0 when y * i = 2kπ (i.e., x * i = 1) and ∂h(x)/∂x i x=x * < 0 when y * i = (2k + 1)π (i.e., x * i = -1), k ∈ Z. Therefore, there exist d coefficients λ i ≥ 0 and λ ′ i ≥ 0 such that the Lagrangian

L (x, λ, λ ′ ) = h(x) + d i=1 λ i (1 -x i ) + d i=1 λ ′ i (1 + x i )
satisfies ∂L (x, λ, λ ′ )/∂x x=x * = 0 with λ i (1 -x i ) = λ ′ i (1 + x i ) = 0 for all i (take λ i = λ ′ i = 0 for i = 1, . . . , q; λ i = ∂h(x)/∂x i x=x * , λ ′ i = 0, when x * i = 1 and λ i = 0, λ ′ i = -∂h(x)/∂x i x=x * , when x * i = -1 for i = q +1, . . . , d). Moreover, for any z ∈ R d such that z i = 0 for i = q +1, . . . , d,

z ⊤ ∂ 2 L (x, λ, λ ′ ) ∂x∂x ⊤ x=x * z = z ⊤ D -1 1...q (y * )
∂ 2 h(y) ∂y∂y ⊤ y=y * 1...q, 1...q D -1 1...q (y * ) z ≤ 0 .

The inequality is strict when the matrix given by the first q rows and q columns of the Hessian matrix ∂ 2 h(y)/(∂y∂y ⊤ ) y=y * is negative definite, which is a sufficient condition for the (local) optimality of x * , see, e.g., Luenberger (1973, p. 235).

  Section 3 presents the algorithm and proves its convergence. A few illustrative examples are given in Sect. 4 where the results obtained with the proposed algorithm are discussed. Section 5 concludes and indicates some remaining open issues. A few technical aspects are collected in Appendix.

  with ψ(•) the inverse function of the transformation a ∈ R + -→ Φ(aI p ). Note that ψ(•) is increasing. We suppose that Φ + (•) is concave. Standard examples are Φ(M ) = log det(M ) with Φ + (M ) = det 1/p (M ), ψ(t) = exp(t/p), and Φ(M ) = -trace(M -1 )

  -y) 2 ω s (max{|x|, |y|})

  i -y i ) 2 ω s (max{|x i |, |y i |}) so that D δ (x) = max y∈X ∩B∞(x,δ) f (y) -f (x) ≤ δ d i=1 ω s (max{|x i + δ|, |x i -δ|})

D

  δ (x) = max y∈X ∩B∞(x,δ) f (y) -f (x) ≤ δ min g(x 2 ) ω s (max{|x 1 + δ|, |x 1 -δ|}) + √ s ω s (max{|x 2 + δ|, |x 2 -δ|}), √ s ω s (max{|x 1 + δ|, |x 1 -δ|}) + g(x 1 ) ω s (max{|x 2 + δ|, |x 2 -δ|}) .Appendix B: Reparameterization for optimization on an hypercubeLet h(•) denote a function to be maximized with respect tox ∈ X = [-1, 1] d , d ≥ 1. We suppose that h(•) is twice continuously differentiable in X .Consider the reparameterization defined by x = x(y) = cos(y) (componentwise), with y ∈ R d , and denote h(y) = h[x(y)], so that

Acknowledgments

The authors thank the referees for useful comments that helped to significantly improve the paper.

* This work was initiated while the authors were invited at the Isaac Newton Institute for Mathematical Sciences, ; the support of the INI and of CNRS is gratefully acknowledged.

Appendix C: measures of performance

In situations where the optimal design is known, to measure the performance of our algorithm we can naturally compare the value of the design criterion φ kmax = φ(ξ kmax ), with ξ kmax the design measure obtained when the algorithm stops, to the optimal value φ * = max ξ∈Ξ(X ) φ(ξ). When there exists a unique optimal measure ξ * ∈ Ξ(X ) such that φ(ξ * ) = φ * , it is also instructive to evaluate how close ξ kmax is to ξ * .

For the results of Sect. 4 we use the Wasserstein metric W p (ξ kmax , ξ * ) with p = 1 and p = 2 and the LévyProkhorov metric L(ξ kmax , ξ * ); both define metrizations for the topology of weak convergence.

For ν 0 and ν 1 two probablity measures on X , the Wasserstein metric W p (ν 0 , ν 1 ) is defined by

with Ω(ν 0 , ν 1 ) the set of all couplings of ν 0 and ν 1 , i.e., of all measures ω on X × X having ν 0 and ν 1 as marginals. When d = 1 we have

where IF -1 i denotes the quantile function IF -1 i (t) = inf{x : IF i (x) ≥ t} with IF i the distribution function for ν i , i = 0, 1, see Major (1978, Theorem 8.1).

The Lévy-Prokhorov metric L(ν 0 , ν 1 ) is defined by

. Note that W p (ν 0 , ν 1 ) and L(ν 0 , ν 1 ) are sensitive to the scaling of X ; in particular, W p (ν ′ 0 , ν ′ 1 ) = α W p (ν 0 , ν 1 ) when ν ′ 0 and ν ′ 1 are obtained from ν 0 and ν 1 by applying the homothetic transformation x ∈ X -→ z = α x, α > 0.

Both metrics are difficult to compute for general probability measures but the fact that the measures ξ * and ξ kmax have finite support makes the situation much simpler. Computational details will be presented elsewhere.