Robust classification via MOM minimization

Introduction

Machine learning has experienced unprecedented growth in recent years with a major societal impact exceeding that it has had in science for several decades. The article is set within the framework of robust learning theory, in the sense that only moment hypotheses are guaranteed on the data and the dataset may contain outliers. Robust learning theory has received particular attention in recent years, the aim sometimes being to construct reliable procedures under weaker assumptions than strong stochastic ones (such as the Gaussian character of noise in regression) and sometimes to detect outliers. This interest can be appreciated, for example, by the challenges recently posted on “kaggle”, the most popular data science competition platform. The 1.5 million dollars problem “Passenger Screening Algorithm Challenge” involves the discovery of terrorist activity from 3D images. The “NIPS 2017: Defense Against Adversarial Attack” consists in building algorithms robust to adversarial data.

Our approach extends Vapnik’s empirical risk minimization (ERM) (cf., ). The ERM is defined as a minimizer of the empirical risk $f\to(1/N)\sum_{i=1}^{N}I_{\{Y_{i}\neq f(X_{i})\}}$ over a class ${\cal F}$ of functions from ${\cal X}$ to $\{\pm 1\}$:

The example in Figure 1 is representative of a general problem. Algorithms approaching ERM do not directly solve the minimization problem (1), which is NP-hard even for simple function classes like half-spaces indicators . Instead, the algorithm seeks to approach a convex relaxation of the problem, in which the considered loss (e.g. hinge or logistic loss) is not necessarily bounded. The advantage of this relaxation is that the function to be minimized is convex and that the minimum can easily be approached using for example gradient descent algorithms. On the other hand, this function becomes unbounded and can be diverted by a single outlier. Minimizing this corrupted criterion then leads to a bad classifier choice.

In learning theory, most alternatives to ERM manage the problem of outliers and heavy tail distributions for outputs only. These solutions are based on the pioneering work of John Tukey , Peter Huber and Frank Hampel , replacing the square loss by a robust alternative like Huber loss or Tuckey’s biweight loss. These methods do not allow to treat the case where the inputs are with heavy tails or corrupted, which is a classical problem in robust statistics also known as the “leverage point problem”, see . In this article, we address this question by considering an alternative to M-estimators, called median-of-means. Several estimators based on this principle have recently been proposed in the literature . To our knowledge, these articles use the small ball hypothesis to treat problems of least square regression or Lipschitzian loss regression. This assumption is restrictive in some classic functional regression frameworks or for problems such as the construction of recommendation system where inputs are sampled in the canonical basis and therefore do not satisfy a small ball condition. In this article, we propose a natural estimator based on the MOM principle, which we call MOM minimizer, and which we study without the small ball hypothesis. Instead we assume an a priori bound on the $L^{2}$-norm of learning functions. We can identify mainly two streams of hypothesis in Learning theory: 1) boundedness with respect to some norm of the class $F$ of functions and the output $Y$, the typical example is the boundedness in $L_{\infty}$ assumption or 2) norm equivalence assumption over the class $F$ (or, more precisely, on the shifted class $F-f^{*}=\{f-f^{*}:f\in F\}$ where $f^{*}$ is the oracle in $F$) and $Y$, the typical example being the subgaussian assumption, i.e. $\left\|f-f^{*}\right\|_{\psi_{2}}\leq L\left\|f-f^{*}\right\|_{L_{2}},\forall f\in F$. The small ball assumption is a norm equivalence assumption between the $L_{1}$ and $L_{2}$ norms and is concerned with the second type of assumptions. Our approach deals with the first type of assumption since we only assume boundedness in $L_{2}$ which can be seen as an significant improvement upon a $L_{\infty}$ boundedness assumption even though the rates of convergence obtained for our estimators turn out to be minimax in the absence of a margin condition even under a $L_{\infty}$ assumption.

The estimation of the expectation of a univariate variable by median-of-means (MOM) is done as follows: given a partition of the dataset into blocks of the same size, the empirical mean is constructed on each block and the MOM estimator is obtained by taking the median of these empirical means (see Section 2.2 for details). These estimators are naturally resistant to the presence of a few outliers in the dataset: if the number of these outliers does not exceed half the number of blocks, more than half of these blocks are made of “clean” data and the median is a reliable estimator.

In this article, we also study a family of algorithms approaching the MOM minimizer. The general idea we wish to promote is to use the MOM principle within algorithms originally intended for the evaluation of ERM associated to convex loss functions. In Section 4, we present the modification of gradient descent algorithms following this philosophy. The general principle of this iterative algorithm is as follows: at iteration $t$, a dataset equipartition $B_{1},\ldots,B_{K}$ is selected uniformly at random and the most central block $B_{\text{med}}$ is determined according to the following formula

In Section 5, the practical performances of this algorithm are illustrated on several simulations, involving in particular different loss functions. These simulations illustrate not surprisingly the gain of robustness that there is to use these algorithms in their MOM version rather than in their traditional version, as can for example be appreciated on the toy-example of Figure 1, cf. Figure 5. MOM estimators are compared to different learning algorithms on real datasets that can be modeled by heavy tailed data, obtaining in each case performances comparable to the best of these benchmarks.

Another advantage of our procedure is that it works on blocks of data. This can improve speed of execution and reduce memory requirements, which can be decisive on massive datasets and/or when one wishes to use non-linear algorithms as in Section 4.3. This principle of dividing the dataset to calculate estimators more quickly and then aggregating them is a powerful tool in statistics and machine learning . Among others, one can mention bagging methods or subagging —a variant of bagging where the bootstrap is replaced by subsampling— . These methods are considered difficult to study theoretically in general and their analysis is often limited to obtaining asymptotic guarantees. By contrast, the theoretical tools for non-asymptotic risk analysis of MOM minimizer have already essentially been developed and the consistency of the algorithm can also be derived from standard arguments. Finally, subsampling by the central block $B_{{\rm med}}$ ensures robustness properties that cannot be guaranteed by traditional alternatives.

Detecting outliers is usually performed in machine learning via some unsupervised preprocessing algorithm that detects outliers outside a bulk of data, see for example or other algorithms like DBSCAN or isolation forest . These algorithms assume elliptical symmetry of the data, a solution for skewed data can also be found in . These unsupervised preprocessing removes outliers in advance, i.e. before starting any learning task. As expected, these strategies work well in the toy example from Figure 1. There are several cases where it will fail though. First, as explained in , this strategy classifies data independently of the risk, it is likely to remove from the dataset outlier coming from heavy-tailed distribution, yielding biased estimators. Moreover, a small group of misclassified data inside a bulk won’t be detected. Our notion of depth, based on the risk, seems more adapted to the learning task than any preprocessing procedure blind to the risk.

The paper is organized as follows. Section 2 presents the classification problem, the ERM and its MOM versions and gathers the assumptions granted for the main results. Section 3 presents theoretical risk bounds for the ERM and MOM minimizers on corrupted datasets. Section 4 deals with theoretical results on the algorithm computing MOM minimizers. We present the algorithm, study its convergence and provide theoretical bounds on its complexity. Section 5 shows empirical performance of our estimators in both simulated and real datasets. Proofs of the main results are postponed to Section 6 where we also added heuristics on the practical choice of the hyper-parameters.

Setting

Consider the supervised binary classification problem, where one observes a sample $(X_{1},Y_{1})$, $\ldots$, $(X_{N},Y_{N})$ taking values in ${\cal X}\times{\cal Y}$. The set ${\cal X}$ is a measurable space and $\mathcal{Y}=\{-1,1\}$. The goal is to build a classifier —that is, a measurable map $f:{\cal X}\to{\cal Y}$— such that, for any new observation $(X,Y)$, $f(X)$ is a good prediction for $Y$. For any classifier $f$, let

Note that $\widehat{f}_{\text{ERM}}$ does not build a classifier. To deduce, a classification rule from $\widehat{f}_{\text{ERM}}$ one can simply consider its sign function defined for all $x\in{\cal X}$ by ${\rm sign}(\widehat{f}_{\text{ERM}}(x))=2(I\{\widehat{f}_{\text{ERM}}(x)\geq 0\}-1/2)$. The procedure $\widehat{f}_{\text{ERM}}$ is solution of a convex optimization problem that can therefore be approximated using a descent algorithm. We refer for example to for a recent overview of this topic and Section 4 for more examples.

2 Corrupted datasets

In this paper, we consider a framework where the dataset may have been corrupted by outliers (or anomalies). There are several definitions of outliers in the literature, here, we assume that the dataset is divided into two parts. The first part is the set of inliers, indexed by ${\cal I}$, data $(X_{i},Y_{i})_{i\in{\cal I}}$ are hereafter always assumed to be independent and identically distributed (i.i.d.) with common distribution $P$. The second one is the set of outliers, indexed by ${\cal O}\subset[N]$ which has cardinality $|{\cal O}|$. Nothing is assumed on these data which may not be independent, have distributions $P_{i}$ totally different from $P$, with no moment at all, etc… In particular, this framework is sufficiently general to cover the case where outliers are i.i.d. with distribution $Q\neq P$ as in the $\epsilon$-contamination model .

3 Main assumptions

As already mentioned, data are divided into two groups, a subset $\{(X_{i},Y_{i}):i\in{\cal O}\}$ made of outliers (on which we will make no assumption) and the remaining data $\{(X_{i},Y_{i}):i\in{\cal I}\}$ contains all data that bring information on the target/oracle

Data indexed by ${\cal I}$ are therefore called inliers or informative data. To keep the presentation as simple as possible, inliers are assumed to be i.i.d. distributed according to $P$ although this assumption could be relaxed as in . Finally, note that the ${\cal O}\cup{\cal I}=\{1,\ldots,N\}$ partition of the dataset is of course unknown from the statistician.

Let us now turn to the set of assumptions we will use to study MOM minimizers procedures. For any measure $Q$ and any function $f$ for which it makes sense, denote by $Qf=\int fdQ$. Denote also, for all $q\geq 1$, by $L^{q}$ the set of real valued functions $f$ such that $\int|f|^{q}dP<\infty$ and, for any $f\in L^{q}$, by

Our first assumption is an $L^{2}$-assumption on the functions in $F$.

For all $f\in F$, we have $\left\|f\right\|_{L^{2}}\leq\theta_{2}$.

Of course, Assumption 1 is granted if $F$ is a set of classifiers. It also holds for the linear class of functions from Example 1 as long as $P\left\|X\right\|_{2}^{2}<\infty$ with $\theta_{2}=\Gamma(P\left\|X\right\|_{2}^{2})^{1/2}$.

The second assumption deals with the complexity of the class $F$. This complexity appears in the upper bound of the risk. It is defined using only informative data. Let

The Rademacher complexity is a standard measure of complexity in classification problems . It can be upper bounded, depending on the situation, by the VC dimension or the Dudley’s entropy integral or the Gaussian mean width of the class $F$, see for example for a presentation of these classical bounds. Our second assumption is simply that the Rademacher complexity of the class $F$ is finite.

The Rademacher complexity of $F$ is finite, ${\cal R}(F)<\infty$.

Theoretical guarantees

Our first result follows Vapnik’s original risk bound for the ERM and shows that $\widehat{f}_{\text{ERM}}^{0-1}$ is insensitive to the presence of outliers in the dataset. Moreover, it quantifies this robustness property since Vapnik’s rate of convergence is still achieved by $\widehat{f}_{\text{ERM}}^{0-1}$ when there are less than (number of observations) times (Vapnik’s rate of convergence) outliers.

Theorem 1 is proved in Section 6.1. It is an adaptation of Vapnik’s proof of the excess risk bounds satisfied by $\widehat{f}_{\text{ERM}}^{0-1}$ in the presence of outliers.

In the last result, one can easily bound the excess risk using ${\cal R}(\mathcal{F})$ instead of ${\cal R}({\cal L}_{\cal F}^{0-1})$ since

The final bound is of similar flavor: for all $K>0$, with probability at least $1-e^{-K}$, we have

As a consequence, when the number of outliers is such that $|{\cal O}|\lesssim N\times\sqrt{p/N}$, Vapnik’s classical “slow” rate of convergence $\sqrt{p/N}$ is still achieved by the ERM even if $|{\cal O}|$ outliers have polluted the dataset. The interested reader can also check that “fast rates” $p/N$ could also be achieved by the ERM in the presence of outliers if $|{\cal O}|\lesssim p$ and when the so-called strong margin assumption holds (see, ). Note also that the previous remark also holds if $F$ is a class with VC dimension $p$ beyond the case of indicators of half spaces.

satisfies an excess risk bound under weak assumptions introduced in Section 2.

Grant Assumptions 1, 2 and 3. Assume that $N>K>4|{\cal O}|$ and let $\Delta=1/4-|{\cal O}|/K$. Then, with probability larger than $1-2\exp\left(-2\Delta^{2}K\right),$ we have

Theorem 2 is proved in Section 6.2. Compared to Theorem 1, $\hat{f}_{\text{MOM},K}$ achieves the same rate $({\cal R}(F)/N)\vee(\sqrt{K/N})$ under the same conditions on the number of outliers with the same exponential control of the probability as for the ERM $f_{\text{ERM}}^{0-1}$. The main difference is that the loss function may be unbounded, which is often the case in practice. Moreover, unlike classical analysis of ERM obtained by minimizing an empirical risk associated with a convex surrogate loss function, we only need a second moment assumption on the class $F$.

These theoretical improvements have already been noticed in previous works . Contrary to tournaments of , Le Cam MOM estimators of or minmax MOM estimators , Theorem 2 does not require the small ball assumption on $F$ but only shows “slow rates” of convergence. These slow rates are minimax optimal in the absence of a margin or Bernstein assumption . Proof of Theorem 2 does not enable fast rates to be obtained. Indeed, the non-linearity of the median excludes the possibility of using localization techniques leading to these fast rates. However, we show in the simulation study (cf. left side picture of Figure 14) that fast rates seem to be reached by the MOM minimizer.

The main advantage of our approach is its simplicity, we just have to replace empirical means by their MOM alternative in the definition of the ERM. Moreover, as expected, this simple alternative to ERM yields a systematically way to modify algorithms designed for approximating the ERM. The resulting “MOM version” of these algorithms are more robustness and faster than their original “ERM version”. Before illustrating these facts on simulations, let us describe algorithms approximating MOM minimizers.

Computation of MOM minimizers

In this section, we present first a generic algorithm to provide a MOM version of descent algorithms minimizing the empirical risk. We show how classical convergence proofs can be adapted to these algorithms.

The general idea is that any descent algorithms such as gradient descent, Newton method, alternate gradient descent, etc. (cf. ) can easily be turned into a robust MOM-version. To illustrate this idea, a basic gradient descent is analyzed in the sequel.

The choice of blocks greatly influences the practical performance of the algorithm. In particular, a recurring flaw is that iterations tend to get stuck in local minima, which greatly slows the convergence of the alogorithm. To overcome this default and improve the stability of the procedure, we modify the partition at each iteration by drawing it uniformly at random, cf. step 2 of Algorithm 1.

Let ${\cal S}_{N}$ denote the set of permutations of $\{1,\ldots,N\}$. For each $\sigma\in{\cal S}_{N}$, let $B_{0}(\sigma)\cup\cdots\cup B_{K-1}(\sigma)=\{1,\ldots,N\}$ denote an equipartition of $\{1,\ldots,N\}$ defined for all j\in\left\llbracket 0,K-1\right\rrbracket by

In Section 5, we use MOM ideas (as in the generic Algorithm 1) to construct MOM versions for various classical algorithms such as Perceptron, Logistic Regression, Kernel Logistic Regression, SGD Classifiers or Multi-layer Perceptron.

2 Convergence properties

We denote by ${\cal D}_{N}=\{(X_{1},Y_{1}),\dots,(X_{N},Y_{N})\}$ the dataset. We consider the following set of assumptions:

There exist $L>0$ and a sequence of decreasing positive numbers $(\eta_{t})_{t\geqslant 1}$, such that

$\sum_{t=1}^{\infty}\eta_{t}^{2}<\infty$,

Recall that the dataset ${\cal D}_{N}=\{(X_{i},Y_{i}):i=1,\ldots,N\}$ may contain outlier, which might be deterministic. We use the terminology “for almost all datasets ${\cal D}_{N}=\{(X_{i},Y_{i}):i=1,\ldots,N\}$ such assumption holds” when the latter holds for $P^{|{\cal I}|}$-almost all informative data $(X_{i},Y_{i})_{i\in{\cal I}}$ and for every outliers $(X_{i},Y_{i})_{i\in{\cal O}}\in({\cal X}\times{\cal Y})^{|{\cal O}|}$.

(where the a.s. convergence is with respect to the randomness on the choice of data partitions at each step).

[MOM gradient descent, SGD and assumption H-6] All assumptions in Theorem 3 are the classical assumption used to prove the convergence of the classical SGD except for H-6 which is specific to MOM algorithms.

The assumption H-5 illustrate why we used a random shuffle of the data at each step of the algorithm, the principle is that the MOM risk in general may present local minima and hence the algorithm without the shuffle would stop in a local minimum. By using the random shuffle at each step, it seems that the problem has a unique minimum, at least experiment seems to say so. For example, in figure 3 we see the plot of the loss function when we change one of the parameters (the intercept) of a linear classifier for the classification problem depicted in figure 4.

where for all blocks $B\subset\{1,\ldots,N\},P_{B}f_{p}=|B|^{-1}\sum_{i\in B}f_{p}(x_{i})$ and the blocks $B(k)(p),k=1,\ldots,K$ are rearranged blocks defined such that $P_{B(1)(p)}(f_{p})\geq\dots\geq P_{B(K)(p)}(f_{p})$. Proposition 1 shows that Assumption H-6 is satisfied in several situations.

where $\partial_{j}$ denote the derivative with respect to the $j^{th}$ coordinate of $p$.

Proposition 1 is proved in Section 6. The function $h$ is bounded for a large panel of probability distributions. For example, an exponential variable with parameter $\lambda$ has a constant hazard function equal to $\lambda$; a Cauchy random variable has also a bounded hazard function. Intuitively, this property of derivation of the median of means is more easily verified when the distribution of $P_{B(k)(p)}(f)$ is heavy tailed (i.e. the mean blocks are far from one another) and in fact we can verify that the hazard function of the Gaussian is unbounded and hence this proposition cannot be used on Gaussian random variables.

3 Complexity of MOM risk minimization algorithms

Let $C(m)$ be the computational complexity of a single ERM gradient descent update step on a dataset of size $m$ and let $L(m)$ be the computational complexity of the evaluation of the loss function on a dataset containing $m$ data. Here the computational complexity is simply the number of basic operations needed to perform a task .

For each epoch, we begin by computing the “MOM empirical risk”. We perform $K$ times $N/K$ evaluations of the loss function, then we sort the $K$ means of these blocks of loss to finally get the median. The complexity of this step is then $O(KL(N/K)+K\ln(K))$, assuming that the sort algorithm is in $O(K\ln(K))$ (like quick sort ). Then we do the ERM gradient step on a sample of size $N/K$. Hence, the time complexity of this algorithm is

For example, if the ERM gradient step and the loss function evaluation have linear complexity – like Perceptron or Logistic Regression – the complexity of the MOM algorithm is $O(T(N+K\log(K)))$ against $O(TN)$ for the ERM algorithm. Therefore, the two complexities are of the same order and the only advantage of MOM algorithms lies in their robustness to outliers.

If on the other hand the complexity is more than linear as for Kernel Logistic Regression (KLR), taking into account the matrix multiplications whose complexity can be found in , the complexity of the MOM version of KLR, due to the additional need of the computation of the kernel matrix, is $O(N^{2}+T(N^{2}/K+K\log(K)+(N/K)^{2.373}))$ against $O(TN^{2.373})$ for the ERM version. MOM versions of KLR are therefore faster than the classical version of KLR on top of being more robust. This advantage comes from the fact that MOM algorithms work on blocks of data instead on the entire dataset at every step. More informations about Kernel Logistic Regression can be found in for example.

In this last example, the complexity comes in part from the evaluation of the kernel matrix that can be computationally expensive. Following the idea that MOM algorithms are performing ERM algorithm restricted to a wisely chosen block of data, then one can modify our generic strategy in this particular example to reduce drastically its complexity. The idea here is that we only need to construct the kernel matrix on the median block. The resulting algorithm, called Fast KLR MOM is described in Figure 2.

In Figure 2, we compute only the block kernel matrices, denoted by $N^{1},\dots,N^{k}$ and constructed from the samples in the block $B_{k}$. We also denote by $N_{i}^{k}$ the $i^{th}$ row in $N^{k}$.

There are several drawbacks in the approach of Algorithm 2. First, the blocks are fixed at the beginning of the algorithm; therefore the algorithm needs a bigger dataset to work well and it may converge to a local minimum. Nonetheless, from the complexity point of view, this algorithm will be much faster than both the classical KLR and MOM KLR (see below for a computation of its complexity) which is important given the growing use of kernel methods on very large databases for example in biology. The choice of $K$ should ultimately realize a trade-off between complexity and performance (in term of accuracy for example) when dealing with big databases containing few outliers.

The complexity of Fast KLR-MOM is $O(N^{2}/K+T(N^{2}/K+K\log(K)+(N/K)^{2.373}))$ against $O(TN^{2.373})$ for the ERM version.

Implementation and Simulations

The algorithms we study are the MOM adaptations of Perceptron, Logistic Regression and Kernel Logistic Regression.

Based on our theoretical results, we know that the number of blocks $K$ has to be larger than $4$ times the number of outliers for our procedure to be on the safe side. The value $K=120$ is therefore used in all subsequent applications of MOM algorithms on the toy dataset except when told otherwise. To quantify performance, we compute the miss-classification error on a clean dataset made of data distributed like the informative data.

For Kernel Logistic Regression, we study here a linear kernel because outliers in this dataset are clearly adversarial when dealing with linear classifiers. The algorithm can also use more sophisticated kernels, a comparison of the MOM algorithms with similar ERM algorithms is represented in figure 5, the ERM algorithms are taken from the python library scikit-learn with their default parameters.

Figure 5 illustrates resistance to outliers of MOM’s algorithms compared to their classical version.

These first results are completed in Figure 6 where we computed accuracy on several run of the algorithms. These results confirm the visual impression of our first experiment.

Finally, we illustrate our results regarding complexities of the algorithms on a simulated example. MOM algorithms have been computed together with state-of-the art algorithms from scikit-learn (we use Random forest, SVM classifier as well as SGD classifier optimizing Huber loss which entail a robustness in $Y$ but not in $X$, see [34, Chapter 7]) on a simulated dataset composed of two Gaussian blobs $\mathcal{N}((-1,-1),1.4I_{2})$ and $\mathcal{N}((1,1),1.4I_{2})$ with label respectively $1$ and $-1$. We sample $20000$ points for the training dataset and $20000$ for the test dataset. The parameters used in the algorithms are those for which we obtained the optimal accuracy, (this accuracy is illustrated in the next section). Time of training plus time of evaluation on the test dataset are gathered in Figure 7.

Not surprisingly, very efficient versions of linear algorithms from Python’s library are extremely fast (results are sometimes provided before we even charged the dataset in some experiments). The performance of our algorithm are nevertheless acceptable in general (around 5 times longer than random forest for example). The important fact here is that non linear algorithms such as SVM take much more time to provide a result. FAST KLR MOM is able to reduce substantially the execution time of SVM with comparable predictive performance.

2 Applications on real datasets

We used the HTRU2 dataset, also studied in , that is provided by the UCI Machine Learning Repository. The goal is to detect pulsars (a rare type of Neutron star) based on radio emission detectable on earth from which features are extracted to gives us this dataset. The problem is that most of the signal comes from noise and not pulsar, the goal is then to classify pulsar against noise, using the $17$ $898$ points in the dataset.

The accuracy of different algorithms is obtained using on several runs of the algorithms each using 4/5 of the datasets for training and $1/5$ for testing algorithms. Boxplots presenting performance of various algorithms are displayed in Figure 8. To improve performance, RBF kernel was used both for KLR MOM and Fast KLR MOM.

3 Outlier detection with MOM algorithms

We apply this idea on the toy dataset with the Logistic Regression MOM algorithm. Results are gathered in a sorted histogram given in Figure 9. Red bars represent outliers in the original datasets.

Quite remarkably, outliers are in fact those data that have been used the smallest number of times. The method targets a very specific type of outliers, those disturbing the classification task at hand. If there was a point very far away from the bulk of data but in the half-space of its label, it wouldn’t be detected.

This detection algorithm doesn’t scale well when the dataset gets bigger as a large number of iterations is necessary to choose each point a fair number of times. For bigger datasets, we suggest to adapt usual outlier detection algorithms . We emphasize that clustering techniques and K-Means are rather easy to adapt in a MOM algorithm and detect points far from the bulk of data. This technique might greatly improve usual K-Means as MOM K-Means is more robust.

Let us now analyze the effect of $K$ on the outlier detection task. The histogram of the $1000$ smaller counts of points of HTRU2 dataset as $K$ gets bigger is plotted in Figure 10.

It appears from Figure 10 that $K$ measures the sensitivity of the algorithm. Severe outliers (as in the toy example) are detected for small $K$ while mild outliers are only discovered as $K$ gets bigger.

It seems therefore that the optimal choice of $K$ in MOM depends on the task one is interested in. For classification, $K$ should be as small as possible to get better risk bounds (but it still should be larger than the number of outliers) whereas for detecting outliers we may want to choose $K$ much larger to even detect an outlier, (but it should also be small enough for the underlying classification to perform correctly). As a proof of concept, for Pulsar database, we got optimal results choosing $K=10$ for classification whereas we only detect a significant amount of outliers when $K$ is around $1000$.

Proofs

It follows from the definition of the ERM that $P_{N}{\cal L}_{\hat{f}}\leq 0$. Therefore, if we denote by $P_{\cal I}$ (resp. $P_{\cal O}$) the empirical measure supported on $\{(X_{i},Y_{i}):i\in{\cal I}\}$ (resp. $\{(X_{i},Y_{i}):i\in{\cal O}\}$), we have

because $|{\cal L}_{\hat{f}}|\leq 1$ a.s.. Then, by the bounded difference inequality [10, Theorem 6.2], since all $f\in{\cal F}$ satisfies $-1\leq{\cal L}_{f}\leq 1$, one has, for any $x>0$,

Furthermore, by the symmetrization argument (cf. Chapter 4 in ),

Therefore, for any $x>0$, with probability larger than $1-e^{-2|{\cal I}|x^{2}}$,

The proof is completed by choosing $x=\sqrt{K/(2|{\cal I}|)}$.

2 Proof of Theorem 2

Let us now control the two expressions in the right-hand side of (7). Let $x>0$. We have

where we recall that ${\cal K}=\{k\in\{1,\cdots,K\}:B_{k}\cap{\cal O}=\emptyset\}\}$.

As $0\leq\phi(\cdot)\leq 1$, by the bounded-difference inequality, for any $y>0$, with probability larger than $1-e^{-2y^{2}K}$,

Since $\phi$ is $1$-Lipschitz and $\phi(0)=0$, by the contraction principle (see [45, Chapter 4] or more precisely equation (2.1) in ),

Thus, for any $y>0$, with probability larger than $1-\exp(-2y^{2}K)$,

Let $\Delta=1/4-|{\cal O}|/K$ and let $y=\Delta$ and $x=8L\max\left(\theta_{2}\sqrt{K/N},4{\cal R}(F)/N\right)$ so

Plugging this result in (7) concludes the proof of the theorem.

3 Proof of Theorem 3

From here, the rest of the proof is identical to the proof of consistency of the SGD. As a proof of concept, we reproduce here the arguments of . In what follows, all results are given conditionally on ${\cal D}_{N}$. We set

By definition of the update of the algorithm, we have

Let $F_{t}$ be the $\sigma$-algebra spanned by $\sigma_{1},\dots,\sigma_{t-2}$, we have

Then, we use Jensen’s inequality, we have

Hence, using hypothesis H-1 and H-5, we get

We plug this convergence in the inequality (10) to get

From the hypotheses H-4 and H-5, we can then conclude that

4 Proof of Proposition 1

We denote by $B(1)(f),\dots,B(K)(f)$ the blocks such that the corresponding block empirical mean $P_{B(i)(f)}(f(X_{1}^{n}))$ are sorted: $P_{B(1)(f)}(f(X_{1}^{n}))\geq\dots\geq P_{B(K)(f)}(f(X_{1}^{n}))$. To simplify the notations we denote by $K/2=\lceil K/2\rceil$ that correspond to the rank of the median when $K$ is odd, we suppose $K$ odd.

To show that, it is sufficient to check that we have

We decompose this difference in three parts,

The two last terms are controlled by the Lipshitz property of $p\mapsto f_{p}$,

Then, we use order statistics and interquantile range to control the first term.

To do that, we use Renyi representation to say that the order statistics of the blocks can be expressed as functions of the order statistics of exponential variables. Let $E_{1},\dots,E_{K}$ be K i.i.d sample with a standard Exponential law, then using the function $U(t)=F_{p}^{\leftarrow}(1-1/t)$ where $F_{p}$ is the c.d.f of the random variable $P_{B_{k}}(f_{p})$ to make the link between the sample $P_{B(1)(p)},\dots,P_{B(K)(p)}$ and the sample $E_{1},\dots,E_{K}$, we have

Then, because the hazard function is the inverse of the derivative of $U$ and is bounded by $1/m_{N,K}$, we have by the mean value inequality,

which we can rewrite as: for all $\delta\in(0,1)$, with probability greater than $\delta$,

Finally, with probability greater than $\delta$,

Which implies that with probability greater than $\exp(-ktL||\varepsilon||_{2}/(2m_{N,K}))$,

Then, taking the union bound for $1\leq k\leq K-1$, we have

and we study the limiting event $\Omega^{j}=\overline{\lim}_{n\to\infty}A_{n}^{j}$.

First, let us note that for all $1\leq j\leq d$, the sequence of set $(A_{n}^{j})_{n}$ is non-increasing, hence

then, for all $1\leq j\leq d$, we can study the $\overline{\lim}_{n\to\infty}(A_{n}^{j})^{c}$ with Borel-Cantelli Lemma. Indeed, we have

which implies that for all $1\leq j\leq d$, $\forall n\geq N$,

Annex

Let us study the behaviour of our algorithms when the number of blocks changes. We plot the accuracy as a function of $K$ averaged on 50 runs to have a good idea of the evolution of the performance with respect to $K$, the result is represented in figure 11.

There is a clear separation around $2|{\cal O}|=60$ that is consistent with the theory. On the other hand the accuracy doesn’t decrease when $K$ gets bigger one would expect. This may be due to the symmetry of the dataset. If we run the same experiment on the real dataset, we get a much more regular plot, see Figure 12.

Figure 12 confirms our predictions on clean datasets, the accuracy getting better as $K$ gets smaller (the MOM minimizer is the ERM when $K=1$ and ERM is optimal in the i.i.d. setup, ). This may be due to the small number of outliers in this dataset.

2 Illustration of convergence rate

In this section, we estimate the rate of convergence of the MOM risk minimization algorithm Logistic Regression on two databases (see figure 13). The first dataset is composed of points located on two interlaced half-circle with a Gaussian noise of standard deviation $0.3$, the two ”moons” are each of a different class. We assume that these moons don’t satisfy the margin property (we checked that the rate was slow for ERM algorithms, using the vanilla logistic regression). The second dataset is composed of two Gaussians $\mathcal{N}((-1,-1),1.4^{2}I_{2})$ and $\mathcal{N}((1,1),1.4^{2}I_{2})$ with respective label $1$ and , we can prove that this dataset verifies the margin property needed to obtain fast rate in ERM

There are no outliers in the datasets because we only want to test the rate of convergence. To illustrate the rates of convergence of our algorithms, we plot the curve $\log\left(\left|\hat{R}^{0-1}(\hat{f}_{K}))-\hat{R}^{0-1}(f^{*})\right|\right)$ as a function of $\log(n)$ where the risk is estimated by Monte-Carlo. The figure obtained for Logistic Regression MOM is represented in figure 14. It seems that MOM minimizers can achieve fast rates of convergence even if we did not prove them.

We used random blocks sampled at each iteration for this application because it is the algorithm that we described earlier but even if we use one partition of blocks for the whole algorithm (as in the theory we developed) we obtain nonetheless fast rate for the Gaussians dataset.

3 Comparison with robust algorithms based on M-estimators.

In this section we compare the algorithm Logistic Regression MOM with two other algorithms based on M-estimators, these algorithms are studied on the toy dataset presented in Section 5.

where $\psi_{c}=\max(-c,\min(c,x))$ is the Huber function, $c>0$. Using this definition of $\hat{\mu}_{f}$, it is then easy to compute the gradient $\nabla\hat{\mu}_{f}$ and then use a gradient descent algorithm. The theory behind this algorithm is studied further in .

The second algorithm uses a ”redescending” loss function, in short we do ERM with a bounded loss function. Here we use Tukey biweight loss function rescaled by MADN scale estimator and IRLS algorithm to optimize the empirical risk.

Figure 15 shows that all algorithms perform similarly on this easy, low dimensional dataset. The situation is quite different in higher dimension. In Figure 16 we used a 200 dimensional dataset and the algorithm using a redescending loss function does not perform well. This may be due to local minima in which the algorithm gets stuck, as local minima are multiplied when the dimension gets higher. The other algorithms don’t suffer this drawback since they use a “projection by the loss function” that makes the problem one dimensional.

The algorithm using redescending loss functions is a simple gradient descent that has linear complexity. The Huber gradient algorithm estimates at each iteration a Huber estimator of location. The complexity of this estimator depends on the algorithm used but for most M-estimators a commonly used algorithm is an iteratively reweighted algorithm whose complexity is linear in the sample size. In practice we can nonetheless notice a great complexity of the Huber estimator in some cases where data are not well spread. In most cases, Logistic Regression MOM is the fastest among these three algorithms and the gradient Huber is the slowest, even though logistic regression may need a lot more iterations than the other algorithms.

References