Jean-Luc Dominique

Xavier

Laurent, Tiphaine Stefan Laetitia

Germán J'

Malgré l'impossibilité de tenir dans une page le vécu de toute cette période, je tâcherai d'exprimer ma gratitude à tous ceux qui m'ont accompagné pendant ces travaux de thèse de la façon la plus exhaustive possible.

Je commence par un grand merci à mes directeurs de thèse, Anne et Marcel. Leur enthousiasme et excellente disposition ont beaucoup contribué à l'accomplissement de cette thèse. Je tiens à remercier en même temps Jean-Baptiste Hiriart-Urruty, qui a suivi de près l'avancement de cette thèse, en apportant toujours de conseils et du soutien, et Felipe Alvarez, qui a guidé mes premiers pas dans le monde de la recherche et dont le soutien a été déterminant pour continuer mes études en France. Pareillement, c'est avec grand plaisir que je salue la convivialité de l'équipe Optimisation et Interactions, et plus particulièrement d'Aude Rondepierre, avec qui j'ai partagé bureau pendant une très agréable première année de thèse.

J'ai également l'honneur de remercier les professeurs Fabio Schoen, Michael Schyns et Roy Welsch pour avoir accepté de rapporter cette thèse et de s'être, par surcroît, déplaces malgré leur emploi de temps chargés pour participer à mon jury de thèse. Toujours sur le plan académique, une pensée amicale pour Matías Salibián-Barrera et Rubén Zamar, qui m'ont accueilli pour un séjour de recherche au Canada. Il a été un grand plaisir de partager avec eux et profiter de leur expérience.

Outre le développement professionnel visé en priorité par le doctorat, j'en ai aussi tiré un enrichissement personnel d'une valeur inestimable. Y ont contribué, tout d'abord, Benjamin, à qui je dois son amitié et patience envers mon français hésitant du début.

Résumé

La statistique robuste est une branche de la statistique qui s'intéresse à l'analyse de données contenant une proportion significative d'observations contaminées avec des erreurs dont l'ampleur et la structure peuvent être arbitraires. Les estimateurs robustes au sens du point de rupture sont généralement définis comme le minimum global d'une certaine measure non-convexe des erreurs, leur calcul est donc un problème d'optimisation globale très coûteux. L'objectif de cette thèse est d'étudier les contributions possibles des méthodes d'optimisation globale modernes à l'étude de cette classe de problèmes.

La première partie de la thèse est consacrée au τ -estimateur pour la régression linéaire robuste, qui est défini comme étant un minimum global d'une fonction nonconvexe et dérivable. Nous étudions l'impact des techniques d'agglomération et des conditions d'arrêt sur l'efficacité des algorithmes existants. Les conséquences de certains phénomènes liés au voisin le plus proche en grande dimension sur ces algorithmes agglomératifs d'optimisation globale sont aussi mises en évidence.

Dans la deuxième partie de la thèse, nous étudions des algorithmes déterministes pour le calcul de l'estimateur de moindres carrés tronqués, qui est défini à l'aide d'un programme en nombres entiers non linéaire. En raison de sa nature combinatoire, nous avons dirigé nos efforts vers l'obtention des bornes inférieures pouvant être utilisées dans un algorithme du type branch-and-bound. Plus précisément, nous proposons une relaxation par un programme sur le cône de deuxième ordre, qui peut être renforcée avec des coupes dont nous présentons l'expression explicite. Nous fournissons également des conditions d'optimalité globale.

iii

Abstract

Robust statistics is a branch of statistics dealing with the analysis of data containing contaminated observations. The robustness of an estimator is measured notably by means of the breakdown point. High-breakdown point estimators are usually defined as global minima of a non-convex scale of the errors, hence their computation is a challenging global optimization problem. The objective of this dissertation is to investigate the potential contributions of modern global optimization methods to this class of problems.

The first part of this thesis is devoted to the τ -estimator for linear regression, which is defined as a global minimum of a nonconvex differentiable function. We investigate the impact of incorporating clustering techniques and stopping conditions in existing stochastic algorithms. The consequences of some phenomena involving the nearest neighbor in high dimension on clustering global optimization algorithms is thoroughly discussed as well.

The second part is devoted to deterministic algorithms for computing the least trimmed squares regression estimator, which is defined through a nonlinear mixedinteger program. Due to the combinatorial nature of this problem, we concentrated on obtaining lower bounds to be used in a branch-and-bound algorithm. In particular, we propose a second-order cone relaxation that can be complemented with concavity cuts that we obtain explicitly. Global optimality conditions are also provided.

Lo terrible se aprende enseguida, y lo hermoso nos cuesta la vida...

Presentation of the statistical problem

Robust statistics is a branch of statistics dealing with the analysis of data that may contain large portions of contaminated observations. The origin of such contamination can be very diverse. In Figure 1, we show an example where the contamination comes from measurement errors [START_REF] Rousseeuw | Robust regression and outlier detection[END_REF]. For each year from 1950 to 1973, the number of outgoing international phone calls from Belgium are plotted. The bulk of the data follows a linear model; nonetheless, there are 6 observations that deviate from the majority. In fact, during the period between 1964 and 1969, there was a change on the recording system, which actually recorded the total duration, in minutes, of the international phone calls instead of the number of calls. A case where a subpopulation acts differently is shown in Figure 2. There are plotted, for 32 chemical compounds, a quantity called the krafft point versus a molecular descriptor called heat of formation. There is a main group that follows a regression line, but there is also, besides some few outliers at right, a second, smaller group forming what seems to be another regression line. The observations in the second group correspond to sulfonates. Therefore, due to diversity of the possible situations, it is impossible to make assumptions on the magnitude nor on the structure of the contamination.

Besides robustness, in a sense to be specified soon, robust regression estimators satisfy statistical properties such as asymptotic normality, square-root rate of convergence and equivariance properties. Though, the use of robust estimators is not as widespread as one could expect; mainly because their computation is very time consuming and, unlike many other problems arising in statistics, the difficult problems involved in computing robust estimators are almost unknown to optimization specialists. The objective of this thesis is to explore the potential improvements that can be obtained in the computing efficiency by taking advantage of recent advances in optimization.

Robust statistics is nowadays a well-established discipline, and most classical statistical procedures have a robust counterpart. However, it is customary to consider separately the case of multivariate analysis from linear regression. In many circumstances the linear model assumption leads to work in lower dimension than in multivariate analysis.

The regression context

The linear regression context is the following:

• Observed sample {x 1 , ..., x n } ⊆ R d ;

• response y = (y 1 , ..., y n ) ∈ R d ;

• linear model assumption

y i = x i β + δ i ,
with δ i independent and identically distributed , E[δ i ] = 0, V ar[δ i ] = σ 2 ;

• The observations are rows of a matrix X ∈ R n×d ,

X =        x 1 x 2 . . . . . . x n        =         x (1) 1 x (2) 1 • • • x (d) 1
x

(1) 2

x

(2) 2 x

(2) n

• • • x (d) n         ;
• for β ∈ R d , we denote by r(β) the vector of residuals r = y -Xβ, with components r i = y i -x i β.

Our objective is to obtain estimators β(X, y) of the parameter β satisfying additionally some equivariance properties:

β(X, λy + Xγ) = λ β(X, y) + γ for all λ ∈ R, γ ∈ R d ,
and for all nonsingular d × d matrix A:

β(XA, y) = A -1 β(X, y). Classical estimators satisfy these properties, but break down in presence of contamination. Let us take as an example the well-known least squares estimator, which is defined as βLS = Arg min where

S LS (r) = n i=1 r 2 i = r 2 2 . (1)
We can see in Figure 3 the rout of the LS estimator on a dataset with 20% of contamination. It has been observed that the quadratic scale in (1) overweights observations with large residuals. Then it seems reasonable to replace the quadratic scale by a scale with linear growth that should result in an estimator less sensitive to outlying observations. Such an estimator is the so-called (with a slight abuse of notation) L 1 estimator, defined as the minimizer of the 1 norm of residuals,

βL 1 = Arg min β∈R d r( β) 1 ,
The example of Figure 4, presented in [START_REF] Clarke | Optimization and nonsmooth analysis[END_REF] as an application of non smooth optimization in statistics, confirms somehow this intuition; however, the 1 regression over the dataset of Figure 3 does not give better results than the least squares regression, as shown in Figure 5. These contradictory signals raise the need for a clear definition of what «robust» means. This will be the subject of the next section.

A measure of robustness: the Breakdown Point

As the previous examples about the L 1 estimator highlight, it is necessary to formalize the idea of robustness. Hereafter we will adopt the robustness notion introduced by [START_REF] Donoho | The notion of breakdown point[END_REF], which is based on the concept of breakdown point (BDP). Roughly speaking, the BDP of an estimator on a sample is defined as the minimum fraction of the observations that need to be replaced by arbitrary ones for Figure 5: 1 regression on the same dataset with a 20% of contamination as in Figure 3.

the estimator to take on arbitrary values. The concept of BDP is consistent with the requirement that no hypothesis can be made on the distribution of the contamination, as it only depends on the fraction of the sample that is to be corrupted. The formal definition of the BDP is the following: consider any sample of n points (X, y) and let T be a regression estimator, T (X, y) = β. Then consider all possible corrupted samples (X , y ) that are obtained by replacing m of the original points by arbitrary values. Let us denote by bias(m; T, X, y) the maximum bias that can be caused by such a contamination bias(m; T, X, y) = sup X ,y T (X, y) -T (X , y ) 2 .

(2)

The breakdown point of the estimator T at the sample (X, y) is defined as * n (T, X, y) = min m n | bias(m; T, X, y) = ∞ .

In order to compare different estimators, one usually considers the asymptotic behaviour of * n (T, X, y): * (T, X, y) = lim n→∞ * n (T, X, y).

So, as * n ( βLS , X, y) = * n ( βL 1 , X, y) = 1/n for any sample (X, y), one says that both the LS and the L 1 estimators have a breakdown point of 0%.

Note that in equation (2) contamination of any type and magnitude is allowed. This is the main difference with the approach of Robust Optimization, which searches for an optimal solution to problems involving uncertain data on the objective function f 0 and/or on the constraints f i , for which one can specify an uncertainty set U : min

x {f 0 (x, ζ) | f i (x, ζ) ≤ 0, i = 1, ..., m. ∀ζ ∈ U}.
In general, robust optimization counterparts to classical estimators will not have a higher breakdown point than the original estimator. However, in problems where the uncertainty can be confidently supposed to belong to well-defined (convex!) uncertainty sets, the robust optimization framework could be preferable, since there exists efficient polynomial-time algorithms for solving those problems. The interested reader is referred to the article by [START_REF] Ben-Tal | Selected topics in robust convex optimization[END_REF] for further details.

The BDP of the L 1 estimator

As the example in Figure 5 suggests, * ( βL 1 , X, y) = 0 if contamination is allowed in X. However, since the L 1 estimator can be efficiently computed, there has been some interest on its BDP robustness in the case of planned experiments, namely, when the matrix of regressors X is deterministic and only y can be contaminated. In this case there exists some characterizations of * ( βL 1 , X, y), which can be positive. Let m * be the largest integer such that for any set S ⊆ {1, ..., n} of cardinality m * ,

max ξ =1 i / ∈S |x i ξ| |x i ξ| ≥ 1 2 ,
then * ( βL 1 , X) = (m * + 1)/n (see [START_REF] He | Tail behavior of regression estimators and their breakdown points[END_REF][START_REF] Mizera | The influence of the design on the breakdown points of 1 -type m-estimators[END_REF]. The combinatorial nature of this characterization makes the BDP of the L 1 estimator very difficult to compute, thus unmanageable. [START_REF] Giloni | The finite sample breakdown point of 1 -regression[END_REF] gave another characterization based on the linear structure of the problem, which involves solving a Mixed-Integer Program (MIP) with 2n integer variables. Nevertheless, the intuition behind the example of Figure 4 is not completely without foundation. [START_REF] Candes | Decoding by linear programming[END_REF] showed the following: assume that we want to recover the vector β from corrupted measurements y = Xβ + e, for a deterministic matrix X and an arbitrary unknown vector of errors e. Then β is the unique solution to min

g∈R d y -Xg 1 ,
provided that the 0 "norm", or cardinality norm of e, e 0 = cardinality({i | e i = 0}) is less than or equal to some η > 0. The number η must satisfy δ η + θ η,η + θ η,2η < 1, where

δ η = max |J|≤η,c∈R |J| F J c 2 c 2 -1 , θ η,η = max |J|≤η,|J |≤η ;c∈R |J| ;c ∈R |J | | F J c, F J c | c c ,
and F is any matrix F such that F X = 0. For an index set J, F J denotes the sub matrix of the rows of F indexed by J. Since the result holds under the hypothesis that uncontaminated observations have null residuals, this framework is rather adapted to signal recovery. It is not surprising that in the example of Figure 4 the majority of the observations pass exactly through a line. Going even deeper into the theoretical basis of this results we can find a result due to [START_REF] Fazel | Matrix rank minimization with applications[END_REF] stating that the 1 norm is indeed the convex envelope of the cardinality norm 0 over the unit ball. In simple words, the 1 norm is the one that gives the least weight to largeresiduals observations while keeping the problem convex. In the following section we will describe a class of estimators with a high breakdown point. As it could be expected, they involve nonconvex minimization problems.

High BDP estimators

Robust estimators are defined similarly to the LS and L 1 estimators, but the quadratic and linear scales of residuals are replaced by robust residual scales, which are not influenced by observations with huge residuals. We present below two such robust scales, M -estimates of scale and L-estimates of scale. The former defines robust estimators whose computation involve continuous nonconvex optimization problems, and the later gives raise to mixed-integer nonlinear programs (MINLPs).

Robust estimators based on L-estimates of scale

a i |r| (i) , or n i=1 a i |r| 2 (i) 1/2 , (3)
where the a i s are nonnegative constants.

The least trimmed squares (LTS) estimator It is defined by minimizing an L-scale as in the second form in (3), with a i = 1 for 1 ≤ i ≤ h and a i = 0 for h < i ≤ n, for a given h. In other words, the LT S estimator is defined as

βLTS = Arg min β∈R d S LT S (r( β)), (4)
where

S LT S (r) = h i=1 |r| 2 (i) 1/2
.

The breakdown point of the LTS estimator depends on h, and the optimal choice is h = n/2 + (p + 1)/2 which leads to a BDP of nearly 50% [START_REF] Rousseeuw | Robust regression and outlier detection[END_REF].

From an optimization viewpoint, problem (4) can be written as the following MINLP: [START_REF] Rousseeuw | Computing lts regression for large data sets[END_REF] describe a stochastic algorithm for approximating the LTS estimator. It has been implemented and incorporated into the robustbase package in R. The only known deterministic approach is that of [START_REF] Giloni | Least trimmed squares regression, least median squares regression, and mathematical programming[END_REF], where the authors study the polytope defined by the constraints w i ∈ {0, 1}, w i ≥ h and show that the relaxed problem attains the minimum at an extreme point, which necessarily is 0 -1. Then the authors propose a heuristic algorithm.

min β∈R d w i ∈{0,1} n i=1 w i r i ( β) 2 1/2 (5) subject to the constraint n i=1 w i ≥ h.
The least median of squares (LMS) estimator The LMS estimator is defined as a minimizer of the median of the residuals:

βLMS = Arg min β∈R d median i=1,...,n |r| (i) = Arg min β∈R d |r| (h) ,
where h = n/2 . It minimizes an L-scale as in the first form of (3), with a i = 0 except for a h which equals 1. The breakdown point of the LMS estimator is the same as that of the LTS estimator. As for the LTS estimator, there exists a stochastic algorithm, implemented in R, for approximating the LTA estimator; [START_REF] Giloni | Least trimmed squares regression, least median squares regression, and mathematical programming[END_REF] conducted a theoretical study of this estimator.

The least trimmed absolute deviations (LTA) estimator The LTA estimator is defined in the same way as the LTS estimator, but using the first form in (3), i.e it minimizes the least trimmed absolute deviations:

βLTA = Arg min β∈R d S LT A (r( β)),
where

S LT A (r) = h i=1 |r| (i) .
The breakdown point of the LTA estimator is the same as that of the LTS estimator. It has the advantage with respect to LTS that by adding (or splitting) the variables it can be written as a linear MIP. [START_REF] Hawkins | Applications and algorithms for least trimmed sum of absolute deviations regression[END_REF] proposed to solve it exactly by simple enumeration, using the following property of the L 1 minimizer: there always exists a subset of size d for which the fit is exact, i.e the residuals are 0. Hence it suffices to enumerate all subsets of size d, instead of the more numerous subsets of size h.

Robust estimators based on M -estimates of scale

M -scales: definition The M -estimate of scale, S M (r), is defined for each r as the solution of

1 n n i=1 ρ r i S M = b,
where the function ρ : R → R + satisfies the following hypotheses:

R1. ρ is even and twice continuously differentiable.

R2. ρ(0) = 0.

R3. ρ is strictly increasing on [0, c], for some 0 < c < +∞.

R4. ρ is constant on [c, ∞).
The constant b is conveniently defined as b = ρ(t)dΦ(t), where Φ denotes the standard normal distribution, in order to obtain consistency for the scale at the normal error model. τ -estimator of regression The τ -estimator of regression, βτ , is defined as

βτ = Arg min β∈R d S 2 M (r( β)) 1 nb 2 n i=1 ρ 2 r i ( β) S M (r( β)) , (6)
where S M (r) is an M -scale:

1 n n i=1 ρ 1 r i S M (r) = b 1 , (7)
ρ 1 and ρ 2 are functions satisfying R1-R4, and b 1 , b 2 are real constants.

The breakdown point of τ estimators equals * = b 1 /ρ 1 (c). Hence, for an adequate choice of the function ρ 1 the maximal breakdown point of 50% can be obtained. Similarly, by adjusting the function ρ 2 , the asymptotic efficiency at the normal distribution can be made arbitrarily close to 1. A related estimator is the so-called S-estimator, defined in the same way as τ -estimators, but replacing the objective function (6) just by β → S 2 M (r( β)). It has the same robustness properties of τ -estimators but a lower efficiency.

PRESENTATION OF THE STATISTICAL PROBLEM β 1 β 2 -4 -3 -2 -1 0 1 2 3 4 5 6 7 -4 -3 -2 -1 0 1 2 3 4 5 6 7 -3 -2 -1 0 1 2 3 4 5
The main difficulty in computing τ -estimators comes from its lack of convexity. In fact, in the context of estimators based on residual scales, convexity and robustness are antagonist concepts. For example, for the simpler S-estimators, the first-order optimality condition reads:

n i=1 ρ r i ( β) S M (r( β)) x i = 0
which is a weighted mean with weights w i = ρ (r i ( β)/S M (r( β))). Therefore, a convex ρ would give higher weights to observations with larger residuals. Note that the weights depend on the residuals, and so they are adaptive. This is a difference with estimators based on L-scales, where the weights have a rigid structure. This is an advantage from the viewpoint of statistics, but it is a drawback from an optimization point of view. Actually, we could not find any structure to exploit, contrarily to the LTS optimization problem, which is known to be a concave minimization problem (which is a type of "structured" nonconvexity). For this reason, we will investigate stochastic algorithms for τ -estimators, that adapt themselves to a wide class of problems. For the more structured LTS we will settle for warranted optimal solutions taking advantage of the structure of the problem.

A word on the multivariate case

For completeness, we will briefly describe here the multivariate context. We have a data set Z = {x 1 , ..., x n } ⊆ R d . We suppose that the n observations are independent realizations of a certain random variable X , which follows an elliptic distribution with induced measure absolutely continuous with respect to the Lebesgue measure, with a density h of the form

h(x) = det(Σ) -1/2 f ((x -µ)Σ -1 (x -µ)), (8)
where µ ∈ R d is called the location, Σ is a d × d positive definite matrix (henceforth Σ 0) called scatter matrix, and f is an arbitrary positive function. Two random variables with the same scatter matrix, up to some multiplicative factor, are said to have the same shape but different size. The expression in the argument of f in ( 8) is called Mahalanobis distance, it is the squared norm of the vector x -µ, for the norm induced by the inner product

•, • Σ -1 = •, Σ -1 • .
In the sequel we shall use the following notation for the Mahalanobis distance:

d 2 (x, µ; Σ) (x -µ) Σ -1 (x -µ).
Our objective will be to estimate the location µ and the scatter matrix Σ from the data set Z, which may be strongly contaminated. We shall search for estimators satisfying the same equivariance properties as their regression counterparts,

μ(AZ + b) = Aμ(Z) + b, Σ(AZ + b) = A Σ(Z)A , for any invertible d × d matrix A and any b ∈ R d .
The maximum likelihood estimators, which are also the LS estimators, are the empirical mean for location,

μ(Z) = 1 n n i=1 x i ,
and the unbiased estimator

Σ(Z) = 1 n -1 n i=1 (x i -μ(Z))(x i -μ(Z)) ,
for the scatter matrix.

The breakdown point for multivariate location is defined in the same way as for regression, namely, the smallest fraction of observations that need to be replaced by arbitrary ones to render the bias unbounded. For the scatter matrix, a usual measure of bias is max(λ d ( Σ), λ 1 ( Σ)), where λ d ( Σ) and λ 1 ( Σ) denote the largest and the smallest eigenvalue of Σ, respectively. Robust estimators for multivariate location are defined similarly as for regression, minimizing a robust scale, but this time the residuals are replaced by the Mahalanobis distances. The τ -estimators for multivariate location and scatter are defined as minimizers of the optimization problem

min µ,Σ det(Σ) n i=1 ρ 1 (d 2 (x, µ; Σ)) d s.t n i=1 ρ 2 (d 2 (x, µ; Σ)) ≤ b, µ ∈ R d , Σ 0, symmetric,
where the functions ρ 1 , ρ 2 play a role analogous to the functions in ( 6) and ( 7), and are required to satisfy the same hypotheses. The S-estimator of multivariate location and scatter is obtained through the same problem but keeping only det(Σ) in the objective function. An estimator analogous to the LTS regression estimator exists for multivariate robust estimation, it is called the Minimum Covariance Determinant (MCD). For a given n/2 < h ≤ n, the MCD is defined as the mean and covariance matrix of a subsample of size h for which the determinant of the covariance matrix is minimum.

For S and τ -estimators there exist undocumented implementations of stochastic algorithms similar to those existing for regression. The literature for the MCD is more extensive; we can cite [START_REF] Rousseeuw | A fast algorithm for the minimum covariance determinant estimator[END_REF] for two-step stochastic algorithms based on subsampling and concentration steps. Recently, [START_REF] Schyns | Relaxmcd: Smooth optimisation for the minimum covariance determinant estimator[END_REF] presented an algorithm based on a new approach consisting of relaxing the variables of the combinatorial problem for dealing finally with a smooth problem in continuous variables. [START_REF] Agulló | Computación de estimadores con alto punto de ruptura[END_REF] proposed a branch-and-bound algorithm for the exactly computing the MCD estimator. Recently, [START_REF] Nguyen | Outlier detection and robust covariance estimation using mathematical programming[END_REF] proposed a method for robust covariance estimation. It is not a high breakdown point estimator, but it is defined through a positive semidefinite program and therefore it can be efficiently computed even for large data sets.

Part I

Stochastic algorithms for approximately solving continuous global optimization problems.

Chapter I.1

Introduction

The first part of this thesis is devoted to τ -estimators for robust regression. This choice was done because the τ -estimator enjoys the best statistical properties: the robustness and the efficiency of the τ -estimator can be tuned simultaneously without trade-off. The computation of the τ -estimator of regression involves a non-convex differentiable optimization problem in continuous variables. One of its major difficulties, besides non-convexity, is that the feasible domain of optimization is the whole Euclidean space. Therefore, exhaustive exploration of the feasible domain is impossible; this rules out most deterministic algorithms for global optimization in continuous variables. At that point stochastic algorithms are one alternative, but the unboundedness of the domain appears as an issue again. However, this inconvenient has already been overcome in existing algorithms for robust regression. Most of them are based on a sampling technique called subsampling. Originally conceived for the LTS estimator, subsampling consists of drawing from a multinomial distribution a subsample of h observations, with d + 1 ≤ h ≤ n, from which a candidate β ∈ R d can be obtained as the LS fit to that subsample. Even if for the LTS optimization problem (4) subsampling boils down to uniformly sampling on the feasible domain {w ∈ [0, 1] n , n i=1 w i = h}, the rationale behind subsampling, namely that by drawing enough subsamples we should have the chance to pick at least one outlier-free subsample, remains appealing for any stochastic algorithm intended to compute robust estimators.

Sampling is the cornerstone of any stochastic optimization algorithm. Once a way of sampling over the region of interest is available, a range of alternatives show up, each one more or less adapted to some specific application.

Most authors classify stochastic global optimization algorithms into (see [START_REF] Schoen | Stochastic techniques for global optimization: a survey of recent advances[END_REF][START_REF] Törn | Global optimization[END_REF], for good surveys on the subject):

• Two-phase algorithms.

• Heuristic algorithms.

• Algorithms based on a stochastic model of the function.

The first two families of algorithms form a subclass to which we will further restrict ourselves, because they permit to take advantage of the differentiability of the objective function to perform local searches. The reader interested in the last family of algorithms is refered to [START_REF] Betrò | Bayesian methods in global optimization[END_REF] for a detailed review.

Let us briefly describe two heuristic algorithms that have shown the best performance in practice: simulated annealing (SA) and tabu search (TS). The basic idea of SA is to move through the search space by randomly choosing a neighbor of the current point that decreases the value of the objective function. However, to avoid being trapped in a local minimum, it may randomly allow a move that increases the objective function. A move from the current point to a new one is accepted with probability min{1, exp(-∆/T )}, where ∆ is the potential increase on the objective function that would result from this move, and T > 0 is a tuning constant called temperature. Note that when T → 0 the probability of accepting a move that does not improve the objective function decreases rapidly. The algorithm starts with a relatively large value of T which is then gradually reduced during the search. Eventually the algorithm freezes in a local minimum because no uphill moves can be accepted. Among all points visited by the algorithm the one with the smallest objective function is chosen as the approximate minimizer. Even if it was originally inspired from the physical process of driving a system to its minimum energy state by a "cooling" process, simulated annealing possesses a theoretical motivation. Let Ω denote the domain of interest. The Boltzmann distribution on Ω, π T , is defined as

π T (β) ∝ e -σ(β) T , β ∈ Ω (I.1.1)
where T > 0 is the temperature. It can be proved that π T converges, as T → 0, to the uniform distribution over the set of global minimizers of the function σ.

Since sampling from (I.1.1) would require the evaluation of σ on every point of Ω, the Metropolis algorithm can be used to obtain random samples from π T ; the Metropolis algorithm consists in simulating a sample path of a Markov process with transition kernel given by

K T (x, y) = e [σ(y)-σ(x)] + T R(x, y)
for any symmetric irreducible matrix R. Unlike SA, TS moves deterministically to the point that has the lowest objective value among all the neighbors. Once the search arrives at a local minimum in this way, the algorithm has to make an uphill move (the algorithm is forced to make a move even if every possible neighbor would yield an increase in objective value). In order to avoid that the algorithm immediately moves back to the local minimum in the next step, and thus become trapped, TS uses a tabu list that forbids moves that could yield to previously visited local minima, even if they would produce a better objective value.

On the other hand, there are two-phase algorithms. All of the algorithms in this class of algorithms proceed in three stages: 1. Sampling.

2. Local optimization.

Check stopping condition.

In [START_REF] Flores | On the efficient computation of robust regression estimators[END_REF], we have shown that most of the existing algorithms for robust regression fit in this class of algorithms. Furthermore, the computational study in [START_REF] Salibian-Barrera | The fast-τ estimator for regression[END_REF] shows that for optimization problems associated to robust regression, two-phase methods outperform heuristic algorithms. Nevertheless, stopping conditions are currently absent in existing algorithms. They also pass by clustering techniques that in some problems have helped to improve efficiency of two-phase algorithms. The impact of stopping conditions and clustering techniques for the approximation of robust regression estimators is the subject of Chapter I.2. This work is published in [START_REF] Flores | On the efficient computation of robust regression estimators[END_REF].

The main conclusion of Chapter I.2 is that stopping conditions are a valuable improvement to existing algorithms for robust regression. On the contrary, the impact of performing clustering is rather disappointing, especially in middling and high dimension. In Chapter I.3, one possible reason of the performance degradation of clustering global optimization in higher dimension is studied: the counter-intuitive behavior of the nearest neighbor in high dimension.

Chapter I.2

On the efficient computation of robust regression estimators S. Flores. Computational Statistics & Data Analysis 54 (12), 3044-3056.

I.2.1 Introduction

Robust regression methods have been introduced to cope with the need left by classical techniques for methods that work well in the presence of contamination in the data. Particular interest has been focused on estimators with a high breakdown point as defined by [START_REF] Donoho | The notion of breakdown point[END_REF]. Most of these methods combine robustness with desirable statistical properties such as consistency and asymptotic normality.Nonetheless, they are defined by means of difficult global minimization problems; hence, their computation is very time consuming. We address the problem of the efficient computation of robust estimators for statistical regression based on M-scales. The principal aim of our work is to investigate to what extent the most recent developments in the field of optimization can help improve the existing computational methods.

Let us consider the classical linear regression model

y i = x i β + ε i , i = 1, ..., n,
with errors terms ε i identically distributed with zero center and independent from the covariates x i . We shall denote by y the vector with components y i , by X the n × d matrix with ith row x i , and by r the vector of residuals r(β) = y -Xβ with components r i = y i -x i β.

Many robust estimators of the regression coefficients β ∈ R d based on n independent observations (y i , x i ) ∈ R × R d can be defined as:

β = Arg min β∈R d σ(r(β)), (P)
where σ is a scale estimator. An important case of (P) is the S-estimator, defined with σ(r) = s(r), where s : R n → R + is an M-estimator of scale, or M-scale (Huber, CHAPTER I.2. ON THE COMPUTATION OF ROBUST ESTIMATORS 1981), defined implicitly through:

1 n n i=1 ρ r i s(r) = b. (I.2.1)
The function ρ : R → R + is required to satisfy the following assumptions:

1. ρ is even and twice continuously differentiable, 2. ρ(0) = 0, 3. ρ is strictly increasing on [0, c], for some 0 < c < +∞,

4. ρ is constant on [c, ∞).
Any function satisfying these assumptions will be called in the sequel a "rho function". The constant b is conveniently defined as

b = E Φ (ρ), (I.2.2)
where Φ denotes the standard normal distribution, in order to obtain consistency for the scale at the normal error model. The S-estimator has very good robustness properties, but it lacks efficiency. Robustness is to be understood in the sense of the breakdown point, which is the minimum fraction of the observations that need to be shifted for the estimator to take on arbitrary values. In fact, there is a tradeoff between robustness and efficiency, and for a breakdown point of 50%, efficiency can be as low as 28.7% (Maronna et al., 2006, pp. 131). M M -estimators were introduced to fill this gap. They are obtained by local minimization of a suitable function using an S-estimator as starting point. In this way, they can combine efficiency with a high breakdown point. See Maronna et al. (2006, Sec. 5.5) for details. However, τ -estimators have lower bias curves, and their computing effort is comparable to computing the S-estimator, which is instrumental in the computation of M M -estimators.

In this paper, we shall focus on τ -estimators, introduced in Yohai and Zamar (1988), which are defined as minimizers of the function

σ(β) = s(r(β)) 2 nb 2 n i=1 ρ 2 r i (β) s(r(β)) , β ∈ R d , (I.2.3)
where ρ 2 is a rho function and b 2 is adjusted to ρ 2 as in (I.2.2). This choice is motivated by the robustness and efficiency properties of this estimator. Indeed, τ estimators have a breakdown point * = b/ρ 1 (c); hence, for an adequate choice of the function ρ 1 , the maximal breakdown point of 50% can be obtained. Similarly, by adjusting the function ρ 2 , the asymptotic efficiency at the normal distribution can be made arbitrarily close to 1. However, computing τ -estimators involves solving a difficult global optimization problem. Roughly speaking, global optimization methods can be classified [START_REF] Archetti | A survey on the global optimization problem: General theory and computational approaches[END_REF] into deterministic methods and stochastic or probabilistic methods. Deterministic methods look for a guaranteed global optimum, while stochastic methods settle for a point that is a global optimum within an allowed margin or with a certain probability. More ambitious, deterministic methods are time consuming, and their range of applicability is limited to very specific classes of problems, or to problems that are of small size. The interested reader should refer to [START_REF] Agulló | New algorithms for computing the least trimmed squares regression estimator[END_REF] for a deterministic robust regression algorithm. Most of the existing methods for computing robust estimators are stochastic and, more specifically, based on random subsampling. These methods operate by computing candidates βs based on subsamples of the observations and then starting local minimizations from each of these candidates. Therefore, they are also called multistart methods. Section I.2.2 is devoted to the local minimization aspects of computing τ -estimators. We should stress that this is the only part of our paper that is estimator-specific. The global part of our discussion is relevant to any objective function σ, provided that a way to perform local minimizations is available. Then, in Section I.2.3, we will briefly describe clustering global optimization methods, which is a class of multistart methods that uses clustering analysis techniques [START_REF] Törn | Global optimization[END_REF] in order to reduce the number of local minimizations needed to find a global minimum. Section I.2.4 discusses stopping conditions for multistart methods. Section I.2.5 shows how the existing methods for robust regression fit into the framework of clustering global optimization. In Section I.2.6, we present a few numerical tests comparing the different methods. In particular, the effectiveness of clustering techniques and stopping conditions is evaluated. We finish with our conclusions in Section I.2.7.

I.2.2 Local minimization issues

As already mentioned in the introduction, we focus on the problem of finding global minima of the function

σ(β) = s(r(β)) 2 nb 2 n i=1 ρ 2 r i (β) s(r(β)) , β ∈ R d ,
where the robust scale s(r(β)) is implicitly defined by:

1 n n i=1 ρ 1 (r i (β)/s(r(β)) = b 1 .
All the global optimization methods that we consider in this paper rely on local minimizations, and moreover, a major part of their computing time is spent in local minimizations. This is why fast and reliable local minimization algorithms are crucial. When the Hessian of the objective function is available, the most efficient local minimization algorithm is the Newton-Raphson method. Nevertheless, due to the flat parts present in the function ρ, the Hessian of the τ -objective function may contain large portions filled with zeros, and therefore, it is ill-conditioned. A good alternative for computing local minima of (I.2.3) is the Iterated Reweighted Least Squares (IRLS) algorithm [START_REF] Salibian-Barrera | The fast-τ estimator for regression[END_REF]. Although it has a slower rate of convergence, in practice it has proven to be quite efficient and stable. At each iteration, the IRLS algorithm solves a weighted least squares problem, which is equivalent to minimizing a quadratic local approximation of the objective function. In this section, we propose to use inexact solutions of the weighted least squares problems at each iteration, and we evaluate the gain in efficiency.

The IRLS step is derived from the necessary condition for optimality:

g τ (β) := ∂ σ ∂β = 0.
An expression for g τ (β) has been obtained in [START_REF] Yohai | High breakdown-point estimates of regression by means of the minimization of an efficient scale[END_REF]:

g τ (β) = -2 n X W (β)r(β).
Here, W (β) denotes the diagonal matrix with entries w j (β), where

w j (β) = ω(β)ρ 1 (e j (β)) + ρ 2 (e j (β)) e j (β) , (I.2.4)
with the following notations:

e i (β) = r i (β) s(β) , ω(β) = n i=1 [2ρ 2 (e i (β)) -ρ 2 (e i (β))e i (β)] n i=1 ρ 1 (e i (β))e i (β)
.

This leads to the matrix form of the optimality condition:

X W (β)Xβ = X W (β)y. (I.2.5)
Disregarding the fact that the matrix W depends (nonlinearly) upon β, the system of equations (I.2.5) are the normal equations associated to the Weighted Least Squares problem with weights W . It is a fixed point equation, so the iterative method

β k+1 = (X W (β k )X) -1 X W (β k )y (I.2.6)
has been proposed to solve it in Salibian- [START_REF] Salibian-Barrera | The fast-τ estimator for regression[END_REF]. These are the IRLS iterations.

For each iteration, IRLS constructs a quadratic approximation of the true objective function, where the minimum of this quadratic approximation is the subsequent iterate. The computational cost of each iteration is equivalent to the cost of computing the weights (I.2.4) and solving the system (I.2.6). Computing the weights (I.2.4) is costly because this requires the computation of an M-estimator of scale. Therefore, [START_REF] Salibian-Barrera | The fast-τ estimator for regression[END_REF] have proposed to replace the M-scale in (I.2.4) with an approximate value whereby the resulting iterations are known as approximated IRLS iterations. Approximated IRLS solves the d × d linear system (I.2.6) for each iteration, which corresponds to finding the minimum of a quadratic local approximation of the true objective function. However, when d is not small, the computational burden of solving a linear system for each iteration can be nonnegligible. Keeping in mind that our objective is to solve (I.2.5), an approximate solution of the instrumental subproblem (I.2.6) should be enough. In fact, replacing y = r(β) + Xβ and setting

B k = 2X W (β k )X/n in (I.2.6), we obtain -B k β k+1 = -B k β k + g τ (β k ).
In this form, the IRLS resembles the so called quasi-Newton methods, whose iterations are of the form [START_REF] Nocedal | Numerical optimization[END_REF]

-B k β k+1 = -B k β k +α k g τ (β k ), for an adequate steplength α k > 0.
It has been proved (Nocedal and Wright, 2006, Subsec. 5.7.1) that for quasi-Newton methods, a "good enough" direction suffices to keep convergence.

In order to evaluate the efficacy of this approach, we compared the computing time required to perform local minimizations by solving approximately (I.2.6), which will henceforth be denoted as the Iterated Inexact Reweighted Least Squares (IIRLS), with the computing time of the approximated IRLS iterations.

In our tests, we used the LSQR algorithm [START_REF] Björck | Numerical methods for least squares problems[END_REF] to calculate approximated solutions to least squares problems. This algorithm can handle problems with non-square matrices, and it is stable and easily available on the Matlab environment.

In Figure I.2.1, we show the results obtained from the IIRLS. The abbreviation IIRLS5 (respectively IIRLS20) stands for the algorithm performing 5 (respectively 20) iterations of the LSQR algorithm for each approximated IRLS iteration. We denoted IARLS as the approximated IRLS algorithm introduced in Salibian- [START_REF] Salibian-Barrera | The fast-τ estimator for regression[END_REF] that uses approximated M-scales for computing the weights. This is subsequently used to solve (I.2.6) using a direct method.

We plot, for different values of d from 5 to 200, the average time spent by each algorithm over 500 local minimizations of the function (I.2.3) with different datasets and different starting points generated as described in Section I.2.6. In all the cases, the IIRLS algorithm found the same local minimum as IARLS, although it usually needed more iterations to converge. Nevertheless, the economy in time of performing inexact iterations compensated for the increase in the number of iterations.

We see in Figure I.2.1 that for d = 5, there is not a great difference among the three methods, but the difference increases with d and becomes of great importance for d in the medium and large range. Considering the fact that the quality of the results is exactly the same, it is worthwhile to use IIRLS for local minimizations.

I.2.3 Clustering methods

This section describes a technique for solving global optimization problems, such as (P), by incorporating clustering analysis techniques. The objective of clustering methods in global optimization [START_REF] Törn | Global optimization[END_REF] is to identify groups of βs such that, if used as an initial point in a local minimization, every member of the same group yields the same local minimum. Thus, it would suffice to perform only one local minimization from each of these groups in order to locate all local minima, and the best of these minima would be a global minimum. A schematic explanation of clustering methods is presented in Figure I.2.2.

It consists in repeating the following steps iteratively, which we shall describe in detail in the rest of the section.

1. Sampling: sample candidates βs. Add them to the candidates sampled in previous iterations.

2. Concentration and/or selection: concentrate the sampled candidates around the minima in order to facilitate the clustering. . . . Let us describe in detail the first three steps: sampling, concentration/selection and clustering. For clarity of the exposition, stopping conditions will be discussed separately in the next section.

β sN → → → . . . → β c 1 β c 2 β c 3 . . . β c sN β c 1 β c 2 . . . β c γsN → → β 1 . . . β 2 . . . β w . . . (sN,w) ----→ Stop Cond

I.2.3.1 Sampling

In absence of additional information, candidates are usually sampled uniformly in the feasible region. However, when the feasible region is unbounded, it is not completely clear how the sampling should be done.

For the particular case of robust regression, [START_REF] Ruppert | Computing s estimators for regression and multivariate location/dispersion[END_REF] proposed a method known as random subsampling: candidates β i , i = 1, ..., N, are constructed by drawing random subsamples of size h ≥ d from the data and letting β i be the least squares fit to the ith subsample. The rationale behind this method is that for a large enough N , at least one outlier-free subsample could be sampled, which should give a good candidate, hopefully in the neighborhood of a global minimum.

I.2.3.2 Concentration and/or selection

The concentration step consists of performing some iterations of a local minimization procedure, usually one, starting from each candidate. In the selection step, a prespecified fraction of the sampled candidates with lowest function value is retained. It has been proposed [START_REF] Törn | Global optimization[END_REF] to do only the concentration step, only the selection step or both. For robust regression, [START_REF] Ruppert | Computing s estimators for regression and multivariate location/dispersion[END_REF] proposed to do selection and then concentration, while Salibian- [START_REF] Salibian-Barrera | The fast-τ estimator for regression[END_REF] explored the use of concentration followed by selection. The term "concentration step" has already been used in [START_REF] Rousseeuw | A fast algorithm for the minimum covariance determinant estimator[END_REF] to denote a particular local minimization procedure for the LTS estimator. In the sequel, we will use this term in the broader context described previously.

I.2.3.3 Clustering

Many ways to perform clustering for global optimization have been proposed [START_REF] Törn | Global optimization[END_REF]. Here, we review only Single Linkage, which is the simplest method (Rinnooy Kan and Timmer, 1987a).

Single linkage

At each iteration k, compute the radius

r k = 1 √ π Γ(1 + d/2) ξ ln(kN ) kN 1/d
, for some ξ > 0. Where Γ denotes the gamma function. The single linkage algorithm consists in iterating the following three steps until all points have been assigned to a cluster 1. Choose a local minimum to be used as seed.

2. Initialize the cluster with the seed.

3. Grow cluster : given a partially constructed cluster, iterate :

(a) find the unclustered point closest to the cluster.

(b) if this point is within distance r k from the cluster, add it to the cluster and repeat (3a). Otherwise, go back to step (1) and start the next cluster.

The theoretical convergence properties of Single Linkage have been proved to minimize a function over a bounded set S, supposing that the initial candidates have been sampled uniformly over S (and that no concentration step has been performed). The proof proceeds by estimating the probability that a local minimization is started from a point a at iteration k. This probability is bounded by the probability that there exists another candidate within distance r k with a lower function value. This is because if the ball with center a and radius r k contains a candidate z with lower function value, then if z is assigned to a cluster, a will be assigned to the same cluster. Moreover, if in step (1) the local minimizations are performed by first considering the candidates with lower function values, then we will not apply a local minimization to a before z is assigned to a cluster.

Therefore, the choice

r k = 1 √ π Γ 1 + d 2 ξVolume(S) ln(kN ) kN 1/d , ξ > 0,
makes the probability of applying a local minimization to any candidate decrease with k, for any 0 < η < 1/2, as O(k 1-ηξ ). See Rinnooy Kan and Timmer (1987a) for the details.

In the interest of improving the effectiveness of the clustering method, we take into account the following facts that could be detrimental to the clustering method:

In Kaufman and Rousseeuw (1990, Ch.5), various techniques of agglomerative clustering are examined. The authors alert about the chaining effect of single linkage, which makes the clusters stick to each other because of the formation of chains, and argue, based on theoretical analysis and practice, that some other techniques such as complete linkage or average linkage would be better suited for most problems. We have incorporated this insight into our numerical experiments, despite the fact that the previous theoretical analysis for adjusting the radius r k at each iteration does not carry over to the case where we consider the maximum distance to the cluster (complete linkage) or the average distance to the cluster (average linkage). In both cases the existence of a point within radius r with a lower function value does not ensure that a local minimization will not be started, unless there is a relationship between r and the (unknown) diameter of the cluster.

In [START_REF] Beyer | When is "nearest neighbor" meaningful?[END_REF], it had been pointed out that, for points sampled from a broad set of distributions, the distance of any of this points to its nearest neighbor becomes very close to the distance to the farthest point as the dimension increases. This essentially means that the notion of nearest neighbor loses much of its meaning in high dimension. Later on, in [START_REF] Aggarwal | On the surprising behavior of distance metrics in high dimensional space[END_REF], the role played in this phenomenon by the norm used to measure distances was devealed. It was shown that, in expectation as d → ∞, the gap goes to 0 for the ∞-norm, tends to a constant for the Euclidean norm, and goes to ∞ for the 1-norm.

In our numerical tests in Section I.2.6 we compare the actual impact of these two factors on the performance of clustering.

I.2.4 A stopping condition

A crucial issue in global optimization algorithms is the tradeoff between solution accuracy and computing time. In practice, this dilemma is the decision about when to stop.

When using random subsampling (Maronna et al., 2006, Sec 5.7.2), the minimization (P) over the whole R d is reduced to a search over a finite set of candidates generated from subsamples. For the probability of picking one outlier-free subsample from a sample with a fraction ε of contamination to be greater than 1 -δ, we need to sample a number N of candidates such that:

N ≥ |log(δ)| |log(1 -(1 -ε) d )| . (I.2.7)
Unfortunately, this approach has some drawbacks. First, the number N in (I.2.7) grows exponentially with d. For instance, if δ = 0.01 and ε = 0.25, we should sample 1450 candidates for d = 20, 25786 for d = 30 and more than 80 millions for d = 50. Furthermore, it depends on the fraction ε of contamination, which is not known in advance.

The main disadvantage of this criterion is its rigidity; it is an a priori criterion. Information about actually sampled candidates is completely disregarded. Thus, the algorithm will continue in the same way after 10 local minimizations as if it has found 8 local minima or only 1.

Adaptive criteria try to estimate the fraction of the search region that has been actually explored, using the observed information about the structure of each particular problem. Then, for a given level of accuracy, the running time of the algorithm will depend on the problem complexity, expressed mainly through the number of local minima that are found given a number of local minimizations.

In the sequel, we will describe an approach to this problem that uses Bayes' theorem to incorporate information gathered during the optimization process in order to decide when to stop. It was introduced in Boender and Rinnooy Kan (1987) and refined in [START_REF] Piccioni | Stopping rules for the multistart method when different local minima have different function values[END_REF]. The framework is the following: a sequential sample is drawn from a multinomial distribution with an unknown number of cells and unknown cell probabilities. In our context, each cell will be associated with one minimum of problem (P) and will be filled with the subsamples whose candidates converge after a local minimization to the minimum associated with this cell. The cell probabilities will be the fraction of subsamples in the cell.

Let us consider a small illustrative example. In Figure I.2.3, we elucidate a case with n = 20 and d = 4 where the objective function in (P) has 4 local minima β 1 , β 2 , β 3 and β 4 . Seven subsamples have been drawn, thus yielding by least squares

β 1 β 2 β 3 β 4 β 1 β 2 β 3 β 4 β 5 β 6 β 7 x 6 x 3 x 16 x 11 • • • x 12 x 20 x 2 x 14
x 15

x 14 x 5

x 20

x 8

x 9

x 17 x 14

x 19

x 12

x 11 x 16

x 3

x 7

x 1 x 16

x 12 x 6

x 20

x 8 seven candidates β 1 , ..., β 7 . If local minimizations were started from these candidates, β 3 would converge towards β 1 . β 2 and β 4 would converge to β 3 . Candidates β 1 , β 5 , β 6 and β 7 would converge to β 4 , and the minimum β 2 would remain undiscovered. Thus, we have observed 3 cells associated with the minima β 1 , β 3 and β 4 , with observed frequencies of 1, 2 and 4.

• • • ••• ••• ••• ••• ••• •••
In general, let β 1 , β 2 , ..., β k , be the local minima. Let ϑ i , i = 1, ..., k, denote the probabilities of each cell. In the Bayesian approach, the unknowns k, ϑ 1 , ..., ϑ k are supposed to be themselves random variables K, Θ 1 , ..., Θ k with realizations k, ϑ 1 , ..., ϑ k for which a priori distributions can be specified. Given the outcome (m 1 , ..., m w ) of a number of local searches, Bayes theorem is used to compute an a posteriori distribution. Following the approach of [START_REF] Piccioni | Stopping rules for the multistart method when different local minima have different function values[END_REF], we will suppose that different minima have different function values, σ(β i ) = σ(β j ) for i = j. In such case, the minima can be ordered according to their function values σ(β 1 ) < σ(β 2 ) < ... < σ(β k ). The same will be done with the minima found in experiments. Hence, we can compute the statistics O l = M J l , l = 1, ..., W m , where M i denotes the observed frequency of the ith minimum ( M 1 = 1, M 2 = 0, M 3 = 2 and M 4 = 4 in our example), J l the index of the lth observed minimum (in the example J 1 = 1, J 2 = 3 and J 3 = 4) and W m the number of (different) observed minima after sampling m candidates. In simple words, O l is the frequency of the lth observed minimum.

A particularly interesting quantity from the stopping criteria viewpoint is H = J 1 -1, which denotes the number of undiscovered local minima with better function value than the observed minima. Supposing a priori that the number of cells K follows an improper uniform discrete distribution on [1, ∞), and that given K = k, the cell probabilities Θ 1 , ..., Θ k are jointly uniformly distributed on the k -1 unit simplex, the following conditional probability for H can be obtained, when m ≥ 2 and w ≤ m -2 (Piccioni and Ramponi, 1990, Corollary 2):

P(H = h|W m = w, O l = m l , l=1,...,w ) = (m -w -1) (m -2)!(w + h -1)! (w -1)!(m + h -1)! .
In particular, the probability that the global minimum has been already discovered is

P m,w := P(H = 0|W m = w) = m -w -1 m -1 , w ≤ m -2. (I.2.8)
Using (I.2.8) we can readily devise a stopping condition, namely, to stop when the probability of having found the global minimum reaches a prespecified threshold.

The probability P m,w is undefined when w = m -1 or when w = m. However, in practice, this only occurs during the first two or three iterations. In this case, because we do not have evidence that the search region has been well explored, we should keep the algorithm running.

For example, let us consider the robust regression problem (P) for the stackloss dataset (Maronna et al., 2006, pp. 381). The Multistart method consist of sampling candidates βs and starting a local minimization from each of them, until the probability (I.2.8) reaches a given threshold. Let us consider the thresholds 0.3, 0.6 and 0.9. The first threshold 0.3 is reached after 4 local minimizations, which give 2 different local minima, none of which is the global one. The threshold 0.6 is reached after 9 local minimizations and 3 minima, one of which is the global minimum. Finally, the threshold 0.9 is reached after 42 local minimizations, producing 4 local minima.

Of course, since it is a random algorithm, its running is unlikely to be the same every time. What should be retained is that higher thresholds give more accurate results in the sense that the search region is more exhaustively explored, and this is done adaptively. Needless to say, a more exhaustive search will take longer than a rougher one. In general, there is not an easy way to guess how long will it take to solve problem (P) within a given accuracy, but one can always impose a time limit, and the value (I.2.8) can be given as information to the user at the end.

We would like to stress the fact that this approach works for any algorithm based on subsampling and local searches, even if the objective is to compute other estimators, such as LT S, or beyond the context of linear regression, such as the location and scatter estimation problem.

Stopping conditions based on (I.2.8) can also be used for algorithms described in Section I.2.3 that try to foresee the result of a local minimization, and to avoid it if it is likely to re-discover an existing minimum. In that case, in (I.2.8) m will be the number of sampled candidates and not the number of local searches actually carried out. The reason is that clustering methods are supposed to give the same outcome one could have obtained by starting local searches from each candidate. Nevertheless, the precision of the stopping condition will be subordinated to that of the clustering method; if it alters the outcome of the algorithm, the current a posteriori probability will be updated with wrong information.

I.2.5 A clustering global optimization point of view of robust regression estimation

We have already mentioned that most existing methods for computing robust regression estimators can be described as clustering methods, as described in Section I.2.3. In this section we shall describe them, and we will see how they perform each of the steps discussed in Section I.2.3.

I.2.5.1 Random resampling

The original version of random resampling was introduced by Rousseeuw in [START_REF] Rousseeuw | Least median of squares regression[END_REF] for computing the least median of squares estimator, and it was further refined and adapted for S-estimators by [START_REF] Ruppert | Computing s estimators for regression and multivariate location/dispersion[END_REF] • Sampling: Use Rousseeuw's random subsampling.

• Concentration-Selection: No concentration step is performed, but a selection of the t best candidates is done.

• Clustering: No clustering is performed. A local minimization is started from each candidate until convergence.

• Stopping condition: There is no stopping criterion. The number of sampled candidates is fixed in advance.

Selection IRLS

β 1 β 2 β 3 . . . β N t best β 1 . . . β t → . . . → β 1 . . .

I.2.5.2 Fast-S, Fast-τ

The "Fast-S" algorithm was presented in Salibian-Barrera and Yohai (2006) for Sestimators, and it was modified to cope with τ -estimators in Salibian- [START_REF] Salibian-Barrera | The fast-τ estimator for regression[END_REF]. This modification was called "Fast-τ ". From the global optimization viewpoint, it extended Random Resampling by adding a concentration step. As for random resampling, we give below a schematic illustration and a detailed description of each step. Given parameters N, ς and t, iterate the following steps:

• Sampling: Use Rousseeuw's random resampling.

• Concentration: Apply ς steps of IRLS to each candidate. Select the t candidates with best objective function.

• Clustering: No clustering is performed.

• Stopping condition: There is no stopping condition; the number of sampled candidates is fixed in advance.

IRLS

Selection IRLS I.2.5.3 Fast-τ with stopping condition.

β 1 β 2 β 3 . . . β N → → → . . . → β ς 1 β ς 2 β ς 3 . . . β ς N t best β ς 1 . . . β ς t → . . . → β 1 . . .
As previously observed, the considered algorithms for computing robust regression estimators are all particular instances of the general clustering procedure depicted in Figure I.2.2. The first one, random resampling, consists only of sampling and selection, Fast-τ adds a concentration step. Unlike clustering, the addition of an adequate stopping condition does not add computational work. For this reason, we introduce here a modification of Fast-τ for including stopping criteria, without doing clustering. It simply executes Fast-τ iteratively, and at the end of each iteration, it adds the (eventually new) minima found on the list of minima encountered in previous iterations, and evaluates a stopping condition. In our case, we require the probability P kN,w defined by (I.2.8) to exceed a given threshold θ.

Our proposition is summarized in Figure I.2.6. For given parameters N, ς and t, and a probability of success θ, Note that under this form, the modified version of Fast-τ does not exactly fit in the framework illustrated in Fig I .2.2 because the selection is performed considering only the candidates sampled in the last iteration and disregarding the ones from previous iterations.

IRLS

Selection IRLS Add New Min. Stop. Cond.

β 1 β 2 β 3 . . . β N → → → . . . → β ς 1 β ς 2 β ς 3 . . . β ς N t best β ς 1 . . . β ς t → . . . → β 1 . . . β t → β 1 β 2 . . . β w → if P sN,w ≥ θ stop Figure I.2.6:
The proposed modification of Fast-τ including a stopping condition.

I.2.6 Numerical tests

We conducted numerical experiments in order to investigate the impact of clustering techniques and stopping criteria. As in [START_REF] Ruppert | Computing s estimators for regression and multivariate location/dispersion[END_REF] and Salibian-Barrera et al. ( 2008), we consider simulated data contaminated with clustered outliers, where contamination is highly concentrated around outlying values.

• (1 -δ)100% of the points follow the regression model y = Xβ + ε , where the d -1 covariates are distributed as N d-1 (0, I d-1 ), x id = 1 is the intercept, ε i ∼ N (0, 1) and β = 0.

• δ100% of the points are "bad" high-leverage points: the covariates now follow a N d-1 (s, 0.1 2 I d-1 ) distribution where s = (100, 0, ..., 0) t , and the response variables y i ∼ N (100l, 0.1 2 ). In the following we will call the parameter l the "contamination slope".

For evaluating the stopping condition, we tested only Multistart, and we compare it with Fast-τ using N = 50 samples. For performance evaluations, we tested all the algorithms described in Section I.2.5, except for Random Resampling (Subsection I.2.5.1). The reason for excluding Random Resampling from those tests is that we included the Fast-τ algorithm, and the numerical tests reported in Salibian-Barrera et al. ( 2008) already compared Fast-τ with Random Resampling, among others, and showed that it outperforms Random Resampling.

The clustering methods were implemented by the author in Matlab and are available upon request. In our test, it was always executed with the same parameters: at each iteration, N = 100 candidates are sampled and added to the sample, then a concentration step consisting of one IRLS iteration (cf. Section I.2.2) and a selection of the best 10% of the concentrated candidates is performed (see Figure I.2.2). Concerning the radius used for the clustering, as the minimization in (P) is done over the unbounded domain R d , we used the formula r k = π -1/2 (Γ(1 + d/2)ξ ln(kN )/(kN )) 1/d , for a predefined ξ > 0. After extensive experiments, we realized that the results are rather insensitive to the choice of the parameter ξ. We set ξ = 40 in our tests. For the tests involving the Fast-τ algorithm, we used the code available from the webpage of Matias Salibian-Barrera. We used the parameters ς = 2 and t = 5 (see Figure I.2.5). The number N of initial candidates changed from test to test and is indicated each time.

Both Fast-τ and Single Linkage clustering were fitted with a stopping condition. In all cases, it consisted in stopping when the probability P m,w defined in (I.2.8) reached a given threshold θ.

Concerning the performance in terms of computing time, for each algorithm we saved the required time T A . Then, we computed the ratio RT = T A /T R , where T R is the time required by an algorithm used as reference. We usually used Fast-τ as reference, with different parameters depending on the test. This ratio is what we call the relative computing time. In this way, the results are machine-independent; and easier to compare.

As it is not feasible to know certainly if the returned solution is the global minimum, we only give the relative τ -scale with respect to a reference algorithm, σA /σ R , where σA and σR are the values of the objective (I.2.3). We call this ratio "relative τ -scale". Similarly to the results reported in Salibian-Barrera and Yohai (2006) for S-estimators, for the particular type of contamination that we are considering, often the coefficient β1 over the 500 samples forms two well-separated groups: one highly concentrated around the contamination slope, and another one more dispersed around 0, the slope of the data without contamination, hereafter referred to as "clean slope". In those cases, we also show the percentage of samples for which the "slope" is around 0. Note that it is not uncommon that minima with the contamination slope had better function values (I.2.3) than minima around the "clean" slope, especially for those with high proportions of outliers.

In Subsection I.2.6.1 we evaluate the impact of the stopping condition. We discuss in detail how it behaves when applied to Multistart with threshold values ranging in a wide range. Our second and third test, presented in Subsection I.2.6.2 and I.2.6.3 , evaluate the performance of the considered algorithms in two different kinds of situations. In I.2.6.2, this is done for a small problem, whose complexity varied only through the change in the dimension d. In the third test, presented in Subsection I.2.6.3, we compare the performance of the considered algorithms, the parameter used to control the complexity of the problems was the contamination slope l.

I.2.6.1 Analysis of the stopping condition

The objective of this Subsection is to scrutinize the impact of the stopping condition presented in Section I.2.4 on the algorithms in Section I.2.5. The datasets were generated as indicated at the beginning of Section I.2.6, with n = 400 observations in dimension d = 15. The contamination fraction was δ = 40% forming two groups of 20% each with contamination slopes l = 2 and l = 4. The objective of using two groups of outliers was to increase the number of local minima, which would better illustrate the effectiveness of the stopping condition, since the complexity of the problems heavily depends on the number of local minima. The algorithms compared were Multistart (MS), which consist of sampling candidates and launching a local minimization from each of them; it is the same as the Random Resampling algorithm (cf. Section I.2.5.1) without selection. and Fast-τ with N = 50 samples (FT50). Multistart stopped as soon as the probability P m,w in (I.2.8) reached a given threshold θ. In this section, we used the thresholds 0.3, 0.6, 0.9 and 0.95.

In Table I.2.1 we show the following:

• The number of local minima, w.

• The quantity of sampled candidates, m.

• The effective value of the probability (I.2.8) at the termination of the algorithm, P m,w .

• The percentage of samples for which the estimate given by the algorithm had the clean slope, around 0, at the row "% Slope ".

• The relative running time with respect to FT50, T M S /T F T 50 .

• The relative τ -scale with respect to FT50, σMS /σ F T 50 .

All the entries, except for the clean slope, are the average over 500 simulations. By examinating the column corresponding to FT50 in Table I.2.1 we see that, in this example with few local minima, most of the time the effective value of (I.2.8) is already around 0.96; if the function used to have 5 local minima, it would be around 0.9 and it would be around 0.8 if the function had 10 local minima. The reason is that candidates are sampled in batches; thus, if each batch is of size N , the only attainable threshold values are of the form (N -w -1)/(N -1), for positive integers w. Since in MS, candidates are sampled one by one, it stops as soon as the threshold is reached so it can be combined with low threshold values. On the contrary, algorithms that do some kind of selection need to sample in batches; for performing local searches, only from promising candidates. In the case of algorithms performing clustering, a harsh selection is needed in order to well separate clusters and prevent sticking to each other, thus the use of small batches is discouraged. As a consequence, for all of the algorithms of Section I.2.5, if the objective function does not have many local minima, only relatively high threshold values will be observed in practice.

However, MS illustrates quite well how the stopping condition works, since we see how the number of discovered local minima and the number of sampled candidates increases as the threshold increases. Even if the percentage of datasets for which the estimate had the clean slope seems to indicate the contrary, the accuracy of the result also improved. This can be seen by examining the relative τ -scale. We see that MS gives worse results than FT50 for low threshold values, but this ratio decreases for higher threshold values, showing the better performance of MS.

I.2.6.2 Many dimensions, many thresholds on a small problem

In our second test, we set the contamination fraction to δ = 0.2 and the contamination slope to l = 2. For each d ∈ {5, 10, 15, 20}, we generated 500 datasets of n = 100 points. We compare four algorithms:

• Fast-τ with N = 500 candidates (FT500).

• Fast-τ with N = 250 candidates (FT250).

• Fast-τ with a stopping condition (SC), described in subsection I.2.5.3. The parameters are (see Figure I.2.6) N = 100, ς = 1 and t = 5.

• Single Linkage clustering as described in Subsection I.2.3.3 (SL), with the parameters indicated at the beginning of this Section. Motivated by the observation at the end of Section I.2.4, we tested Single Linkage using the 1 and the Euclidean norm. They are denoted as SL1 and SL2, respectively.

For those algorithms incorporating a stopping condition (SC and SL), the threshold parameter θ took the values 0.95, 0.97 and 0.99. These results are shown in Table I.2.2 for the algorithms with stopping condition, and in table I.2.3 for Fast-τ with N = 250 and N = 1500 samples.

For algorithms with stopping condition, a curious situation appears, as the algorithms performing clustering always perform worse than those without clustering, but in low dimension, they find quite often a better solution than the optimal one, in the sense that they find a non-global minimum with the clean slope. As the dimension increases, however, their performance degrades both in terms of objective function and in the percentage of times that they find the clean slope. Single Linkage using the 1 norm generally had better function values than using the Euclidean norm, but the difference is negligible. In preliminary tests we also tried Complete Linkage and Average Linkage clustering, but the results were essentially identical, so we only used Single Linkage.

Overall, the Fast-τ algorithm with stopping condition using thresholds 0.95 and 0.97 achieves a good tradeoff between computing time and quality of the solution, as it gives results almost as good as FT500 or even FT1500, but within a fraction of time. By raising the threshold to 0.99, the results are similar to those of FT1500, but with a larger computing time.

I.2.6.3 Varying complexity for a fixed dimension

In the next test, we examine the behavior of the algorithms for a fixed dimension, but also for various proportions of outliers and contamination slopes. The closer the contamination slope is to the clean slope, the more difficult it becomes to identify the clean one. We tested FT250, FT500, FT1500, SL and SC in dimension d = 10 with different contamination slopes. The number of observations was fixed to n = 400, and the proportion of outliers was δ = 10%, 15% and 20%. For SL and SC, the stopping condition uses the threshold θ = 0.95; because the preliminary test did not show significative improvements when we raised the threshold in these problems. Similarly to the previous test, the use of different variants of clustering did not significantly change the results; therefore we only tested Single Linkage clustering, with the usual Euclidean norm.

The performance was measured as in the previous test. Namely, we compare the objective function value with respect to a reference algorithm, and we keep the percentage of the samples from which convergence occurred to a minimum with a slope around 0. These two quantities were the measure of quality of the solution. We also record the time relative to the time spent by FT500, as explained in the previous section.

The quality of the solutions obtained by Fast-τ was identical to SC; thus, we do not include it in the tables. The results of these test are in Table I As in the previous test, the algorithm with clustering has a very good running time, about 20% of the time spent by Fast-τ , and in some cases, it gives a worse solution to the minimization problem (P) than algorithms without clustering. Thus, it does not compute the τ -estimator, but the result it gives has the clean slope. We do not have a clear explanation for that, and we think this phenomenon bears closer examination; in another article. We see once again that the stopping condition permits to find the global minimum with a high probability in a proper computing time.

I.2.7 Conclusions and future work

We have investigated the effectiveness of the usage of clustering techniques and stopping conditions for global optimization in the particular case of robust regression. Our viewpoint is that of an user of robust regression who wants to compute robust estimators without having to adjust parameters that depend upon the details of the chosen algorithm.

The integration of a stopping condition is completely justified both by the quality of the results and by the performance in terms of computing time. Additionally, it is very simple to implement in new and existing software. It should be incorporated in algorithms not only for computing τ -estimators, but also in any algorithm based on subsampling and concentrations steps.

A threshold between 0.95 and 0.97 achieves a good compromise between efficiency and quality of the results. A higher threshold value ensures a very good solution; at an extra cost in terms of computing time. For routine utilization, a threshold of 0.96 or 0.97 should give good results at a competitive computing time, which will adapt by itself to the complexity of the problem. Additionally, it possesses a very appealing interpretation in terms of probability of finding the global optimum, independently of the particular problem under consideration.

In their present states, the existing clustering techniques for global optimization do not seem to fit the needs of robust estimation problems. However, as shown in Section I.2.5, they are a natural extension of the existing methods, and their computing times and behavior in some difficult problems suggest that they deserve further investigation.

Chapter I.3

High dimensional issues in clustering (global optimization)

I.3.1 Introduction

While evaluating the use of clustering techniques to improve existing global optimizations methods particularly suited for robust regression problems [START_REF] Flores | On the efficient computation of robust regression estimators[END_REF], we were very disappointed by the poor performance of some clustering techniques in dimension as moderate as 10 or 15. Trying to understand the breakdown of (single linkage) clustering, we learned that the same problem had been highlighted by computer science researchers working on high dimensional indexing [START_REF] Beyer | When is "nearest neighbor" meaningful?[END_REF][START_REF] Aggarwal | On the surprising behavior of distance metrics in high dimensional space[END_REF], and by people working with other clustering algorithms [START_REF] Schoen | Global optimization methods for high-dimensional problems[END_REF]. In fact, the common point between clustering and indexing is the nearest neighbor subproblem they rely on. Therefore, the problem is unlikely to be restricted to single linkage clustering as, quoting from Locatelli and Schoen (1999, pp. 380), « All of the papers dealing with clustering methods base this decision on some sort of nearest neighbor statistics -that is, the decision of starting a local search is in some way connected with the fact that no point in a properly defined subset of the sample is too near to the current one».

The phenomenon described in [START_REF] Beyer | When is "nearest neighbor" meaningful?[END_REF]; [START_REF] Aggarwal | On the surprising behavior of distance metrics in high dimensional space[END_REF] is: under rather mild hypothesis, the distance to the nearest point approaches the distance to the farthest point as the dimension increases. This essentially means that the notion of nearest neighbor loses much of its meaning in middling and high dimension.

The objective of this chapter is to go into the consequences of the blurring of the notion of nearest neighbor for clustering algorithms in high dimension, and to propose amendments permitting to uphold the efficacy of clustering algorithms as the dimension increases.

I.3.1.1 Clustering global optimization

We focus on the continuous global optimization problem,

min β∈Ω σ(β), (I.3.1)
where Ω ⊆ R d is a compact set, and σ : Ω → R is a continuous function. We will assume that a local minimization procedure is available.

Algorithm 1 Multilevel Single Linkage (MLSL) Choose N > 0, ξ > 0 and γ ∈ (0, 1]. Set k := 0, X = ∅ and iterate:

1. Let k := k + 1; 2. let r = r k,N,ξ = π -1/2 Γ 1 + d 2 vol(Ω)ξ log(kN ) kN 1/d ; (I.3.2)
3. generate a uniform random sample of size N in Ω, add it to X ; 4. sort X by increasing function value, and select the γkN best of them, call it X k ; 5. launch a local search algorithm from each point in X k except if

• the distance of the point from the boundary is less than a fixed threshold.

• the point is "too near" from a critical point Or,

• if there exists another point with better function value at a distance which is less than or equal to the radius r k ;

6. Check if some stopping condition is satisfied, otherwise repeat from step 1.

Rinnooy Kan and Timmer (1987a,b) introduced a stochastic global optimization method for solving (I.3.1) using clustering techniques. It is an improvement of the multistart (MS) method, which proceeds by sampling points uniformly on Ω, followed by a local search from each of the sampled points. The MS method has the evident disadvantage of potentially wasting time by launching several local searches yielding the same minimum. The method proposed in the cited articles, called multi-level single linkage (MLSL, Algorithm 1), aimed at overcoming this drawback by avoiding the local search if it is likely to find a minimum already known. MLSL enjoys the following asymptotic properties as the number of iterations goes to infinity: P1. The best observed value converges to the global optimum with probability 1.

P2. The probability of starting a local search decreases to 0.

P3. The algorithm performs a finite number of local searches with probability 1, even if it runs forever.

Despite its strong theoretical properties, the necessity of revising at each iteration the decisions about local searches and some other technical details left some room for improving MLSL. This was done by [START_REF] Locatelli | Simple linkage: Analysis of a threshold-accepting global optimization method[END_REF], where the authors propose a family of algorithms for global optimization satisfying properties P 1, P 2 and P 3 without the computational overhead of MLSL. In a later article [START_REF] Locatelli | Random Linkage: a family of acceptance/rejection algorithms for global optimisation[END_REF], a particular member of this family, called simple linkage (SL, Algorithm 2) was analyzed in detail and it was shown to make almost the same decisions as MLSL.

Algorithm 2 Simple linkage (SL)

Choose ξ > 0 and ε > 0. Set k = 1, generate β 1 uniformly and iterate:

1. Let k = k + 1; 2. let r k be defined as in (I.3.2), with N=1; 3. generate a single uniform random point

β k ; 4. launch a local search from β k except if ∃j < k, β j -β k 2 ≤ r k : σ(β j ) ≤ σ(β k ) + ε; (I.3.3)
5. check stopping condition, if it is not met, repeat from step 1.

In [START_REF] Locatelli | Random Linkage: a family of acceptance/rejection algorithms for global optimisation[END_REF] it was also proved that SL enjoys property P2 if and only if lim k→∞ k 1/d r k = ∞, and that if ξ > 2 d /2d then P3 holds as well.

A key difference between MLSL and SL is that MLSL never starts local searches from points within a prescribed distance from the boundary, while SL allows local searches to be started from any point in the feasible domain. This difference is particularly important in high dimension. For instance, given a threshold distance ν, the probability of sampling from a uniform distribution within distance ν from the boundary of the unit cube is 1 -(1 -2ν) d . Moreover, if local searches from points at a distance from the boundary less than a fixed threshold ν are inhibited in SL, then property P3 holds if ξ > 1. Even if the practical interest of property P3 is not clear, this is a first evidence that if radius r k,ξ is used, the parameter ξ should vary with the dimension d.

Shortly after, [START_REF] Schoen | Global optimization methods for high-dimensional problems[END_REF] pointed out some problems affecting clustering global optimization algorithms in high dimensions, and proposed some solutions. At that time there was already some evidence indicating that the use of formula (I.3.2) forced algorithms to behave much like best start (BS), a quite inefficient random sampling method that samples a point and launches a local search only if the new point has the overall best function value.

The radius in equation (I.3.2) has been derived in [START_REF] Boender | A stochastic method for global optimization[END_REF]; Rinnooy Kan and Timmer (1987a) by approximating the distribution of the nearest neighbor statistic within a set of uniformly distributed points, and has been used afterward in many clustering algorithms. In fact, r k was chosen in such a way that

π d/2 Γ(1 + d/2) r d k = vol(B r k ) = vol(Ω)ξ log(k) k ,
where B r denotes a ball of radius r for the 2 norm. The asymptotic development for the Γ function, as d → ∞ for fixed k, shows that [START_REF] Schoen | Global optimization methods for high-dimensional problems[END_REF])

r k ∼ ξ log(k) k 1/d (d/2) 1/2+2 -d .
That number grows thus as √ d for large d, as the diameter of the 2 unit ball does, forcing the radius (I.3.2) to be too big at the first iterations. For these reasons, [START_REF] Schoen | Global optimization methods for high-dimensional problems[END_REF] proposed to use the ∞ norm instead, for which the diameter of the unit ball is constant in any dimension. The resulting Algorithm, called ∞-Simple Linkage (∞-SL), follows the lines of Algorithm 2, using the ∞ norm in (I.3.3), and replacing radius r k,ξ by

s k = s k,ξ = (ξ log(k)/k) 1/d . (I.3.4)
Algorithm ∞-SL satisfies property P3 for any ξ > 1.

Nevertheless, [START_REF] Beyer | When is "nearest neighbor" meaningful?[END_REF] pointed out that the difference between the distance from any point to its nearest neighbor and the distance from the same point to the farthest point in the sample, hereafter simply called the gap, shrinks to 0 as the dimension grows. Later on, in [START_REF] Aggarwal | On the surprising behavior of distance metrics in high dimensional space[END_REF], the role played in this phenomenon by the norm used to measure distances was revealed. It was shown, for a wide range of distributions, that in expectation as d → ∞, the gap goes to 0 for the ∞ norm, tends to a constant for the 2 norm, and goes to ∞ for the 1 norm (see the left part of Figures I.3.4, I.3.3 and I.3.2).

In the next sections, using the results of [START_REF] Beyer | When is "nearest neighbor" meaningful?[END_REF] and [START_REF] Aggarwal | On the surprising behavior of distance metrics in high dimensional space[END_REF] on the behavior of the nearest neighbor in high dimensional spaces, we show that the problem of choosing the radius is indeed very delicate. We propose an alternative method intended to be independent of the dimension, the feasible domain and the sampling distribution.

I.3.2 Nearest neighbor behavior in high dimension

In this section we shall show how the poor performance of clustering algorithms in high dimension, particularly their excessive resemblance to MS or BS, is related to the above-mentioned "gap phenomenon".

Let us focus on condition (I.3.3) at step 3. of Algorithm 2 (the following also holds for MLSL). Denote by D k min and D k max respectively the minimum and maximum distance from the last sampled point β k to the other points β j , j = 1, ..., k -1 in the sample. Three scenarios are possible at this iteration: i) r k < D k min . A local search will be started, independently of the function values of the rest of the sample. At this iteration the algorithm will make the same decision as Multistart, so we call this iteration a 'Multistart iteration'. iii) D k max ≤ r k . A local search will be started only if β k has the best overall function value. Thus, the algorithm will make the same decision as that of best start, and we shall talk of a 'best start iteration'.

We see that for small r, we are always in case (i), and the algorithm behaves like multistart. On the contrary, for large r, the function value of all the points in the sample are taken into account before performing a local search, hence the algorithm behaves like best start. It is precisely for r k in the range between D k min and D k max that more complex algorithms can outperform basic methods like multistart and best start. In Figure I.3.1, we illustrate the situation with Ω = [0, 1] 35 . During 500 iterations, we show the maximum (top), minimum (bottom) and the median (between) of the distances from the points sampled at previous iterations to the last sampled point. For a given curve r k , at the ks for which the curve lies on the upper region, the algorithm will behave like BS, and at those k for which the curve lies on the lower region, the algorithm will make the same decision as MS. It results clear that if the space between the top and the bottom curves wipes out and the choice of the radius does not follows accordingly, the algorithm will behave either as MS or BS.

Going into the details, Aggarwal et al. (2001, Lemma 2) show, for a large set of distributions, that for a sample of size N ,

C(p) ≤ lim d→∞ E D max -D min d 1/p-1/2 ≤ (N -1)C(p), (I.3.5)
where C(p) is a constant depending on p, the distances being measured according to the p-norm, and D max (resp. D min ) denotes the maximum (resp. minimum) distance from any point to the other points in the sample. Note that D k min and D k max depend also on the dimension and the norm. A similar result holds for "fractional norms", i.e for the maps (x 1 , ..., x d ) → (x p 1 + • • • + x p d ) 1/p with p ∈ (0, 1). The result (I.3.5) announces that:

• for p > 2, the gap tends to 0 as the dimension d increases, promoting subite switching from best start regime to a multistart regime;

• for p = 2, the gap tends, for a fixed number of sampled points, to an unknown constant;

• for 1 ≤ p < 2, the gap increases to ∞ as the dimension raises.

As in (I.3.5) the limit as d → ∞ is taken for a constant sample size, and in SL-type optimization algorithms the current sample at iteration k consists of k points, we simulate the stage-wise sampling of algorithm 2 in dimensions d = 5, 10, 15, 20, 25, ..., 105. For k = 2, ..., 500, we sample uniformly a point β k in [0, 1] d , then we compute the distance for the p-norm, β k -β j p for j = 1, ..., k -1 and save the α-quantile q k α such that

α = P( β k -β j p ≤ q k α )
for α = 0, 0.25, 0.5, 0.75 and 1. In particular, the 1-quantile q k 1 coincides with D k max , and the 0-quantile q k 0 coincides with D k min . For fixed d and α, most of the plots of the distances against the iteration number k look much as Figure I.3.1. So a first observation is that the influence of changing number of points k is not significant. Much more interesting are the plots of the gap q 1 -q 0 and the inter-quartile difference q .75 -q .25 for varying d and fixed k, shown in Figures I.3.4, I.3.3 and I.3.2 for the ∞ , 2 and 1 norm, respectively.

As predicted by the theoretical results, for p = ∞ (Figure I.3.2) the gap between the curves depicting the closest point and the farthest point to β k vanishes as the dimension increases. The augmentation of the number of points does not seem to be significant enough to change this trend. For p = 2 (Figure I.3.3), the variability due to the change in dimension is not significant. The situation is better when using the 1 norm, as shown in Figure I.3.4. As predicted, the gap increases with dimension. We can also observe that even in the lowest considered dimension, d = 5, the gap is of the same magnitude as for the 2 or ∞ norms. We conclude from these observation that the 1 norm is the best alternative for distinguishing points in any dimension. For p < 1, even if the maps (x 1 , ..., x d ) → (x p 1 +• • •+x p d ) 1/p are not norms, theoretical results say that the gap should increase faster than with p = 1, and this is actually the case. Nevertheless, besides the embarrassment at interpreting the outcome of functions that are not norms, they present the inconvenient that in low dimension the gap is quite close to 0. The gap q 1 -q 0 (left) and the inter-quartile difference q .75 -q .25 (right) for the ∞ norm. The gap q 1 -q 0 (left) and the inter-quartile difference q .75 -q .25 (right) for the 2 norm. The gap q 1 -q 0 (left) and the inter-quartile difference q .75 -q .25 (right) for the 1 norm. in (I.3.4) with ξ = 2 (dotted line), and the ∞ distances from β k to its nearest point in the sample, the farthest point in the sample and the median of the ∞ distances from β k to the previously sampled points, as in Figure I.3.1. We see how, as the dimension increases, both the curve s k and the curve of the median go up, but the narrowing of the gap makes s k to get closer and closer to the curve of the farthest distances, indicating that the algorithm behaves like BS with respect to deciding to start (or not) a local search.

In

However, the numerical results reported in [START_REF] Schoen | Global optimization methods for high-dimensional problems[END_REF] seemed rather good. This is probably because of another modification they proposed: make the radius depend on the sampled point. Again, this modification was introduced through the parameter ξ, which is no longer fixed in advance as in Algorithm 2, but is computed once β k has been sampled, and depends both on d and β k : ξ = ξ d (β k ).

I.3.3 Replacing radius by quantiles

The direct application of the nearest neighbor distribution to estimate the clustering radius does not give good results during the first iterations. Moreover, finding a useful functional form for r = r(k, d, ξ, β k , p) for more general domains than the unit cube seems to be out of reach, even for uniformly distributed points.

A dual viewpoint is that of quantiles. At the beginning of iteration k, k -1 points β 1 , ..., β k-1 have already been sampled. Then, we sample β k and every point from β k . Hence, by taking a radius r ≥ r 1 we are considering all points, in other words, we are considering the 1-quantile. On the other extreme, by taking r ≤ r 0 the local search is launched without considering any function value, or equivalently, taking the 0-quantile for the threshold radius.

β i , i = 1, .., k -1 is within distance r 0 = min i β k -β i p and r 1 = max i β k -β i p
In the same way, for intermediate radius r 0 < r < r 1 we are considering the α r -quantile, with α r = P( β k -• p ≤ r | β 1 , ..., β k-1 ). Conversely, if we consider the α-quantile, 0 < α < 1, we can associate to this quantile a radius r α obtained as the distance from β k to the farthest point within this quantile. Note that for α = 0 or α = 1 we cannot associate a unique radius, but only the intervals [0, r 0 ] and [r 1 , ∞] respectively.

Algorithm 3 Quantile Global Optimization (QGO) Choose a norm p, a sampling distribution F and a sequence α k . Set k=1, randomly generate β 1 ∼ F and iterate:

1. Let k = k + 1; 2. randomly generate β k ∼ F and compute the p -distance to the previously sampled points; 3. let r k be the α k -quantile of the distances computed in 2; 4. launch a local search from

β k except if ∃j < k, β k -β j p ≤ r k and σ(β j ) ≤ σ(β k ); (I.3.6)
5. check stopping condition, if it is not met, repeat from step 1.

Based on the later observation, we propose and evaluate Algorithm 3, called quantile global optimization (QGO). It is formulated for an arbitrary norm p and sampling distribution F , and depends on the sequence {α k } of quantile orders.

Some key aspects should be taken into account for choosing the sequence α k :

• At the first iteration, the probability of finding a new minima is 1, hence we perform an MS iteration.

• During the fist iterations, the probability of finding new minima is largest, so we would like to be 'close to MS'.

• As the iterations go, the probability of rediscovering known minima raises, then we would like to start leaving from MS and getting closer to BS.

These guidelines suggest taking increasing sequences α k such that α 0 = 0 and α k → 1 as k → ∞. For such a quantile order sequence, the following holds:

Theorem 1 Consider the QGO algorithm described hereabove, with 1 ≤ p ≤ ∞ and a sampling distribution F whose support contains the whole feasible domain Ω. Then the following properties hold, as k → ∞:

1. The best observed function value converges to the minimum of the function σ over Ω with probability 1.

2. If the sequence α k converges to 1, the probability of starting a local search converges to 0.

Proof.

The first affirmation follows from the fact that since Ω ⊆ support(F ), then the random sampling tends to cover the whole domain. Let us denote as σk the best observed value at iteration k. In order to prove 2, let us rewrite step 3 of Algorithm 3 as: launch a local search if

∀j < k, σ(β j ) > σ(β k ) or β j -β k p > r k ,
then, if we denote by E A k the event «Algorithm A starts a local search at iteration k », we have

P(E QGO k ) = P j<k {σ(β j ) > σ(β k )} j<k β j -β k p > r k (I.3.7) ≤ P(σ(β j ) > σ(β k ), ∀j < k) + P( β j -β k p > r k , ∀j < k) = P(σ(β j ) > σ(β k ), ∀j < k) + (1 -α k ). Next we rewrite {σ(β j ) > σ(β k ), ∀j < k} as {σ k > σ(β k )} to conclude that P(σ(β j ) > σ(β k ), ∀j < k) = P({σ k > σ(β k )}) → 0, because σk → inf Ω σ.
Since 1 -α k → 0 as well, the second assertion follows.

In order to illustrate how the QGO algorithm works, in Figure I.3.6 we display an hyperbolic tangent-shaped sequence α k → 1, and the corresponding radii obtained using quantiles in the same situation as that of Figure I.3.1. We can observe that at the beginning the algorithm behaves like MS, and tends asymptotically to behave like BS, and this, independently of the dimension, the norm used to measure distances, or the domain Ω.

A notable such sequence is α k = 1-p k , where p k denotes the probability of finding a new minimum. This proposal follows closely the previous reasoning, since at the first iteration the probability of finding a new minimum is 1, because no local minima are known, so we have α k = 0, which means performing a multistart-type iteration. Then, as the iterations go, p k → 0, and α k → 1. In practice the probabilities p k are unknown. One possibility is to estimate them with the information gathered from previous iterations, as [START_REF] Boender | Bayesian stopping rules for multistart global optimization methods[END_REF] or [START_REF] Piccioni | Stopping rules for the multistart method when different local minima have different function values[END_REF]. In [START_REF] Piccioni | Stopping rules for the multistart method when different local minima have different function values[END_REF], the probability of having found a global optimum at iteration k, provided that in the previous iterations the algorithm has found w k local minima, is:

α k = 1 -p k = k -w k -1 k -1 , then α k → 1 as k → ∞,
and therefore the probability of launching a local search decreases to 0. However, there is still a lot of freedom for choosing the sequence α k . The main difficulty for imposing conditions on the sequence α k comes from the fact that in QGO we are not just proposing to replace distance thresholds by quantiles, but also to use local information during the first iterations and more and more global information as the number of iterations increases. This behavior has as consequence that showing that the probability of launching a local search tends to 0 is quite easy, since asymptotically QGO behaves as BS. However, obtaining the rate at which this convergence occurs is much harder than for methods using a threshold that goes to 0, which use local techniques such as Taylor expansions for asymptotic analysis (see eg. Rinnooy Kan and Timmer, 1987a, Lemma 7).

r k = q k α k D k max = q k 1 D k min = q k 0 Iteration k 2 -norm distances

Conclusion

In this first part of the thesis, we analyzed some possible approaches to the efficient approximation of τ -estimators. The bibliographic review, both on global optimization and on computational robust statistics, suggested that two-phase stochastic algorithms coupled with appropriate stopping conditions, could help improving the efficiency of state-of-the-art algorithms.

Our numerical tests confirmed that stopping conditions are a valuable tool, permitting to shorten the computing time of existing algorithms while keeping the quality of the solutions thereof obtained.

Clustering techniques proved to be useful in low dimension, but their efficacy rapidly deteriorate in middling and high dimension. A deeper study of clustering global optimization in higher dimensions shed some light on the problem. The approach of quantiles besides the advantage of being adaptive to the dimension and the current sampled point, represents also a first step towards non-uniform sampling.

However, further research would be necessary to get more concluding results. There may be many other factors influencing the efficiency of clustering methods in high dimension. Just to mention one, let us cite a result of [START_REF] Ledoux | Probability in Banach spaces : isoperimetry and processes[END_REF] about concentration of measures: let λ denote the uniform measure on the cube, and σ a Lipschitz function with Lipschitz constant L and median M . Then

λ(|σ -M | > t) ≤ ce -t 2 /2s 2
where c is a constant, and s 2 = (2π) 1/2 L 2 . This could suggest that uniform sampling is not a good strategy, since sampling far from the median, thus eventually close to minima, is unlikely. This fact has been noticed by [START_REF] Schoen | Random and quasi-random linkage methods in global optimization[END_REF], who proposed to replace uniform random sampling by deterministic, well distributed, sample points. Nevertheless, the simplicity of the uniform distribution permits to have many explicit results, particularly useful for obtaining usable stopping conditions. The theoretically ideal sampling distribution is given by (I.1.1) as T → 0. Unfortunately, implementable version of this sampling strategy, notably simulated annealing, does not perform better in practice than two-phase algorithms based on uniform sampling.

All the experience collected from numerical experiments seems to indicate that stochastic algorithms are able to provide "good" solutions to global optimization problems. Nevertheless, the quality of the outcome is quite variable and difficult to foresee even for particular problem classes with loose computational-time constraints.

Introduction

The second part of this dissertation is devoted to the study of deterministic algorithms for exact calculation of robust regression estimators. In particular, we shall focus on robust regression estimators defined through L-estimates of scale.

We begin by recalling the definition of L-scales: Let |r(β)| (1) ≤ ... ≤ |r(β)| (n) be the ordered absolute values of residuals (see model on page 3 ). L-estimates of scale are defined as weighted 1 or 2 norms of |r(β)

| (i) , n i=1 a i |r| (i) , or n i=1 a i |r| 2 (i) 1/2 , (II.1.1)
where a i are nonnegative constants. The restriction to this class of estimators is motivated by the fact that the rigidity of scales (II.1.1) gives to the problem a combinatorial structure which permits to solve it to optimality. Besides the exhaustive enumeration approach, that is feasible only for very small problems, there exist many more sophisticated techniques for solving global optimization problems; some of them are suited to a particular subclass of problems, and others are applicable to a wide range of problems. The method to be used here strongly depends on the characteristics of the particular problem under consideration. Therefore, the first step will be to examine the problem of minimization of L-scales from an optimization point of view.

Let us concentrate on the problem

min β∈R d n i=1 a i |r(β)| p (i) (II.1.2)
for p equal to 1 or 2, and where the a i s are such that a i = 1 if i ≤ h, and a i = 0 if i > h for a certain h ∈ {1, ...., n}, as it is the case for the LTS and the LTA estimators (see pages 7 and 8). Problem (II.1.2) can be written as a mixed integer nonlinear program using the fact that for arbitrary where Ch = {w ∈ {0, 1} d , n i=1 w i = h}. It is not difficult to check that we have an analogous formulation for the sum of the largest residuals, replacing minimum by maximum.

r ∈ R n , if r (1) ≤ r (2) ≤ • • • ≤ r (n) denote its ordered components, then h i=1 r (i) = min n i=1 w i r i w ∈ {0, 1} d ,
n i=n-h r (i) = max n i=1 w i r i w ∈ {0, 1} d , n i=1 w i = h = -min - n i=1 w i r i w ∈ {0, 1} d , n i=1 w i = h ,
which for h = 1 amounts to minimizing the ∞ -norm of the residuals.

Hereafter we consider the problem of computing the LTS estimator,

min n i=1 w i r i (β) 2 , β ∈ R d , w ∈ Ch , (II.1.5)
which is a MINLP, quadratic in β and linear in w.

Since the constraint w ∈ {0, 1} n can be equivalently expressed as w 2 i -w i = 0, problem (II.1.5) fits into the polynomial programming framework. [START_REF] Lasserre | Global optimization with polynomials and the problem of moments[END_REF] introduced an innovative approach to polynomial global optimization that consists in approximating the problem by a hierarchy of efficiently solvable relaxations. The relaxation is exact for a sufficiently high order. Though, the large size of the relaxations for moderate-sized problems urged researchers on to exploit sparsity in order to diminish the size of the relaxed problems. [START_REF] Lasserre | Convergent SDP-relaxations in polynomial optimization with sparsity[END_REF] and [START_REF] Waki | Sums of squares and semidefinite program relaxations for polynomial optimization problems with structured sparsity[END_REF] went in that direction by putting forward lightweight relaxations under structured sparsity assumptions. Nevertheless, even if sparsity is present, the relaxation of order r of problem (II.1.5) involves (d + 1) 2 n 3+2r 2r variables; so increasing by one the number of observations will result in (d + 1) 2 3+2r 2r new variables added to the relaxation of order r. To get an idea of the sizes, in dimension d = 10 adding one observation would increase the number of variables of the first 5 relaxations by 1210, 4235, 10164, 19965 and 34606 respectively. The methodology of Lasserre ( 2001) is intended to deal with a very wide class of problems; moreover it is a generalization of the «lifting »procedure. Our objective will be to gain efficiency by exploiting the characteristics of our particular problem.

Let us start by cutting problem (II.1.5) off as

min{ v(w) | w ∈ Ch }, (II.1.6)
where v(w) is the marginal value of problem (II.1.5), obtained by minimization over

β ∈ R d for fixed w ∈ Ch , v(w) = inf β∈R d n i=1 w i r i (β) 2 . (II.1.7)
Problem (II.1.6) is a combinatorial optimization problem over the discrete set Ch . In order to avoid the exhaustive enumeration of all of its elements, [START_REF] Agulló | New algorithms for computing the least trimmed squares regression estimator[END_REF] proposed a branch-and-bound (BB) algorithm. A BB algorithm consists in enumerating in a tree the elements of Ch , disregarding unpromising ones. In Figure II.1.1 we depict the BB tree for a small example with the LTS estimator keeping 3 observations out of 6. at the right is associated to the subset of observations 6, 4 and 3, and that at left to observations 5, 4 and 3. In terms of the optimization variable w, they are associated to the points (0, 0, 1, 1, 0, 1) and (0, 0, 1, 1, 1, 0) respectively. The BB algorithm proceeds by examining the branches from right to the left, from top to the bottom; keeping track of the best solution obtained up to that moment, called incumbent solution and denoted by w in the sequel. At each node the algorithm estimates the least function value that can be attained by a leaf belonging to the subtree rooted at that node. If this lower bound exceeds the function value of the incumbent solution, then we can safely discard that subtree, since we are sure that it will not improve the incumbent solution.

A subtree rooted at a node in the BB tree is associated to a set of partially prescribed points of Ch . For instance, the branch in Figure II.1.2, rooted at the node (0, 0, 1, 0, 0, 0) , consists of all the points in {0, 1} 6 with (0, 0, 1) at the first coordinates. Let us write these node constraints as N w = ψ. Note that N is a projection operator; in our example N w = (w 1 ; w 2 ; w 3 ) and ψ = (0, 0, 1) . Therefore, that branch can be ignored if the value of the problem

min v(w) s.t N w = ψ, n j=1 w j = h, w ∈ {0, 1} n , (II.1.8)
is greater than the function value of the incumbent solution, which is an upper bound for the optimal value of Problem (II.1.6). However, obtaining lower bounds for a problem like (II.1.5) is not easy, because of the integer variables and the nonlinearity/nonconvexity. In practice, the value of problem (II.1.8) is essentially as difficult to compute as the value of the original problem, this is why we only settle for lower bounds. The efficiency of a BB algorithm is intimately related to the tightness of the lower bounds.

The algorithm of Agulló ( 2001) is based on a monotonicity property of the LTS problem: in the LTS tree (see Figure II.1.1) every node has greater function value than its father; this is because adding observations to the ordinary least-squares regression can only increase the sum of the squared residuals. Thus the function value at a node serves as lower bound for the whole subtree rooted at it. At each node the algorithm compares the function value of its children to that of the incumbent solution, hereafter denoted U B, and decides whether to prune it or not. If w + is one of those child nodes, and |I + | is the number of nonzero entries in w + , then:

1. If |I + | < d, since in dimension d we can always find an hyperplane that passes through |I + | points, then v(w + ) = 0 < U B. Thus, it cannot be discarded and the analysis continues with its children.

If |I

+ | ≥ d, evaluate v(w + ).
(a) If v(w + ) < U B and w + is a leaf, then a new incumbent solution have been found; set U B = v(w + ) and w = w + . Otherwise repeat the analysis with its children.

(b) If v(w + ) > U B discard w + and all its children.

The monotonicity lower bound presents two flaws:

• It does not provide lower bounds for nodes at the d highest levels.

• It does not look ahead. At each node the increase in function value due to the inclusion of the children is exactly computed, but the contribution of the rest of the observations added in the path to the leaves is underestimated by 0. For instance, if n = 10 and h = 6 it says "the sum of the squared residuals of a regression over six observations comprising observations 2, 4, 5 and 8 will be greater than the sum of the squared residuals of the regression over observations 2, 4, 5 and 8 themselves" without quantifying the minimum increase in sum of squared residuals possible due to the incorporation of two more observations.

The objective of the next two chapters is to obtain alternative lower bounds for problem (II.1.8) when far from the leaves, in particular at the d upper levels. In Chapter II.2 we introduce cuts for reducing the feasible domain of problem (II.1.8), thus tightening the lower bounds. In fact, since the objective function in (II.1.7) is linear in w, the function v is the pointwise minimum of linear functions

v(•) = min β∈R d l β •, with l β = (r 1 (β) 2 , r 2 (β) 2 , ..., r n (β) 2 ),
and consequently it is a closed concave function (Hiriart-Urruty and Lemaréchal, 1993a, Prop. IV.2.1.2). A first consequence of the concavity of the objective function v is that minimization over Ch = {w ∈ {0, 1} n , e w = h} is equivalent to minimizing over the polytope C h = {w ∈ [0, 1] n , e w = h}. The reason is that the global minimum of a concave function over a convex set is always attained at an extreme point; in the case of a polytope, at a vertex. The cuts obtained in Chapter II.2 are based on the concavity of the function v, and moreover, they are called concavity cuts. If we represent the cuts under the form of a linear system Cw ≤ b, then the value of the problem (II.1.6) is the same as the value of

min v(w) s.t Cw ≤ b, n j=1 w j = h, w ∈ [0, 1] n ,
(II.1.9) and the least value at a BB node given by problem (II.1.8) is the same as

min v(w) s.t N w = ψ, Cw ≤ b, n j=1 w j = h, w ∈ [0, 1] n . (II.1.10)
The reason is that concavity cuts eliminate portions of the feasible domain that cannot host a global minimum. Note that the values of problems (II.1.6) and (II.1.8) equal the value of problems (II.1.9) and (II.1.10) respectively, but relaxations of problems (II.1.9) and (II.1.10) will have in general a higher value than relaxations of (II.1.6) and (II.1.8).

Our approach for obtaining lower bounds for (II.1.9) is to cast it as a multiquadratic problem, in order to take advantage of techniques specifically conceived for quadratic problems. Using a by-now classic technique, we lift the problem into a space of matrices, to obtain a linear problem with second order cone constraints, which is a relaxation of (II.1.9). By the way, under the same formulation and using a surrogate of a rank constraint that can be be decomposed as the difference of two convex functions, we shall obtain global optimality conditions for problem (II.1.9).

Notations

In the sequel, we shall always use lowercase letters to name vectors and uppercase letters to denote matrices. Vertical concatenation will be indicated by ';' and horizontal concatenation by ','. For instance, (0; 1) is a column vector in R 2 , and (0, 1) = (0; 1) is the row vector with the same entries. Similarly, if v1 and v2 are two column vectors of the same length, then (v1, v2) is the matrix containing the vectors v1 and v2 as columns, and (v1 ; v2 ) is the transpose of that matrix, containing v1 and v2 as rows. The column vector of length n with all components equal to 1 will be denoted as e, it will be used notably for expressing the linear constraint n j=1 w j = h in matrix form as e w = h. For 1 ≤ p ≤ ∞, B p denote the unit ball for the p-norm. The space of real matrices of size s 1 × s 2 will be denoted by R s 1 ×s 2 , and the space of square symmetric matrices of size s × s by S s . The space R s 1 ×s 2 will be equipped with the inner product A, B = trace(B A), where trace denotes the sum of the diagonal elements, as usual. Given a vector v, D(v) will denote the diagonal matrix with diagonal terms D ii = v i .; reciprocally, for a matrix A, d(A) is the vector whose components are the diagonal elements of A, d(A) i = A ii . A matrix of zeros of size s 1 × s 2 will be denoted by 0 s 1 ×s 2 ; if s 1 = s 2 we just write 0 s . We will make use later of two well-known objects from convex analysis. For a compact set E, δ(x | E) = sup y∈E x y is the support function of the set E. The indicator function of an arbitrary set E, denoted as I E , is the function defined as

I E (x) = 0 if x ∈ E +∞ if x / ∈ E.
For a function f finite at a point x, and ε ≥ 0, the ε-subgradient of f at x, denoted ∂ ε f (x), is the set of vectors s such that

f (y) ≥ f (x) + s (y -x) -ε for all y.
For the particular case ε = 0 we shall use the notation ∂f (x) rather than ∂ 0 f (x).

Chapter II.2

Concavity cuts for LTS

II.2.1 Tuy's cuts for concave minimization over polytopes

Concavity cuts or Tuy's cuts, introduced by [START_REF] Tuy | Concave programming under linear constraints[END_REF], are a domain-reduction technique for concave minimization over a polytope C,

min w∈C v(w). (II.2.1)
Given a vertex of the polytope and a level γ, the concavity cut thereof obtained discards a portion of C with a function value not better than γ. Lower values of γ give raise to deeper concavity cuts, therefore γ is usually set as low as possible in such a way that the concavity cut does not discard vertices with lower function value than the incumbent solution.

Figure II.2.1 illustrates the construction of a concavity cut. The first step consists in finding a sub-optimal vertex w 0 to eliminate, and a level γ such that v(w 0 ) > γ (Fig. . Then we find the hyperplane defining the concavity cut as follows. Let w 0 be the vertex of C to be eliminated and let us denote by u 1 , u 2 , ..., u n ∈ R n the directions of the edges of C emanating from w 0 , and by cone(u 1 , ..., u n ) their conical envelope { n i=1 λ i u i , λ i ≥ 0} . Then,

K(w 0 ) := w 0 + cone(u 1 , ..., u n )
is the smallest cone vertexed at w 0 which contains C . To derive a concavity cut we consider the level set L(γ) := {w ∈ R n |v(w) ≥ γ}, which is commonly assumed for the sake of simplicity to be closed and bounded. Note that since v is concave

L(γ) is convex (Fig. II.2.1-a). Let z i ∈ C
and τi be such that z i = w 0 + τi u i are the intersection points of the cone edges

E i (τ i ) := x 0 + τ i u i , τ i ≥ 0 with the boundary of L(γ) (Fig. II.2.1-b).
Then we determine the hyperplane described by the equation c (w -w 0 ) = 1, which intersects the edges of K(w 0 ) in L(γ), and contains at least n of the points z i (Fig. . This hyperplane is given by

c = 1 τ 1 , 1 τ 2 , ..., 1 τ n U -1 (II.2.2) w * w L(v) u 1 u 2 w 0 a) w * w z 1 z 2 w 0 b) w * w z 1 z 2 w 0 c) w * w z 1 z 2 w 0 d)

II.2.2 Concavity cuts for the LTS problem

In this section we discuss the construction of concavity cuts for minimization of the function v defined in (II.1.7) over the polytope C h = {w ∈ [0, 1] n , e w = h}. There are two issues that hinder obtaining concavity cuts for problem (II.2.1),

1. The polytope C h is degenerate.

2. The level sets of v are unbounded.

The polytope C for n = 3 and h = 2 is depicted in Figure II.2.2, where each vertex corresponds to a subset of 2 observations. In general, the polytope C lies on the hyperplane H defined by the constraint e w = h, and any vertex has (n -h)h neighbors. Instead of making the effort to obtain n-1 directions to define a flat cone on the hyperplane H, we will ignore the constraint e w = h and we will construct a concavity cut on the unit cube. The resulting cut constructed on the cube will intersect H unless H and the cut are parallel, which is impossible as we will see later. Incidentally, the BB tree nodes, other than the leaves, correspond to vertices of the unit cube not belonging to the hyperplane e w = h. Let w 0 be an extreme point of [0, 1] n and I 0 = {j | (w 0 ) j = 1} be the index set of the components equal to one.

Let us define the directions u 1 , ..., u n as u j = -e j if j ∈ I 0 e j otherwise, (II.2.3) where e j is the jth euclidean unit vector (whose jth component is equal to 1 and 0 otherwise). We have that w 0 + cone(u 1 , ..., u n ) is the smallest cone pointed at w 0 that contains the unit cube. In Figure II.2.3 we see an example of cut on the unit cube for the situation of Figure II.2.2 with n = 3 and h = 2. The next step to construct the concavity cut on the cube is to determine the intersection point of the rays w 0 + τ j u j , j = 1, ..., n with the boundary of L(γ), for γ < v(w 0 ). In other words, we look for τj , j = 1, ..., n such that v(w 0 + τj u j ) = γ.

For j ∈ I 0 , by moving along the direction -τ j e j we will downweight observation j , and we will eventually meet the boundary of L(γ). On the contrary, if j / ∈ I 0 we will move along the direction τ j e j , which amounts to add an observation with weight τ j , and we will always increase the objective value, for any τ j > 0. This means that the set L(γ) in unbounded in the direction u j for j / ∈ I 0 ; in this case we will take τ j = +∞. In fact, we are looking for the largest conical set vertexed at w 0 , completely contained in L(γ) and possessing a concise description. Of course, when the level sets are bounded, such a polyhedron can be obtained as described hereabove.

The remaining computations in (II.2.2) are:

• Obtaining the 'step sizes' τj for j ∈ I 0 .

• Inverting the matrix U that contains the directions u j as columns.

As we work on the unit cube, U is simply proportional to the identity matrix, thus inverting that matrix is not an issue. So we turn next to the computation of τj for j ∈ I 0 .

Given that the function v is defined in (II.1.7) as the value of the convex minimization problem inf The next step is to give a formula for v(w 0 + τ j u j ), j ∈ I 0 . Let us compute v(w -τ j e j ). From the definition of matrices M and M we have that M (w -τ j e j ) = M (w) -τ j x j x j , M (w -τj e j ) = M (w) -τ a j a j , where a j = (x j ; y j ) ∈ R n+1 . Then,

v(w -τ j e j ) = det( M (w))(1 -τ a j M (w) -1 a j ) det(M (w))(1 -τ x j M (w) -1 x j ) = v(w) 1 -τ a j M (w) -1 a j 1 -τ x j M (w) -1 x j
and using that (Agulló, 2001, p. 429)

a j M (w) -1 a j = r j (w) 2 /v(w) + x j M (w) -1 x j ,
where r j (w) 2 = (y j -x j β w ) 2 is the square of the jth residual, we obtain v(w 0 -τ j e j ) -v(w 0 ) = -τ j r j (w 0 ) 2 1 -τ j x j M (w 0 ) -1 x j .

(II.2.7)

Finally, using the notation ∆ γ -v(w 0 ), we obtain an explicit formula for the step size τj such that v(w 0 -τj e j ) = γ τj = -∆ ∆x j M (w 0 ) -1 x j -r j (w 0 ) 2 .

(II.2.8)

Combined with the general theory presented in Section II.2.1, these calculations yield to the following result:

Proposition 1 Let γ be a positive real number, and let w 0 be a point in C h such that γ < v(w 0 ) and det(X D(w 0 )X) = 0. Then the vector c ∈ R n defined as

c j = r j (w 0 ) 2 -(γ -v(w 0 ))x j (X D(w 0 )X) -1 x j γ -v(w 0 ) (II.2.9)
if (w 0 ) j = 1, and c j = 0 otherwise, defines a valid γ-cut for the function v over C h , i.e

C h ∩ {c (w -w 0 ) ≤ 1} ⊆ {v ≥ γ}.
We would like to stress the fact that a closed-form expression for the concavity cut as in (II.2.9) is quite unusual. It is crucial in our calculations to work on the unit cube, firstly because we do not need to inverse the matrix U , and also because for obtaining the step-sizes τj we have used the Sherman-Morrison formula for the inverse of a rank-one perturbation of a matrix. We see from (II.2.9) that the only possibility to obtain a cut parallel to the hyperplane e w = h is to take w 0 = e, which will never be the case in a BB algorithm, since all nodes have at most h non-zero components. In general the situation will be that of Figure (II.2.3).

II.2.3 Numerical experiments

In order to get a first idea of the efficiency of the concavity cuts for the LTS problem, we perform the following experiment. We randomly generate a standard normal sample X of size n in R 2 , and a response y i = θ 0 + Xθ + ε with θ 0 = 1, θ = (1; 1) and ε following a standard normal distribution. The number n is chosen small enough to be able to generate all the vertices of C h , corresponding to the number of subsamples of size h. Then, we derive some cuts and keep track of the percentage of vertices that have been eliminated so far.

In our implementation, we store the cuts as a linear system by defining

C =    -c 1 . . . -c nc    , b =    -(1 + c 1 w 1 ) . . . -(1 + c nc w nc )    ,
where w 1 , ..., w nc are the vertices used to construct the cuts. With this notation, points w ∈ C such that Cw > b are eliminated by concavity cuts. First of all, we need a "good" initial solution, whose function value will be used as γ level for generating concavity cuts. In order to obtain such a "good" initial point, we randomly generate a number of candidates and apply to each of them a "concentration step", that consists in repeating iteratively the following steps [START_REF] Rousseeuw | Computing lts regression for large data sets[END_REF] Given a candidate w:

• Compute β as the LS fit to the observations selected by w, the residuals vector r and the ordered residuals

|r π(1) | ≤ |r π(2) | ≤ • • • ≤ |r π(n) |, where π is the ordering, |r π(i) | = |r| (i) .
• Set wπ(1) = wπ(2) = • • • = wπ(h) = 1 and 0 otherwise.

The best of the generated candidates after the previous procedure becomes the incumbent solution w. Alternative methods for obtaining the improved candidate are presented in [START_REF] Agulló | New algorithms for computing the least trimmed squares regression estimator[END_REF] and [START_REF] Schyns | Relaxmcd: Smooth optimisation for the minimum covariance determinant estimator[END_REF]. Here, we described the simplest one.

Iterate

• Derive a concavity cut from w k , add it to C, b.

• Let s C be the row vector of the sums of each column of C. Set the component j of w k equal to 1 if j is one of the h smallest components of s C .

• Test whether Cw k ≤ b, if so iterate, otherwise break.

The results are shown in Tables II.2.1 and II.2.2. In Table II.2.1 we show the results obtained when fixing h = 17 and varying the number of observations from n = 21 to n = 25. In Table II.2.2 the parameter h changes as a function of n as h = round(3n/4), for n in the same range. These results show that for the tested problem the bundle of cuts obtained from the very simple Procedure 4 rules out in all cases more than a half of the subsamples. This evidence suggests that the feasible domain of problem II.1.10 (with cuts) would be significantly narrower than that of problem II.1.8 (without cuts).

II.2.4 Some final comments

In the context of a BB algorithm, concavity cuts will be derived from the leaves that are visited and from nodes that are pruned by the lower bound criterion. The distance from a vertex w 0 to the hyperplane c (w -w 0 ) = 1 equals 1/ c , and is often used as a measure of the depth of the cut. Therefore, "sparse" vertices would result in deeper cuts and it is expected that cuts from pruned nodes would yield deeper cuts than cuts from leaves.

In large trees, adding one cut each time a node is pruned or a leaf is explored can result in a very large number of cuts. Numerical experiments reported by [START_REF] Alarie | Concavity cuts for disjoint bilinear programming[END_REF] suggest that better results are obtained by keeping only a fraction of the cuts. It is suggested to drop cuts that are nearly collinear to existing ones.

Chapter II.3

The bilinear formulation and SDP-type relaxations

In this chapter we reformulate the MINLP (II.1.9), whose optimal value equals the sum of trimmed squares of the residuals of the LTS estimator, with its domain reduced by concavity cuts, as a multi-quadratic problem. The objective of those relaxations is to obtain lower bounds for the optimal value of (II.1.9), or lower bounds for the least value of a branch in the BB tree when they are complemented with node constraints.

Recall from (II.2.5) that, for w such that det(X D(w)X) = 0, v(w) = y D(w)r(β w ), for the unique β w satisfying X D(w)r(β w ) = 0. As a consequence, the function v can be equivalently expressed, using that y D(w) = w D(y), as v(w) = min Therefore, problem (II.1.9) is equivalent to the quadratically constrained quadratic problem (also called multi-quadratic problem, or even

Q 2 P ) min w,β w D(y)r(β) s.t X D(w)r(β) = 0, Cw ≤ b, e w = h, w ∈ {0, 1} n , β ∈ R d .
(II.3.2) Problem (II.3.2) involves only linear terms in w and bilinear terms of the form w k β j . It is possible, using a technique dating back to [START_REF] Glover | Improved linear integer programming formulations of nonlinear integer problems[END_REF], to linearize crossed products of integer and continuous variables by introducing a new variable R kj along with a set of linear constraints to enforce that R kj equals w k β j at all binary realizations of w. This technique was used in [START_REF] Adams | Mixed-integer bilinear programming problems[END_REF] to treat mixed-integer bilinear programs under decoupled constraints w ∈ P 1 ∩ {0, 1} n , β ∈ P 2 for two bounded polyhedral sets P 1 and P 2 . Nevertheless, those linearizations ask the user to provide explicit bounds on β, which are then incorporated in the linearized problem. If such bounds are not available, which is the case of general regression problems, the problem can be solved with a "big box" constraint, but the introduction of such large numbers results in weaker continuous relaxations (c.f. [START_REF] Adams | Comparisons and enhancement strategies for linearizing mixed 0-1 quadratic programs[END_REF]. Note that in our problem, obtaining such a bound before computing the robust regression coefficient would mean that we can obtain reliable information about the "true" or "clean" regression coefficient β from the eventually contaminated data X, y. For these reasons, in the next section we explore relaxations that suit better the structure of our problem.

II.3.1 SDP and SOC relaxations

In semidefinite programming (SDP) one minimizes a linear function, possibly under linear constraints, over the cone of positive semidefinite matrices. Although semidefinite programs are a pretty general class of problems, they can be solved very efficiently by interior-point methods. Semidefinite programming has drawn a lot of attention in the last years, mainly because of its success in obtaining tight lower bounds for hard combinatorial optimization problems. [START_REF] Vandenberghe | Semidefinite programming[END_REF] give an excellent survey on the theory and applications of semidefinite programming. In this section we analyze the SDP relaxation of problem (II.3.2), and then we shall propose a derived relaxation which is better adapted to our problem.

We start by putting problem (II.3.2) in quadratic form. The objective function w D(y)r(β) and the constraint X D(w)r(β) = 0 can be written as z Q0 z + q 0 w, and z Ql z + q l w = 0 l = 1, .., m, respectively, with z = (β, w) ∈ R ζ , ζ := n + d; and using, for l = 0, 1, ..., d, the block matrix

Ql =      0 d M l M l 0 n     
, where the n × d blocks M l and the n × 1 vectors q l equals M 0 = -1 2 D(y)X, q 0 = y 2 , for l = 0, and for l = 1, ..., d,

(M l ) ij = X il X ij , (q l ) i = y i X il .
For problem (II.3.2) in the form:

min w,z z Q0 z + q 0 w s.t z Ql z + q l w = 0, l = 1, .., m, Cw ≤ b, e w = h, w 2 j = w j j = 1, ..., n, z ∈ R ζ , (II.3.3)
the SDP relaxation is obtained by lifting problem (II.3.3) in a matrix space, using the fact that z Qz = trace(z Qz) = trace(Qzz ) = Q, zz .

Note that since the variable z is partitioned as z = (β, w), then the rank-one matrix zz has the following block structure:

zz = ββ βw wβ ww . (II.3.4)
Next we replace the product of vector variables zz by the matrix variable

Z ∈ S ζ , Z = B R R W (II.3.5)
with B ∈ S d , R ∈ R n×d and W ∈ S n , and we add the constraint Z = zz . The constraints w 2 j = w j become W jj = w j , which can be written compactly as d(W ) = w. Thus, the multi quadratic problem with cuts (II.3.3) is equivalent to

min w∈R n , Z∈R n×d Q0 , Z + q 0 w s.t Ql , Z + q l w = 0, l = 1, .., m, Cw ≤ b, e w = h, d(W ) = w, Z -zz = 0. (II.3.6)
Let us introduce, for l = 0, 1, ..., d, the matrices

Q l = 0 q l /2 q l /2 Ql =        0 d+1 q l M l q l M l 0 n        .
The SDP relaxation (which coincides with the Lagrangian relaxation) is obtained by relaxing the non-convex constraint Z -zz = 0 to Z -zz 0 , which can be expressed equivalently as Z 0, with a matrix Z ∈ S ζ+1 structured as

Z =   1 β w β B R w R W   .
In this way we end up with the SDP relaxation of problem (II.3.2):

(SDP )

                     min w∈R n , Z∈S ζ+1 Q 0 , Z s.t Q l , Z = 0, l = 1, .., m, Cw ≤ b, e w = h, d(W ) = w, Z 0.
(II.3.7)

Since the data matrices Q l involved in problem (II.3.7) are sparse, the d × d block B and the off-diagonal elements of the n × n block W appears only in the semidefinite constraint. Thus, the solution to problem (II.3.7) cannot be unique, since for a given feasible pair (w, Z), adding a positive-definite d × d matrix to the block B of Z will keep the pair feasible without changing the objective value of the pair. Besides of introducing a number of unnecessary variables, the non-uniqueness is undesirable from a computational point of view. [START_REF] Kim | Second order cone programming relaxation of a positive semidefinite constraint[END_REF] introduced a transparent approach for dealing with sparsity. It consists in replacing the SDP constraint by second order cone (SOC) constraints that are implied by the SDP constraint. The second order cone K u ⊆ R u is defined as K u = {z ∈ R u : z 1 ≥ (z 2 ; z 3 ; ...; z u ) }; a constraint of the type z ∈ K u is called a SOC constraint. The condition Z 0 implies that any principal submatrix of Z is positive semidefinite. Saying that each 1 × 1 principal submatrix is positive semidefinite amounts to imposing the positivity constraints Z ii ≥ 0, i = 1, ..., ζ. The 2 × 2 principal submatrices are positive semidefinite if and only if (Z kj ) 2 ≤ Z kk Z jj for any pair of indices (k, j). Together, they are equivalent to the SOC constraints

Z kk -Z jj 2Z kj ≤ Z kk + Z jj , or   Z kk + Z jj Z kk -Z jj 2Z kj   ∈ K 3 .
In our problem, for k = 1, d + 2 ≤ j ≤ ζ + 1 we obtain w 2 j ≤ W jj , which, combined with the constraint d(W ) = w, will imply that w 2 j ≤ w j , then 0 ≤ w j , W jj ≤ 1. Note that the introduction of the matrix variable W becomes useless, since it only appears due to the constraint w 2 i -w i = 0, which will be relaxed to w

i ∈ [0, 1]. If we choose 2 ≤ k ≤ d + 1 and d + 2 ≤ j ≤ ζ + 1, we obtain R 2 kj ≤ B kk W jj .
This constraint expresses the logical implication W jj = 0 ⇒ R kj = 0, bringing into play the diagonal of the matrix variable B. Ignoring inequalities involving the offdiagonal elements of the blocks B and W , we obtain the SOC relaxation of problem (II.3.2):

(SOC)                        min 2 M 0 , R + q 0 w s.t 2 M l , R + q l w = 0, l = 1, .., m, Cw ≤ b, e w = h, ( k -w j ; 2R kj ) ≤ k + w j , k = 1, ..., d; j = 1, ..., n, w j ≤ 1, j = 1, ..., n, w ∈ R n , R ∈ R n×d , ∈ R d . (II.3.8)
This time we keep only the diagonal of the block B in the vector variable ∈ R d ; we completely drop the block W , and we introduce directly the constraint w j ∈ [0, 1] instead. Problem (II.3.8) is a convex problem, and can be efficiently solved using interior points methods (see [START_REF] Alizadeh | Second-order cone programming[END_REF], for a comprehensive exposition of SOC programming). The only annoying point with it is the non-uniqueness of the variable . This inconvenient can be eliminated by adding a term of the form µe , for small µ > 0 to the objective function in order to achieve uniqueness of without perturbing the solution of the problem. Note that if we had a bound for the variable β, say

|β k | ≤ M k , we could eliminate the variable replacing the SOC constraint R 2 kj ≤ b k w j by R 2 kj ≤ M 2 k w j .
This is in fact one of the linear inequalities added by the above-mentioned linearization proposed in [START_REF] Adams | Mixed-integer bilinear programming problems[END_REF].

The SOC relaxation can be seen as a compromise between the SDP relaxation that gives tight bounds at a strong computational cost, and linear programming relaxations, which are cheaper from a computational point of view but give much weaker bounds.

II.3.2 Global optimality conditions

The objective of this section is to describe a reformulation of problem (II.3.2) that keeps only the relevant crossed terms w j β k in bilinear programs. It is based on a convenient writing of the condition R = wβ as a rank-one constraint, followed by a formulation of that constraint as the difference of two convex functions.

Let us reconsider problem (II.3.8), with the SOC constraints replaced by the exact definition of R as wβ , d+1) ,

min 2 M 0 , R + q 0 w s.t 2 M l , R + q l w = 0, l = 1, .., m, Cw ≤ b, e w = h, R = wβ , 0 ≤ w j ≤ 1, j = 1, ..., n, w ∈ R n , R ∈ R n×d , β ∈ R d .
min A 0 , Z s.t A l , Z = 0, l = 1, ..., m, CZ(2 : n + 1, 1) ≤ b, n+1 j=2 Z j1 = h, Z 11 = 1, 0 ≤ Z j1 ≤ 1, j = 2, ..., n + 1, rank(Z) = 1, Z ∈ R (n+1)×(
where A l =   0 0 1×d q l 2M l   , l = 0, 1, ..., m, and Z(2 : n + 1, 1) denotes the n last elements of the first column of the matrix Z. We introduce two additional matrices A m+1 and A m+2 to express the constraints n+1 j=2 Z j1 = h and Z 11 = 1 respectively,

A m+1 =     0 0 1×d e 0 n×d     , A m+2 =     1 0 0 0 ,     ,
and the matrices C i , i = 1, ..., nc, to express the constraint CZ(2 : n + 1, 1) ≤ b as C i , Z ≤ b i , i = 1, ..., nc. Each matrix C i has only zeros, except for its first column that equals (0; c i ), where c i is the ith concavity cut (stored in the matrix C).

Let A 1 : R (n+1)×(d+1) → R m+2 and A 2 : R (n+1)×(d+1) → R nc be defined as

A 1 Z =( A 1 , Z ; ...; A m , Z ; A m+1 , Z ; A m+2 , Z ), A 2 Z =( C 1 , Z ; ...; C nc , Z ).
and f ∈ R m+2 be the column vector defined by f i = 0, i = 1, ..., m and f m+1 = f m+2 = 1.

Then we can rewrite problem (II.3.9) in compact form as: II.3.10) where ∆ ⊆ R (n+1)×(d+1) is the convex set

           min A 0 , Z s.t A 1 Z = f, rank(Z) = 1, Z ∈ ∆, (
∆ := {Z | A 2 Z ≤ b, 0 ≤ Z j1 ≤ 1 for j = 2, ..., n + 1}.
Note that problem (II.3.10) is an exact reformulation of problem (II.3.2) with the constraint w j ∈ {0, 1} relaxed to w j ∈ [0, 1]. It is a linear problem in the matrix variable Z except for the rank-one constraint. The rank function is non-convex and not even continuous; in fact it is integer-valued and constant along rays. For these reasons, the study of optimization problems involving the rank function is quite recent. The Eckart-Young Theorem is possibly the only 'old' optimization result involving the rank function. It gives an explicit solution to the problem of minimizing the distance from an arbitrary matrix to the set of matrices of a given rank. A significant breakthrough in the study of the rank function has been made by [START_REF] Fazel | Matrix rank minimization with applications[END_REF], by computing the convex envelope of the rank function restricted to a ball. This result is of great importance for problems involving the rank in the objective function. In the remaining of this section we shall focus rather on relaxations of optimization problems under rank constraints, paying particular attention to the rank-one constraint.

II.3.2.1 The rank-one constraint as the difference of convex functions

Let a given n × d matrix A have the singular value decomposition

A = U D(ς)V ,
where U and V are orthogonal matrices of size n × d and d × d respectively, and ς ∈ R d contains the singular values in decreassing order:

ς 1 ≥ • • • ≥ ς d ≥ 0.
Let φ be a norm in R d , then we can define an associated matrix norm in R n×d as A φ = φ(ς).

(II.3.11)

By putting φ = • 2 we obtain the norm associated to the Frobenius inner product

A 2 F = trace(A A) = d i=1 ς 2 i .
A very important case in the sequel will be φ = • 1 , which originates the nuclear norm

A * = d i=1 |ς i |.
The following function will play an important role in the study of the rank-one constraint:

Definition 1 Let us define the continuous surrogate of the rank as the function g : R n×d → R given by .3.12) This result permits to put the rank-one constraint in the context of dc programming. In this way we have at hand global optimality conditions, local minimization algorithms that exploits that structure and, of course, the advantage of working with continuous, differentiable functions instead of the integer-valued rank function.

g(A) = A 2 * A 2 F , if A = 0 0, if A = 0. ( II

II.3.2.2 Rank constraints as a reverse-convex constraint

In this subsection we shall give yet another formulation of the rank constraint as a reverse convex constraint. We say that a constraint is reverse convex if it can be expressed as g(A) ≥ 0 for a convex function g, or equivalently as h(A) ≤ 0 for a concave function h. This formulation has the advantage of describing the set {A | rank(A) ≤ k} for any 1 ≤ k ≤ d. Its weak point is that it is primarily formulated for square matrices.

Let us suppose that A is a square positive semidefinite matrix, and let ς 1 ≥ • • • ≥ ς d ≥ 0 be the eigenvalues of A. Then (Horn and Johnson, 1985, p. 191)

ς k+1 + • • • + ς d ≤ trace(U AU )
for any U ∈ U k , where U k denotes the set of n×(n-k) matrices such that U U = I n . Moreover, equality holds if and only if the columns of U form an orthonormal basis of eigenvectors of the matrix A. In other words, we have

ς k+1 + • • • + ς d = min U ∈U k trace(U AU ).
(II.3.15)

As the function A → trace(U AU ) is linear for any U ∈ U k , then the function

h k (A) = min U ∈U k trace(U AU )
is concave, and by virtue of (II.3.15), we have proved Proposition 4 Let A be a d × d symmetric positive semidefinite matrix. Then

rank(A) ≤ k ⇐⇒ h k (A) ≤ 0. Note that the function ς → ς k+1 + • • • + ς d is concave (use (II.1.
3) to write it as a minimum of linear functions) as a function of the singular values, but the function A → ς k+1 (A) + • • • + ς d (A), where ς(A) denote the vector of singular values of A, is not necessarily a concave function of the matrix A, even for square A (cf. Dacorogna and Marcellini, 1999, Theorem 7.8). For an optimal control problem, [START_REF] Kim | Structurally constrained H 2 and H ∞ control: A rankconstrained LMI approach[END_REF] use the characterization above to devise a penalty method for the rank constraint, since if U ∈ U k then rank(A) > k ⇒ trace(U AU ) > 0, thus a term of the form αtrace(U AU ) for α > 0 and U ∈ U k acts as penalization of the rank constraint violation.

II.3.2.3 Global optimality conditions using ε-subdifferentials

In this section, using the dc formulation of the rank-one constraint (II.3.14), we obtain global optimality conditions for problem (II.3.10).

Let us dualize the rank-one constraint, introducing a multiplier α > 0

min A 0 , Z + α( Z * -Z F ) s.t A 1 Z = f, Z ∈ ∆,
which we rearrange as:

min ( A 0 , Z + α Z * ) -α Z F s.t A 1 Z = f, Z ∈ ∆.
(II.3.16)

Then, for each α > 0, the value of problem (II.3.16) gives a lower bound for the original problem. A matrix Z ∈ ∆ is a local optimum of problem (II.3.16) if there exists vectors µ ∈ R m+2 , λ ∈ R n and ν ∈ R nc such that (II.3.17) where Λ is a matrix of zeros, except for the first column that equals (0; λ). The multipliers λ and ν must also satisfy the complementarity conditions Next we recall the characterization given in [START_REF] Watson | Characterization of the subdifferential of some matrix norms[END_REF] for the subdifferential of norms defined as in (II.3.11), at a matrix A = U D(ς)V , ∂ A φ = conv{U D(s)V , s ∈ ∂φ(ς)}.

A 1 Z = f, A * 1 µ + A * 2 ν + Λ + α∂ Z F ⊂ {A 0 } + α∂ Z * ,
λ i ∈      {0} if 0 < z i < 1, (-∞, 0] if z i = 0, [0, ∞) if z i =
this result extends naturally to ε-subdifferentials

∂ ε A φ = conv{U D(s)V , s ∈ ∂ ε φ(ς)}, ε ≥ 0. (II.3.19)
Similarly, since each q norm is the support function of the unit ball for the dual norm q, with 1/q + 1/q = 1,

• 1 = δ(•|B ∞ ), • 2 = δ(•|B 2 )
we have (see Hiriart-Urruty and Lemaréchal, 1993b, pp. 97)

∂ ε δ(d, C) = {s ∈ C : δ(d, C) ≤ s d + ε}, ε ≥ 0.
Then, we obtain the expression for the ε-subgradient of a norm:

∂ ε d q = {s ∈ B q : d q ≤ s d + ε}.

Note that by Hölder inequality s d ≤ s q d q . Thus, for s ∈ B q, s d ≤ d q , and in the particular case ε = 0 we have an equality:

∂ d q = {s ∈ B q : d q = s d}.

An interesting fact about ε-subgradients of norms, emanating directly from the later characterization, is that, for any ε ≥ 2 ς 1 ,

∂ ε ς 2 = B 2 , and ∂ ε ς 1 = B ∞ .
This will be useful since for dc. programs there exist global optimality conditions which are stated precisely in terms of ε-subgradients.

Proposition 5 A feasible matrix Z ∈ ∆ is a global optimum of problem (II.3.16) if, for every 0 ≤ ε ≤ 2 ς 1 , where ς is the vector of singular values of Z, there exist ε 1 , ε 2 satisfying ε 1 + ε 2 = ε, and vectors λ ∈ N ε 2 , µ ∈ R m+2 and ν ∈ R nc such that, AZ = 0, (II.3.20)

A * 1 µ + A * 2 ν + Λ + α∂ ε/α Z F ⊂ {A 0 } + α∂ ε 1 /α Z * , (II.3.21) D(ν)(b -A 2 Z) ≤ ε 2 , (II.3.22)
where N ε 2 is the set of all vectors ξ ∈ R n whose components verify: (II.3.23) and Λ is as before.

ξ j ∈      {0} if 0 < Z j1 < 1, (-∞, ε 2 ] if Z j1 = 0, [-ε 2 , ∞) if Z j1 = 1,
Proof. Let us denote by H the kernel of the linear operator A 1 . Note that H is a subspace of R (n+1)×(d+1) . Then, rewrite the constraints of problem (II.3.16) using indicator functions, defined as The final form follows by using that ker(A 1 ) ⊥ = Im(A * 1 ), and that ∂ ε 2 I ∆ (Z) = Λ + A * 2 ν.

I E (x) = 0 if x ∈ E, +∞ if x / ∈ E,
Note that in obtaining such a "closed" form of the global optimality conditions we have greatly benefited from the linear structure of our problem, since ∂ ε is a singleton only for affine functions (cf. Hiriart-Urruty and Lemaréchal, 1993b, prop. 1.2.4).

Just to illustrate the mechanism, let us consider the toy example of pure rank minimization over non-zero matrices:

min R∈R n×d \{0} R * -R F ,
whose solution set is obviously the ray of rank-one matrices. The global optimality condition reads:

∂ ε Z F ⊂ ∂ ε Z * .
Let Z = U D(ς)V be a non-optimal candidate, then there exists ε > 0 such that

∂ ε Z F \ ∂ ε Z * = ∅, this can occur only if ∂ ε ς 2 \ ∂ ε ς 1 = ∅, and s ∈ ∂ ε ς 2 \ ∂ ε ς 1 =⇒ s ∈ B 2 , ς 2 ≤ s ς + ε ≤ ς 1
but we know that in general ς 2 ≤ ς 1 with equality only if ς has at most one non-zero component. Proposition 5 is not directly related to the original problem (II.1.5), since it gives a characterization of the global optimum of a transformation of the same problem with the integer constraint relaxed. Its interest comes mainly from a global optimization viewpoint, because it gives global optimality conditions for bilinear problems, for which the SDP relaxation does not apply directly.

Conclusions and perspectives

We have studied some global optimization problems involved in the computation of high breakdown point regression estimators. Due to their advantageous statistical properties, we consider robust estimators based on M -scales and robust estimators based on L-scales. Robust estimators based on M -scales, and more particularly τestimators of regression, can be tuned to have the desired breakdown point and also the desired efficiency, therefore they are preferable from a statistical point of view. A closer look shows that the efficiency of τ -estimators is obtained by giving a positive weight to every observation with small scaled residuals. In contrast, estimators based in L-scales trim a fraction of the observations independently of the value of their residuals. However, in practice the rigidity of the L-scales gives to the optimization problem a combinatorial structure that permits, time considerations aside, to find a global optimum. Concerning robustness, for the LTS estimator the number of leaves in the branch-and-bound tree equals n h , and is maximal for h ∼ n/2, which gives the maximal breakdown point. In general, better estimators involve more difficult global optimization problems.

The first part of the thesis is devoted to the approximation of τ -estimators for robust regression. Considering factors such as the non existence of a characterization for the global optimum of the associated optimization problem and the unboundedness of the feasible domain that impedes an exhaustive inspection, we decided to study two-steps stochastic algorithms based on random subsampling and local searches that permit to take advantage of the differentiability of the problem. In [START_REF] Flores | On the efficient computation of robust regression estimators[END_REF], after reviewing existing algorithms for approximating τ -estimators and observing that most of them can be seen as restricted versions of clustering global optimization algorithms, we investigate the impact of incorporating clustering techniques and stopping conditions. The main conclusions of these extensive numerical tests are that stopping conditions improve the efficiency of existing algorithms, while clustering techniques work well in low dimension, but they harm the efficacy of the algorithm in middling and high dimension. The disappointing behavior of clustering techniques in moderate dimension led us to further investigate the reasons thereof in Chapter I.3. Since most clustering global optimization algorithms rely in some way on the nearest neighbor, we payed particular attention to the gap phenomenon affecting the notion of nearest neighbor in high dimension. We show how the quantile order of distances to points in the sample acts as an index in [0, 1], equal to 0 for multistart and equal to 1 for best start, and that the gap phenomenon pushes algorithms to take extreme values and behave therefore as MS or BS. We propose as an alternative to replace distance thresholds by quantiles. The main strength of this approach, that motivates further research, is that it is independent of the dimension, the feasible domain and the sampling distribution.

Concerning this last point, it is worth to note that clustering global optimization algorithms performing a concentration step (that consists in applying a few iterations of a local minimization routine to each sampled point, see [START_REF] Törn | Global optimization[END_REF]), that performed very well in our tests of Chapter I.2, have been very little studied; mainly because the sample after the concentration steps is no longer uniformly distributed.

The choice of the sequence of quantile orders that characterizes the QGO algorithm remains as an open question. Guidelines for choosing such sequences are needed; the main difficulty in doing so is that our proposal consists not only in replacing radius by quantiles, but also in using more and more global information as the number of iterations increases. As a consequence, local techniques such as Taylor expansions will not be useful for the asymptotic analysis of the algorithm.

The key points of the first part are:

• We position the existing algorithms in the context of clustering global optimization.

• We introduce inexact iterations for accelerating the search for local minima.

• We perform extensive numerical tests in order to evaluate the impact of stopping conditions and of clustering techniques in existing algorithms for approximating the τ -estimator.

• We investigate the consequences of the gap phenomenon in clustering global optimization algorithms.

• We prove the convergence of a new clustering-type algorithm intended to be independent of the dimension, the optimization domain and the sampling distribution.

The main research direction emanating from the work in this part is undoubtedly the search for some criteria permitting to find "good" quantile order sequences for the QGO algorithm. At a longer term, what is needed is a better understanding of the consequences of the interactions between sampling distribution and feasible domain, such as the concentration of measures mentioned at the end of Chapter I.3.

The second part of the thesis is devoted to deterministic algorithms for computing the least trimmed squares regression estimator, which is defined through a nonlinear mixed-integer program. Due to the combinatorial nature of this problem, we concentrated on obtaining lower bounds to be used in a branch-and-bound algorithm. It is known that the associated problem can be written as a concave minimization problem over a polyhedral domain. Under this formulation, we obtain a closed-form expression for concavity cuts, that can be derived at a negligible computational cost. We also perform some preliminary numerical tests showing that concavity cuts eliminate a considerable fraction of the vertices of the feasible polyhedron, which could considerably tighten lower bounds obtained from relaxations of 88CHAPTER II.3. THE BILINEAR FORMULATION AND SDP-TYPE RELAXATIONS the problem. Then, we show that the same problem can also be cast as a bilinear program with bilinear constraints, and study the associated positive-semidefinite and second-order cone programming relaxations. The last section deviates a few from the original motivation, and is devoted to global optimality conditions for the bilinear problem lifted into a matrix space. Before that, we introduce and study a continuous minorant of the rank function that deserves attention by its own, and we derive alternative formulations of some constraints involving the rank function as difference-of-convex (dc) or reverse convex constraints. Finally, we use the dc formulation of the rank-one constraint to obtain global optimality conditions for the bilinear problem.

The key points of the second part are:

• We obtain explicit concavity cuts.

• We give a reformulation of the original mixed-integer nonlinear problem as a bilinear problem with bilinear constraints.

• We propose a second-order cone programming relaxation to obtain lower bounds to be used in a branch-and-bound algorithm, reinforced by concavity cuts.

• We conduct a study of rank constraints, obtaining dc and reverse convex formulations of them.

• We obtain global optimality conditions for the bilinear problem.

We would have loved to test the tightness of the second-order cone programming relaxation in practice through numerical experiments, but unfortunately this (timeconsuming!) task will remain at the place number one of the ToDo list. We expect to continue with the research lines that led to the reformulations of rank constraints, as we hope to contribute to the rapidly growing topic of rank optimization.

Let us conclude by mentioning the possible extensions of the work in this thesis to the computation of robust estimators for multivariate analysis. The experience with two-steps algorithms of Chapter I.2 will be very valuable for devising approximation algorithms for τ -estimators of multivariate location and scatter. In fact, some advances in this direction have been made during a stay of the author at the Department of statistics of the University of British Columbia in Vancouver, Canada. That piece of work (still in progress), in collaboration with Matias Salibian-Barrera and Ruben Zamar, led us to identify some issues specific to the multivariate case.

Figure 1 :

 1 Figure 1: International phone calls from Belgium in the period 1950-1973.

Figure 2 :

 2 Figure 2: Krafft point for 32 chemical compounds.

Figure 3 :

 3 Figure 3: LS regression on a dataset with 20% of contamination.

Figure 4 :

 4 Figure 4: L 1 and LS regressions on a toy example.

Figure 6 :

 6 Figure 6: Two admissible ρ functions.

Figure 7 :

 7 Figure 7: The τ objective function (top), and the τ -estimator fitting (bottom), for the same data as in Figure 3.

  Figure I.2.1: Average times needed to find a solution to (I.2.5) with respect to d, using approximated (IARLS) and inexact (IIRLS5 and IIRLS20) iterations. Sampling Concentration-selection Clustering Stopping

Figure I. 2 . 2 :

 22 Figure I.2.2: Schematic representation of clustering methods, at the sth iteration.

Figure

  Figure I.2.3: An example of subsampling with n = 20 and d = 4. Seven subsamples have been sampled among the 4845 possible subsamples, coming from three of the four cells.

  Figure I.2.4: Schematic representation of Ruppert's version of random resampling.

  Figure I.2.5: Schematic representation of the Fast-S and Fast-τ methods.

  Figure I.3.1: For each iteration k, the upper curve plots the maximum distance from the last sampled point to the points sampled in previous iterations. Similarly, the lowest curve indicates the minimum distance; between them we show the median of the distances.

  Figure I.3.2: The gap q 1 -q 0 (left) and the inter-quartile difference q .75 -q .25 (right) for the ∞ norm.

  Figure I.3.3: The gap q 1 -q 0 (left) and the inter-quartile difference q .75 -q .25 (right) for the 2 norm.

  Figure I.3.4: The gap q 1 -q 0 (left) and the inter-quartile difference q .75 -q .25 (right) for the 1 norm.

Figure

  Figure I.3.5: Behavior of s k in dimensions 20, 35, 50 and 75. Excepting s k , all the plots have the same meaning as in Figure I.3.1

52CHAPTER

Figure I. 3

 3 Figure I.3.6: Above: a sequence of quantile orders α k ; below: radii obtained from the sequence α k .

∅Figure

  Figure II.1.1: The BB tree for n = 6 and h = 3.

Figure

  Figure II.2.1: Construction of a concavity cut.

Figure

  Figure II.2.2: The LTS polytope for n = 3 and h = 2.

Figure

  Figure II.2.3: A concavity cut on the unit cube.

  ) D(w)r(β), then we have, when det(X D(w)X) = 0, v(w) = r(β w ) D(w)r(β w ), for the unique β w satisfying X D(w)Xβ w = X D(w)y. (II.2.4) Since r(β) = y -Xβ, v(w) =r(β w ) D(w)r(β w ) =(y -Xβ w ) D(w)r(β w ) =y D(w)r(β w ) -β w (X D(w)r(β w )), and rewriting (II.2.4) as X D(w)r(β w ) = 0 we conclude that v(w) = y D(w)r(β w ). (II.2.5) Furthermore, as β w = (X D(w)X) -1 X D(w)y, v(w) = y D(w)(y -Xβ w ) = y D(w)(y -X(X D(w)X) -1 X D(w)y), which motivates the introduction of the matrices M (w) = X D(w)X and M (w) = X D(w)X X D(w)y y D(w)X y D(w)y = M (w) X D(w)y y D(w)X y D(w)y ;then by the Schur complement formula for determinants of block matrices we obtaindet( M (w)) = det(M (w)) det(y D(w)y -(y D(w)X)M (w) -1 (X D(w)y)) = det(M (w)) det(y D(w)[y -XM (w) -1 X D(w)y]) = det(M (w)) det(y D(w)[y -Xβ w ]) = det(M (w))y D(w)r(β w ) = det(M (w))v(w), (II.2.6)which supplies us, provided that det(M (w)) = 0, with a convenient expression for the objective function v,

  The cuts are generated by Procedure 4. Procedure 4 Procedure for generating concavity cuts. 1. Draw an initial point w 0 randomly, apply the concentration steps described above to obtain the incumbent solution w. Set w k = w 0 , γ = v( w) and C, b = [].

(

  II.3.9) 78CHAPTER II.3. THE BILINEAR FORMULATION AND SDP-TYPE RELAXATIONS Then we use the fact that R = wβ if and only if rank 1

  (b -A 2 Z) ≤ 0, respectively. Needless to say, if the solution of problem (II.3.16) is a rank-one matrix, then we have obtained the global solution to problem (II.3.10).

  Z * + I H (Z) + I ∆ (Z)] -[α Z F + -A 0 , Z ] .This is a dc program without constraints, whose global optimality conditions are∂ ε [α Z F + -A 0 , Z ] ⊆ ∂ ε α Z * + I {f }+H (Z) + I ∆ (Z)84CHAPTER II.3. THE BILINEAR FORMULATION AND SDP-TYPE RELAXATIONS by using Hiriart-Urruty andLemaréchal (1993b, Prop. XI.1.3.1.(ii) and (vi)) we see that the left-hand of the inclusion equals -A 0 + α∂ ε/α Z F . For the right-hand side, we have from Hiriart-Urruty andLemaréchal (1993b, Prop. XI.1.3.2), that forZ ∈ {f } + H ∂ ε ( Z * + I {f }+H (Z)) = ∂ ε Z * + H ⊥ .Then, since dom( • * + I {f }+H ) ∩ dom(I ∆ ) = ∅, we can apply Hiriart-Urruty andLemaréchal (1993b, Theo. XI.3.1.1.) to obtain,∂ ε [α Z * + I H (Z) + I ∆ (Z)] = ε 1 +ε 2 =ε α∂ ε 1 /α Z * + ∂ ε 2 I ∆ + H ⊥ .

  . Rousseeuw's original version introduced the sampling technique used by most existing algorithms, and his algorithm consisted only of sampling and choosing the candidate with the least scale. The modification of Ruppert, illustrated in Figure I.2.4, included selection and local minimization applied, without clustering, to the best candidates. The details of each step are described below.

Given parameters N , t, do once:

Table I .

 I 2.1: Details of the execution of Multistart (MS) and Fast-τ with 50 samples (FT50). The datasets consisted of n = 400 observations in dimension d = 15, with a fraction of contamination of δ = 40%.

				MS		FT50
	θ	0.3	0.6	0.9	0.95
	w	1.2	1.2	1.7	2.4	2.2
	m	3.19 4.69 19.08 50.24	50
	P m,w	0.47 0.66	0.9	0.95	0.96
	% Slope	95.4 94.6 86.4 77.2	68.6
	T M S /T F T 50 0.17 0.24 0.96 2.52	1
	σMS /σ F T 50 1.012 1.0116 1.0083 1.0045	1

Table I .

 I 2.2: Quality of the solution and computing effort for SC, SL1 and SL2. The parameters of the contamination were l = 2 and δ = 0.2. The relative time and the relative τ -scale are with respect to FT500

			% Clean Slope	Relative τ -scale	Relative Time
	d		95% 97% 99% 95% 97% 99% 95% 97% 99%
		SC	58	58 58.2	1	1	1	0.34 0.35 0.74
	5	SL1 88	88 88.2 1.019 1.019 1.02 0.15 0.16 0.49
		SL2 88.4 88.4 88.4 1.021 1.021 1.021 0.1	0.1 0.29
		SC 64.4 64.4 64.6	1	1	1	0.466 0.51 1.44
	10	SL1 78.8 78.8 85.8 1.025 1.025 1.020 0.21 0.23 0.63
		SL2 78	78	85 1.024 1.024 1.021 0.12 0.13 0.33
		SC 68.8 69.2 69.8 1.003 1.002 0.999 0.48 0.64 2.45
	15	SL1 32.8 33 52.4 1.104 1.103 1.069 0.23 0.24 0.65
		SL2 36.4 36.4 53.8 1.090 1.090 1.064 0.14 0.15 0.4
		SC 64.8 64.8 65.2 1.016 1.014 0.998 0.51 0.72 3.46
	20	SL1 7.6	7.8 19.8 1.193 1.192 1.149 0.25 0.26 0.62
		SL2 8.4	8.4 14.8 1.206 1.206 1.192 0.18 0.18 0.40
	Table I.2.3: Quality of the solution and computing effort for FT250 and FT1500.
	The parameters of the contamination were l = 2 and δ = 0.2. The relative time and
	the relative τ -scale are with respect to FT500
		d			% Clean Slope Relative τ -scale	Relative Time
			FT250	58.2			1.000	0.57
		5	FT1500	58.4			1.000	2.8
			FT250	64.8			1.000	0.62
		10	FT1500	64.6			0.999	2.63
			FT250	71.4			1.001	0.64
		15	FT1500	72.2			0.998	2.57
			FT250	65.4			1.007	0.65
		20	FT1500	69			0.996	2.6

  .2.4. TableI.2.4: Percentage of samples where the minimum found had the clean slope, and time spent relative to FT500. For d = 10.

	δ	10%		15%			20%
	Slope	1.1 1.4 1.7 2	1.1 1.4 1.7 2	1.1 1.4 1.7 2
					% Clean Slope
	SC	62 100 100 100	0 26.4 98 100	0	0 3.2 61.8
	SL	99 100 100 100 44.2 94.6 100 100 2.8 21.4 77.4 96.4
					Relative Time
	SC	0.34 0.33 0.34 0.35 0.37 0.38 0.40 0.39 0.39 0.41 0.44 0.39
	SL	0.18 0.21 0.20 0.20 0.22 0.17 0.20 0.21 0.24 0.21 0.19 0.14
	FT250 0.57 0.59 0.59 0.60 0.60 0.59 0.59 0.61 0.61 0.62 0.61 0.57
	FT1500 2.71 2.70 2.69 2.68 2.58 2.69 2.62 2.57 2.57 2.56 2.60 2.78
					Relative τ -scale
	SL	1.01 1	1	1	1.07 1.04 1	1	1.01 1.04 1.07 1.01

Remerciements

Part II

Deterministic algorithms for mixed-integer bilinear programs.

This function has a close relationship with the rank function, as the following proposition shows:

Proposition 2 The continuous surrogate of the rank, defined in (II.3.12) satisfies the following properties:

G3. If all the non-zero singular values of A are equal, then g(A) = rank(A).

G4. rank(

Proof. Property G1 follows directly from the homogeneity of norms. As before, let

and since all the singular values are non-negative ς 1 = d i=1 ς i . As a consequence ς 2 ≤ ς 1 and g(A) ≥ 1 for any A = 0. Note also that equality holds in (II.3.13) if and only if ς has at most one non-zero component, which is precisely what G4 states. For the upper bound we use that rank(A) = k if and only if

, where e k is the vector with ones in the first k entries and zeros elsewhere. Then, by the Cauchy-Schwarz inequality, ς 1 ≤ e k 2 ς 2 = √ k ς 2 , hence g(A) ≤ k = rank(A) and G2 is proved. For proving G3 it suffices to consider the case

Properties G1 and G2 show that the function g is a minorant of the rank function with the same homogeneity. We will be mainly interested in property G4, which has been proved in Malick (2007, Theorem 1) for the particular case of square symmetric matrices with only ones on the diagonal, which is a corner of the • * -ball of radius n. It can be directly obtained from G3 as well. Reciprocally, the problem with fixed Frobenius norm has been studied by [START_REF] Fazel | Matrix rank minimization with applications[END_REF], where they proved that, for any M > 0, • * /M is the convex envelope of the rank function on the set {A ∈ R n×d | A F ≤ M }. This result has provided a theoretical background to the so-called trace heuristic, that consists in replacing the rank function by the nuclear norm in rank minimization problems. Properties G1 and G4 say that, for A = 0, the rank-one constraint can be expressed as a difference-of-convex (dc) functions constraint, Proposition 3 Let A ∈ R n×d \ {0}, then rank(A) = 1 ⇐⇒ A * -A F ≤ 0.

(II.3.14)