Robust regression via error tolerance

Real-world datasets are often characterised by outliers; data items that do not follow the same structure as the rest of the data. These outliers might negatively influence modelling of the data. In data analysis it is, therefore, important to consider methods that are robust to outliers. In this paper we develop a robust regression method that finds the largest subset of data items that can be approximated using a sparse linear model to a given precision. We show that this can yield the best possible robustness to outliers. However, this problem is NP-hard and to solve it we present an efficient approximation algorithm, termed SLISE. Our method extends existing state-of-the-art robust regression methods, especially in terms of speed on high-dimensional datasets. We demonstrate our method by applying it to both synthetic and real-world regression problems.


Introduction
In practically all analyses of real-world data we encounter outliers, i.e., data items that do not follow the same patterns as the majority of the data. Such items are problematic, since they may negatively influence modelling of the data. This is observed, for instance, in ordinary least-squares (ols) regression where already a single outlier may lead to arbitrarily large errors (Donoho and Huber 1983). It is, hence, important to consider robust methods that effectively avoid the influence of outliers.
Robust regression methods can be used as almost drop-in replacements for linear regression, which is still widely used because of the inherent interpretability and simplicity. Linear regression is also often used as a part of other machine learning or data mining algorithms, e.g., in explainable artificial intelligence (Ribeiro et al. 2016;Björklund et al. 2019). Furthermore, robust regression can be used to search for outliers by investigating the data items that do not adhere to the robust model.
In this paper we present a sparse robust regression method, termed slise (Sparse Linear Subset Explanations), that achieves the highest possible theoretical robustness and outperforms many existing state-of-the-art robust regression methods in terms of scalability on large datasets. Specifically, we consider finding the largest subset of data items that can be represented by a linear model to a given accuracy.
To illustrate the need for robust regression methods, consider the dataset shown in Fig. 1 containing a cluster of outliers in the top right corner. Here ordinary leastsquares linear (ols) regression clearly fails due to the influence of the outliers. In contrast, our slise method yields high robustness by ignoring the outliers and finds a linear model that approximates a (large) subset of the data items.
The two main goals of the slise algorithm can be stated in brief as follows. The primary goal of the algorithm is to find a maximum subset of points that can be modelled by a linear model to a given maximum error. This subset is visualised in Fig. 1 by the shaded area around the orange slise line. The width of this "corridor" shows the (adjustable) error tolerance. The secondary goal of slise is to minimise the residuals of the items within the subset marked by the corridor. Fig. 1 Outliers can cause ordinary least-squares regression to give unusable results (green dashed line). Robust regression mitigates the effect of the outliers. slise (orange solid line) accomplishes this by finding a maximal subset of items which can be modelled by a linear model to a given maximum error. This maximum error is illustrated as the lightly shaded "corridor" around the slise regression line (Color figure online) models for a subset of the items. The difference is that lts requires the size of subset to be fixed and specified a priori while slise defines the subset based on a maximum tolerable error. This requires a new algorithm that actually scales better with regards to the number of dimension (see Sect. 4.2) without sacrificing any theoretical robustness (see Sect. 2.2). The slise algorithm uses graduated optimisation to find solutions, which makes it both fast and robust against noise and outliers.

Contributions
We present a novel robust regression method, slise, by considering the problem of finding the largest subset that can be approximated by a sparse linear model to a given accuracy (Problem 1). We show that this can yield the best possible breakdown value (Sect. 2.2), but that the problem is NP-hard (Theorem 2). We present an approximative algorithm for solving it (Algorithm 1) and demonstrate that it also works empirically using both synthetic and real-world datasets (Sects. 5-7). Compared to other robust regression methods slise yields high robustness and good scalability to high-dimensional datasets.
The initial version of the slise algorithm was presented in aconference paper (Björklund et al. 2019). That paper has a greater focus on one specific use case for the robust regression method, namely, explaining outcomes from black box supervised learning models. (The word "Explanations" in the name of the slise algorithm stems from this application.) When applied to the problem of explaining models, the idea is to find a simple interpretable linear model that approximates a more complex supervised learning model in the neighbourhood of a data item of interest. For this task slise must be able to find good solutions even for small error tolerances. The advantage of using slise for explaining supervised learning models is that no resampling of datasets is required and that the explanations (linear models) obey constraints imposed by the data. In the present paper, however, we focus solely on robust regression.
Compared to Björklund et al. (2019), the discussion on the problem characteristics and numeric approximation is substantially extended with additional proofs in Sects. 2 and 3. We also provide novel initialisation schemes for the slise algorithm in Sect. 4, and evaluate their effect, and the effect of all the other parameters, on the stability and performance of the algorithm in Sect. 6. We find that one of the new schemes is more robust than the lasso initialisation used previously, and provide recommendations for suitable default values for the other parameters. Furthermore, we perform a more thorough empirical comparison to related methods in Sect. 7.

Organisation
In Sect. 2 we formalise our robust regression problem, and show its complexity and breakdown value. We then discuss the practical numeric optimisation and its approximation properties in Sect. 3. The algorithm, with different initialisation schemes and asymptotic complexity, is presented in Sect. 4. The empirical evaluation, which consists of the default parameter selection and the comparison to related methods, follows in Sects. 5-7. We end with conclusions in Sect. 8.

Problem definition
Let (X , Y ), where X ∈ R n×d and Y ∈ R n , be a dataset consisting of n pairs where we denote the ith d-dimensional item (row) in X by x i ∈ R d (the predictor) and similarly the ith element in Y by y i ∈ R (the response). We use the shorthand [n] = {1, . . . , n}.
Our goal is to develop a linear regression method that is robust to outliers by finding the largest subset of data items that can be described by a sparse linear model to a given precision, as exemplified in the introduction.
We now state the main problem in this paper: Problem 1 Given X ∈ R n×d , Y ∈ R n , the error tolerance ε ∈ R ≥ 0, and the regularisation strength λ ∈ R ≥ 0; find the regression coefficients α ∈ R d minimising the (ε, λ)-loss function where the residual errors are given by is the Heaviside step function satisfying H (u) = 1 if u ≥ 0 and H (u) = 0 otherwise, and α 1 denotes the L1-norm. In the remainder of the paper we will use the L1-norm given by α 1 = d i=1 |α i |, even though the theoretical results would be valid for L2 or other norms as well. The benefits of L1 or Lasso regularization include that it leads to sparse solutions, which is beneficial for interpretability and explanations (Björklund et al. 2019). If necessary, the data matrix X can be augmented with a column of all ones to accommodate the intercept term of the model. Alternatively, the Lagrangian term λ α 1 in Eq. (1) can be replaced by a constraint α 1 ≤ t for some t. Now we can rewrite Eq. (1) as under the constraint α 1 ≤ t, and where Here S is the subset of data items assumed to be non-outliers and the complement, S c = [n] \ S, is the subset of potential outliers. As can be seen in Eq. (2), we only consider the non-outliers when fitting the linear model.
At the limit of ε → ∞ it follows that S = [n] and Problem 1 becomes equivalent to lasso (Tibshirani 1996). When ε is small Problem 1 becomes a combinatorial problem in disguise, where we try to find a maximal subset S, due to the subtraction of ε 2 in Eq. (2).

Theorem 1 Minimising the loss in Eq.
(2) maximises the size of the subset S in Eq. (3).

Proof
The subset size term in Eq. (2) is − i∈S ε 2 = −|S|ε 2 , while for the residual term it holds that i∈S r 2 i /n ≤ ε 2 . This means that any change in the subset size has at least as large an impact on the loss as all the residuals combined, and for any α and α * satisfying Loss(ε, λ, X , Y , α) ≤ Loss(ε, λ, X , Y , α * ), it then follows that |S| ≥ |S * |.

Complexity class
Due to the combinatorial nature, finding an exact solution to Problem 1 is difficult. By showing that Problem 1 is a generalisation of a known NP-hard problem we can give a lower bound for the complexity class.
Theorem 2 Problem 1 is NP-hard and hard to approximate.
Proof We prove the theorem by a reduction to the maximum satisfying linear subsystem problem (Ausiello et al. 1999, Problem MP10), which is known to be NP-hard . In maximum satisfying linear subsystem we are given the system X α = y, where X ∈ Z n×m and y ∈ Z n and we want to find α ∈ Q m such that as many equations as possible are satisfied. This is equivalent to Problem 1 with ε = 0 and λ = 0. Additionally, the problem is not approximable within n γ for some γ > 0, according to Amaldi and Kann (1995).

Breakdown value
The robustness of robust regression methods is often measured using the breakdown value (Donoho and Huber 1983;Rousseeuw and Hubert 2011), which is defined as the theoretical minimum fraction of (adversarial) outliers that can cause an arbitrarily large deviation in the model. This can, e.g., be measured with (Alfons et al. 2013): where v is the fraction of outliers and α v is the linear model that fits the dataset (X v , Y v ) where v of the items have been replaced by items with arbitrary values (outliers).
Non-robust regression methods, such as ordinary least-squares, have a breakdown value of 1/n (Hubert and Debruyne 2009), i.e., a single outlier is enough to cause a breakdown. The practical upper limit for the breakdown value is 0.5, since any value larger than that cannot be guaranteed, without prior information.
Theorem 3 The breakdown value of Problem 1 is 0.5.
Proof Following the definition, the breakdown value can be found as follows. We start from an uncorrupted (i.e., no outliers) dataset (X 0 , Y 0 ) with n items obeying a linear model parametrised by α 0 and where all data items are within the corridor defined by ε. A fraction v of the data items are then changed into adversarial outliers, yielding the corrupted dataset (X v , Y v ).
The subset given by Problem 1 on (X 0 , Y 0 ) is S 0 and |S 0 | = n. With a finite ε the breakdown occurs when the subset S v for (X v , Y v ) contains corrupted items. This requires |S v | ≥ |S 0 |/2, since otherwise the uncorrupted subset is larger, |S 0 \ S v | > |S v |, and is selected. Thus, the breakdown value is v = |S v |/n = (|S 0 |/2)/n = 0.5.
In Sect. 7.2 we empirically validate the breakdown value.

Numeric approximation
In order to solve the NP-hard Problem 1 efficiently in the general case, we relax the problem by replacing the Heaviside function with a sigmoid function σ (u) = 1/(1 + e −u ) and a continuous and differentiable rectifier function φ(u) ≈ min (0, u). This allows us to compute the gradient and find α by minimising where the parameter β determines the steepness of the sigmoid and the rectifier function φ is parametrised by a small constant ω > 0 such that It is easy to see that Eq. (4) is a smoothed variant of Eq.

Graduated optimisation
We perform the optimisation of Eq. (4) using graduated optimisation (Mobahi and Fisher 2015). Graduated optimisation iteratively solves a difficult optimisation problem by progressively increasing the complexity. A natural parametrisation for the complexity of our problem is via the β parameter, since β = 0 corresponds to a convex optimisation problem equivalent to lasso and when β → ∞ the problem becomes equivalent to Problem 1. This yields an iterative optimisation strategy. At each step we use the previous optimal value of α as a starting point for minimisation of Eq. (4), and then increase the value of β. It is important that the optima of consecutive α:s are "close enough". This is why we derive an approximation ratio between solutions with different values of β. We observe that our problem can be rewritten as a maximisation of −β-Loss(ε, λ, X , Y , α). The choice of β does not affect the L1-norm and we omit it for simplicity (λ = 0).
, and the inequality from the theorem holds.
We use Theorem 4 in the graduated optimisation to choose the sequence of increasing β values (β 1 , β 2 , . . . , β i > β i−1 ), so that the approximation ratio, defined as K , stays constant.

Stopping criterion
Iterating until β → ∞ is not possible in practice, so we need a stopping criterion for the algorithm. The iterations should be stopped when σ (β(ε 2 − r 2 )) ≈ H (ε 2 − r 2 ), i.e., the stopping criterion is dependent on the shape of the sigmoid function. The sigmoid shape is determined by both β and ε. However, ε is expected to change often, depending on both the dataset and the task, so a stopping criterion independent of ε would be preferable.
Proof Assume that there exists a pair of values β c and ε c . We say that a sigmoid function parametrised by some β max and ε has the same relative shape if and only if for every value of p ∈ R. Since the sigmoid function is strictly increasing, it can trivially be removed from both sides of the equation and hence where c = β c ε 2 c is a constant. We empirically determine a good default value for c in Sect. 6.5.

Approximation ratio for the numeric approximation
By using the approximation ratio between two β-Losses (Theorem 4) we can derive a new approximation ratio between the losses of Eqs. (4) and (1) (the problem definition and the numeric approximation).
We set β 2 → ∞ and ω → 0 . This leads to a new approximation ratio K ε * derived in the following lemma.

Lemma 1 The approximation ratio between the losses parametrised by
Proof The proof is omitted since it exactly mirrors Theorem 4 with the observation We call α * 2 the matching solution, since it is the optimum for Eq. (1) closest to α 1 . Note that α * 2 has a different ε (namely ε * ) that we need to specify in order to fully define the matching solution.
Due to the non-continuity of the Heaviside function, the maximum can be found at ε * = r j for some j ∈ [n]. We can further assume, without loss of generality, that the residuals are sorted so that r 2 1 ≤ r 2 2 ≤ · · · ≤ r 2 n , which means that n i=1 H (r 2 j −r 2 i ) = j. With a large enough n, so that 1/n ≈ 0, Eq. (7) can be simplified to Now, if the data is subsampled to a constant size, then Eq. (4) has a constant approximation ratio for the matching solution.

Theorem 5 The matching solution
and we can derive, by rearranging the terms in the inequality, Assuming that the residuals are sorted so that r 2 1 ≤ r 2 2 ≤ · · · ≤ r 2 n , then and further following the definition in Theorem 1 which, by rearranging, lets us derive Inserting this into G 1 yields and when combined with K ε * from Theorem 1 we have which completes the proof.

The SLISE algorithm
In this section we describe an approximate numeric algorithm, slise, for solving Problem 1. We start by introducing the algorithm, and then continue by discussing different initialisation schemes. Finally we demonstrate the asymptotic complexity of the slise algorithm. The slise algorithm (Algorithm 1) takes as input the data and the optimisation parameters. The algorithm starts by selecting initial values for the linear model α and the sigmoid steepness β (line 3). There are several potential initialisation schemes, and we will present and discuss four alternatives later in this section (see Sect. 4.1 and Algorithm 3).
The main part of slise consists of graduated optimisation (lines 4-7), that optimise the values for α and β. This is done by alternating between optimising α, and increasing β until we reach β max . At each iteration, we need to find the α that minimises the β-Loss in Eq. (4), using the current value of β (line 5). This optimisation is done with owl-qn (Schmidt et al. 2009) which is a quasi-Newton optimisation method with built-in L1-regularisation. We then increase β (line 6) such that the approximation ratio between consecutive steps, as defined in Theorem 4, equals K .
The pseudocode for the approximation ratio calculation is provided in Algorithm 2. In Eq. (4) we use the rectifying φ function to ensure negativity. This function requires a constant ω and its value can be chosen to be smaller than machine epsilon. Hence, in the approximation ratio calculation (line 4), φ(u) is numerically equivalent to min(0, u). As a minor computational side-effect of this, we have to make sure not to divide by zero if all φ(r i ) = 0 (lines 8-11).

Initialisation schemes
In Algorithm 3 we introduce four alternative initialisation schemes, and later in Sect. 6.3 we empirically compare the proposed approaches.
The first initialisation scheme (lines 2-5) is to use the non-robust counterpart to slise, i.e., lasso-regression, by setting β = 0. Since lasso is a convex problem, the initial α does not matter, and here we use ordinary least squares regression to obtain the initial α.
With β > 0 the problem is no longer convex and the choice of initial α becomes important. The next two initialisation schemes (lines 6-9 and lines 10-13) offer two different alternatives. In the first scheme the initial α (line 7) is given by non-sparse ols, while in the second scheme α (line 11) is a super-sparse constant vector of only zeros. In both cases the approximation ratio (Algorithm 2) is used to select an initial value for β.
The final scheme (lines 14-27) is inspired by the initialisation used in fastlts (Rousseeuw and Van Driessen 2000) and ransac (Fischler and Bolles 1981), which are related robust regression methods. The idea is to generate u init initial candidates, and heuristically select the best one. The candidates are generated by drawing random subsets (X S , Y S ) of size d + 1, i.e., X S ∈ R d×d+1 and Y S ∈ R d , and fitting linear models to them (using ols, lines 20-22). Note that the probability that at least one of the subsets is free from outliers is given by p where o is the fraction of outliers, d is the number of dimensions, and u is the number of candidates. If d is large, then u would also have to be large to compensate. To alleviate this potential issue we use pca to temporarily reduce the number of dimensions when d is larger than a threshold t pca (lines 17-19). Finally, the best candidate α is the one that minimises the β-Loss and β is updated to match the currently best α (line 23-26).

Asymptotic complexity
The evaluation of the loss function (and its gradient) has a time complexity of O(nd), due to the multiplication between the linear model α and the data-matrix X . The approximation ratio calculation also has a complexity of O(nd) for the same reason. This means that the complexity of the initialisation is dominated by the complexity of ols, which is O (min(nd 2 , d 3 )).
The optimisation consists of owl-qn and graduated optimisation. owl-qn is a variant of lbfgs and has a time complexity of O(md) (Schmidt et al. 2009 (2) error tolerance ε, (3) target approximation ratio K , (4) number of candidates u init , (5) treshold for using pca t pca Result: α, β Algorithm 3: Schemes for initialising α and β.
is the size of the "memory" for approximating the inverse Hessian. Additionally, owlqn requires the loss value and gradient to be calculated every iteration. In practice, the number of iterations can be limited by a constant upper bound (see Sect. 6.4).
The graduated optimisation only adds the approximation ratio calculation, the cost of which is negligible. The total asymptotic time complexity of slise is a combination of the initialisation and the optimisation complexities: O(min(nd 2 , d 3 )+ndp), where p is an upper bound for the total number of optimisation iterations. However, in most cases p d and the ols term becomes vanishingly small (see Sect. 7.1) and in cases where d > p the exact ols solution can be replaced by an optimisation, e.g., using lbfgs (O(ndi max )), since the initialisation does not have to be exact. In this case the complexity of slise becomes O(ndp), or rather O(nd) since p is constant. Also note that any given dataset might require fewer iterations than p, but that is linked to the difficulty of finding the largest linear subset, rather than the size of the dataset.
The memory complexity of slise is at least O(nd), i.e., the memory required to store the dataset itself. The loss function adds O (max(n, d)), the approximation ratio O(n), and owl-qn O(md). This makes the total asymptotic memory complexity of slise O(nd).

Experimental setup
We divide the empirical evaluation into three sections. In this section we describe the datasets and environment used for the experiments. In Sect. 6, we empirically determine suitable default values for the parameters of slise, and which initialisation method to recommend. Later, in Sect. 7 we compare slise to competing methods by demonstrating that (i) slise scales better on high-dimensional datasets than competing methods, (ii) slise is very robust to noise, and (iii) the solutions found using slise are optimal.
The experiments were run using R (Microsoft and R Core Team 2019, v. 3.5) on a high-performance computer cluster (FGCI 2021), using 4 cores from an Intel Xeon E5-2680 at 2.4 GHz, and reserving 16 GB of RAM. Our implementation of the slise algorithm, the code to run all the experiments, and the data processing are released as open source .

Datasets
For the empirical evaluation we use both real and synthetic datasets. An overview of all the datasets is shown in Table 1. The real datasets are three regression datasets from the UCI Machine Learning Repository (Dua and Graff 2019), namely student (Cortez and Silva 2008), air quality (De Vito et al. 2008), and superconductivity (Hamidieh 2018).
As additional real datasets we use some classification datasets for which we create regression tasks with the help of complex classifiers. The datasets and classifier pairs are handwritten digits (Cohen et al. 2017, emnist) with a convolutional neural network, movie reviews (Maas et al. 2011, imdb) with a support vector machine, and particle jets (HIP CMS Experiment 2019, physics) with a neural network. The predictions from the classifiers are probabilities, which we turn into real values using the logit function: y i = log(y/ (1 − y)). emnist has ten classes (digits), which yields ten different regression tasks. Whenever one of these tasks are used, we randomly subsample the dataset such that 50% of the items are of the "correct" digit and 50% are from all the other digits. This creates datasets with lots of potential outliers.
Synthetic datasets are generated as follows. The data matrix X ∈ R n×d is created by sampling from ten normal distribution with randomised means and variances (μ ∼ N (0, 2) and σ 2 ∼ U (0, 1)). The response vector Y ∈ R n is created by y i = a k0 + a k x i + e, where e is normal noise with zero mean and unit variance, and a k ∈ R d+1 is one of seven linear models with the coefficients drawn from a uniform distribution  Björklund et al. (2021) between −1 and 1. Each of the seven models a k is used to create 10% of the Y -values, except one that creates 40% of the Y -values. The code for creating the synthetic datasets is available in the repository ) for full reproducibility.

Pre-processing
slise uses lasso regularisation both for introducing sparsity and regularisation. Since the lasso penalty sums the absolute values of the coefficients, it is important that the variables have roughly the same magnitude (Tibshirani 1996). Normalising the variables also makes it easier to compare the relative importance of values when interpreting the results. Thus, the variables in all datasets, except imdb and emnist, have been centred around (subtracted) the median and scaled (divided) by the mad (median absolute deviation) since this is a more robust form of normalisation than using means and standard deviations (Rousseeuw and Hubert 2011). Furthermore, we add an intercept column to the datasets after the normalisation. For all datasets we also normalise the target (Y ) in the same way, to make the ε:s comparable.
In the student dataset we remove the grades for the first and second period, since these are very correlated with the target (the grades for the third period), as noted by Cortez and Silva (2008).
For the emnist dataset the targets (Y ) are robustly scaled as described above, whereas the pixel values of the input images are just linearly scaled from [0, 255] to [−1, 1]. Some of the pixels have the same value in all of the images (i.e., the corners), which is problematic for some of the comparison methods (in Sect. 7), so these are removed and the images are flattened to vectors of length 672.
The text data in the imdb dataset is transformed into real-valued vectors by using a bag-of-words model, after case normalisation, removal of stop words, removal of punctuation, and stemming. The obtained word frequencies are divided by the most frequent word in each review to adjust for different review lengths, yielding real-valued vectors of length 1000.

Parameter experiments
The slise algorithm presented in Sect. 4 has multiple parameters that must be set. In addition, we also presented four different initialisation schemes. In this section we find good default values for the parameters, and determine which initialisation scheme to recommend. The two most important parameters for slise are the error tolerance ε and the sparsity coefficient λ. These depend on both the use case and the dataset. Therefore, they should be manually adjusted whenever slise is used. The default values for the other parameters, and which initialisation scheme to recommend, are selected based on empirical evidence. Specifically, we base the selection on the value of the loss function and the running time.
All experiments are run ten times per dataset (with different random seeds) in order to better capture the variance (this means that emnist is run a hundred times, due to the ten different tasks/digits). Furthermore, any dataset with n > 10,000 is randomly subsampled to n = 10,000. If nothing else is mentioned, the values in Table 2 are used as default values for the parameters in the experiments.

Error tolerance
The role of the error tolerance is to enable detection and ignoring of outliers. Depending on the dataset or task there might be some obvious limits for the error which should be used when available. The question is hence how to select the value for ε without this knowledge?
The value of ε directly affects the size of the subset that fits the resulting model. Ideally, we want a value of ε such that the subset contains a large amount of nonoutliers and no outliers. In order to better understand how to set the value for ε we consider, in Fig. 2, how ε corresponds to subset sizes for the datasets.
To be able to guarantee the maximum possible robustness we could proceed as Rousseeuw (1984) and select the ε:s that correspond to subset sizes of 50%. However, since none of the datasets considered here contain that many outliers another possibility  is to be less strict and choose subset sizes of, e.g., 75% (Alfons et al. 2013). However, decreasing ε makes the task more difficult and we want to test parameters under adverse conditions. Hence, in this paper we use 30% for all datasets. Furthermore, in the previous paper on slise (Björklund et al. 2019) the main use case is explaining outcomes from black box models, where smaller values of ε might be preferable. The corresponding ε values for all datasets can be seen in Table 3. Another way to investigate the choice of ε is to plot the distribution of errors relative to the size of ε. When ε is small only the peak fits within [−ε, ε] and when ε is large parts of the tails are included. It is natural for distributions to have tails, and not all items in the tails are outliers, but this could nonetheless be used as another criterion for selecting ε, which can be seen in Fig. 3.

Regularisation
The regularisation coefficient λ is dependent on the use case and dataset. With lassotype methods it is common to scale the regularisation by the number of items n (or Fig. 3 The distributions of errors relative to ε can give information on the the effect of ε. When ε becomes large some of the tails of the error distributions fit inside [−ε, ε]. This is a sign that ε might be too large. For example, ε > 0.8 with the physics dataset include tails within the interval (selected subset) by dividing the rest of the loss by n). For slise the size of the subset |S| would be a better choice but that is not known in advance. Furthermore, the rest of the loss function is proportional to ε 2 , so selecting λ ∝ nε 2 makes the transition between dataset sizes and parameter values easier. For the purpose of the experiments, we only use a minimal regularisation by setting λ = 10 −4 · nε 2 , since we are not looking for a specific sparsity.

Initialisation
In Sect. 4 we present four different schemes for selecting the initial values for the linear model α and sigmoid steepness β. Figure 4 shows the results from comparing the four initialisation schemes. No particular method stands out, which indicates that the combination of graduated optimisation and owl-qn yields good performance overall. Furthermore, there are no major differences in running time, according to Table 4. However, slise can only be guaranteed to find a local optimum (in contrast to finding the global optimum), so we need to consider possible failure modes.
Both lasso and ols are non-robust so even a single outlier can lead to arbitrarily large deviations (Donoho and Huber 1983). This means that the starting points might be heavily influenced by outliers. Starting from a vector of zeros is good for sparsity, but with a large enough λ it becomes a local optimum. It is, however, easy to detect when the optimisation has failed to escape this kind of local optimum, by checking if α 1 ≈ 0.
Another problem with using a fixed starting point (which is what the lasso, ols, and zeros initialisation schemes do) is that if the starting point is bad, then there is no way to detect and recover from it. This can be seen in Fig. 5, where the fixed starting α:s are all very close to a cluster of outliers, which creates a local optimum that the optimisation cannot escape. The candidates initialisation is designed to detect  and discard these bad local optima early, by selecting a candidate based on the best β-Loss.
The candidates initialisation scheme is non-deterministic, and requires the number of candidates u init as a parameter. A larger number increases the likelihood that a good candidate is found, but also increases the running time. However, the results from Fig. 6 show that the number of candidates only has a small impact, and also that the time differences are negligible. fast-lts (Rousseeuw and Van Driessen 2000) has a similar notion of candidates, and by default they use u init = 500, which seems to be a reasonable choice also in this case.
With u init = 500 fixed, we can find a value for the threshold t pca for using pca , that is independent of any dataset. The formula for the probability of at least one candidate having no outliers is where o is the fraction of outliers, d = t pca , and u = u init = 500. A larger threshold leads to less information lost in the pca while a smaller value increases p. The curves for different t pca can be seen in Fig. 7. If t pca = 10 then p = 0.38 when o = 0.5, which should be sufficient in most scenarios, since including one outlier does not automatically lead to an inescapable and bad local optimum.
Based on the results in Fig. 4, Table 4, and Fig. 5, we choose candidates as the recommended initialisation scheme. It is the only stochastic initialisation scheme, Fig. 5 Example showing when some of the initialisation schemes fail. The data is constructed such that the starting α:s from lasso, ols, and zeros all pass by a cluster of outliers. These outliers create a local optimum which only the candidates initialisation scheme is able to avoid. slise is used with ε = 0.2 and λ = 0.01

Fig. 6
Losses and running times for different number of initial candidates in the initialisation. For a given dataset, the losses and running times are for all practical purposes equal for different number of initial candidates. The results from multiple runs are aggregated using the mean. Lower is better for both time and loss which means that it does not have a fixed (potentially bad) starting point. Alternatively, the failure mode of zeros is easy to detect, and the initialisation is fast, so it is our second choice if a deterministic algorithm is desired.

Iterations
slise incorporates two iterative optimisation methods, owl-qn and graduated optimisation. Increasing the number of iterations leads to better results for both methods, but beyond a point there are clear diminishing returns. The number of iterations in graduated optimisation is determined by the target approximation ratio K (where K > 1). A larger value results in fewer iterations. However, we require the steps to be small enough that consecutive iterations have similar optima. The number of iterations in owl-qn can be controlled by several different convergence criteria, but the simplest one is to simply limit the number of iterations i max . This has the advantage of ensuring an upper bound on the number of iterations in the worst-case scenario. Additionally, since owl-qn is run multiple times on similar problems there will be a lot of wasted resources if it is forced to fully converge each graduated optimisation iteration.
As can be seen in Fig. 8, K and i max complements each other, in that a decrease in K can be compensated for by a decrease in i max while preserving time and loss values. The choice of values for these parameters can have a massive impact on the running time, while the impact on the loss at times is minimal. Based on the results, the combination of K = 1.15 and i max = 300 is a good default trade-off between time and loss. Furthermore, the last optimisation of owl-qn is not an intermediate step and is, therefore, allowed extra time to converge, by multiplying i max with four (i max = 1200 when β = β max ).

Stopping parameter
It is sufficient that the stopping parameter β max is large enough to make the sigmoid function essentially equivalent to a Heaviside function. As shown in Sect. 3.2, in order to make the shape of the sigmoid only depend on β max it has to be defined as β max = c/ε 2 . The results in Fig. 9 show that c = 20 is sufficiently large, and any larger value merely increases the running time.

Robust regression experiments
In this section we compare slise to five robust regression methods: sparse-lts (Alfons et al. 2013), smdm (Koller and Stahel 2017), conquer (Fernandes et al.  (Qin et al. 2017), and ransac (Fischler and Bolles 1981). lasso (Tibshirani 1996) is also included as a non-robust baseline. To make the comparison maximally useful we compare against implementations found in established software libraries. See Table 5 for an overview of all methods.
All algorithms have been used with default settings, with the exception of sparselts, which has been used with a subset size of n/2 for maximal robustness, and  Tibshirani (1996) glmnet (R) ransac, where we use 20,000 trials and the same error threshold as for slise. In the case of slise, the parameter values are the same as above and can be found in Table 2.

Scalability
First, we investigate the scalability of the methods. Many of the methods have the same asymptotic time complexity, O(nd 2 ), but almost all are iterative methods and the number of iterations are not accounted for in the complexity. We empirically determine the running time on the synthetic datasets with (i) a fixed number of dimensions (d = 100) with an increasing number of items, and (ii) a fixed number of items (n = 10,000) with an increasing number of variables. The results are aggregated from ten different runs. We also limit the calculations to 1000 s. The results are shown in Fig. 10. We observe that slise performs comparable to the other robust regression methods when the number of items increases (the left plot of Fig. 10). However, when the number of variables increases (the right plot of Fig. 10) slise is faster than most robust regression methods. The only exception is conquer which is almost as fast as lasso. The scalability experiment only tests running times on one, synthetic , dataset. To get a broader perspective we also evaluate the running times on the real datasets. The results in Fig. 11 show that slise is comparable or only a couple of seconds slower than most methods on the small datasets. However, on the larger datasets, superconductivity, imdb, and emnist, slise is much faster than the other robust methods, except for conquer. On the larger datasets slise is actually faster than Fig. 10 would suggest, which demonstrates how the difficulty of any given dataset or task affects the running time.

Robustness
Next, we empirically compare the methods' robustness to outliers. To do this we corrupt datasets by replacing a fraction of the items with outliers. We utilise two types of outliers commonly found in literature (Rousseeuw and van Zomeren 1990;Alfons et al. 2013): vertical outliers and leverage points. The dataset we use is a variant of the synthetic dataset, where the Y values are from only one model (so that we can be sure that there are no inherent outliers).
All methods are trained on datasets corrupted by outliers and evaluated on the uncorrupted datasets. If a method is robust to a certain fraction of noise then the residuals for the uncorrupted data will be small. The results are shown in Fig. 12. The breakdown value is the point where the curves start trending upwards, and at high outlier fractions all methods are expected to break down. Vertical outliers are outliers where the target value is different from the rest. We create a vertical outlier by taking a non-outlier item i and replace y i with y i ∼ N (10, 1). As we can see in the left plot of Fig. 12, slise and sparse-lts are the most robust ones. However, vertical outliers are generally considered to be an easier type of outliers (Alfons et al. 2013).
Leverage points are outliers where the variable values are unusual. Here we create a leverage point by taking a non-outlier item i and replace x i with x i ∼ N (10, 1). When the fraction of leverage points is high, most of the correlation between the predictors and target is broken, so regression methods tend to converge towards constant predictions, rather than breaking down. Nonetheless, in the right plot of Fig. 12 we can see that lasso and mte-lasso break down immediately, while conquer and ransac follow soon thereafter. The absolute errors on the clean data from smdm, slise, and sparse-lts stay low for larger fraction of outliers indicating that they are more robust towards leverage points.
Since the robustness experiment is performed on a rather strict dataset we also consider the robustness to outliers on the real datasets. In Fig. 13 we can see how adding vertical outliers affects the behaviour of the regression methods on real datasets. ransac fares well on the low-dimensional datasets, physics and air quality, but fails on even moderately sized datasets. This is because the chance of randomly finding a set of non-outliers shrinks exponentially with the number of dimensions (Fischler and Bolles 1981), so even the high number of trials (20,000) is not enough. On the contrary, slise consistently achieves a breakdown value of at least 0.5.

No outliers
Robust regression methods should also work in situations where there are no outliers. To evaluate this we perform 10-fold cross validation on the real datasets (with no added Fig. 13 Robustness to outliers on real datasets. The x-axis shows the fraction of outliers and the y-axis the mean absolute error on the clean dataset. Some methods were not evaluated for the imdb and emnist datasets, due to their time requirements. Lower error is better outliers). As a baseline we include a dummy model that always predicts the mean yvalue from the training data. In Fig. 14 we see that most robust regression methods (including slise) perform about as well as the non-robust lasso, clearly better than the mean model, the exception being ransac on high-dimensional datasets.

Optimality
Finally, we demonstrate that the solution found by slise optimises the loss of Eq. (1). The slise algorithm is designed to find the largest subset such that the residuals are upper-bounded by ε. To investigate if the model found using slise is optimal, we determine regression models (i.e., obtain the coefficient vectors α) using each algorithm. We then calculate the value of the loss-function in Eq. (1) for every model with varying values of ε.
The results are shown in Fig. 15. All loss-values have been normalised with respect to the absolute median for each value of ε and dataset. slise consistently reaches the smallest loss in the region around ε used for training, as expected. For superconductivity and emnist the loss curves for sparse-lts is very close to the curves for slise, but slise and sparse-lts should actually give equally good, or even identical, results if the error tolerance ε in slise happens to match the subset size h in sparse-lts.

Fig. 14
Cross validation (10-fold) on the real datasets with no added outliers. sparse-lts, slise, and ransac use subset sizes/error tolerances corresponding to 50% of the data. All values have been divided by the corresponding mean absolute error for the reference (mean) model. Some methods were not evaluated for the imdb and emnist datasets, due to their time requirements. Lower error is better Fig. 15 Finding the best solution to Problem 1. The loss-values are normalised by dividing by the absolute median loss per dataset and ε value. The ε used for training slise and ransac is marked with a vertical line. Some methods were not evaluated for the imdb and emnist datasets, due to their time requirements. Lower values are better

Conclusions
This paper refines the slise algorithm for robust regression. slise introduces a novel way of detecting and discarding outliers; find a subset of non-outliers where the error is less than an adjustable error tolerance (ε), for fitting the regression model. This flexible subset size (based on ε) is in contrast to other methods (primarily the lts family) where the subset size is fixed. Additionally, slise yields sparse solutions through built-in regularisation.
Although finding an exact solution to the problem definition (Problem 1) is NPhard (see Sect. 2.1), the combination of graduated optimisation with a quasi-Newton optimiser (owlqn ) yields an effective approximation (see Sects. 3 and 7). Adding a stochastic initialisation further mitigates the risk of unfortunate starting conditions (see Sect. 6.3).
When comparing to other robust regression methods, slise is able to achieve robustness levels that are among the best possible, which we show both theoretically (see Sect. 2.2) and empirically (see Sect. 7.2). Furthermore, slise scales better to highdimensional data than many alternative methods (see Sects. 4.2 and 7.1).
Future work could investigate how the choice of λ affects the sparsity, especially during the graduated optimisation. Another direction would be to try changing the balance between maximising the subset and minimising the residuals, or to introduce different weights for the data items.
In an earlier paper (Björklund et al. 2019) we show how slise can be used to explain outcomes from black box models in a way that respects constraints in the data. Along this line we could further investigate the utility of selecting the data used for the explanations in order to answer specific questions in an interactive manner. This would give better insight into the learned models and their behaviour.
The explanations given by slise are local, i.e. for specific outcomes, and an interesting follow up would be to combine these local explanations into one global explanation. Furthermore, slise is a robust regression method and, therefore, quite generic, which means that it can readily be integrated into, or combined with, other explanation methods.
Our implementation of the slise algorithm is released under an open source license. It is available in both Python (Björklund 2021) and R , which also includes the code for running all the experiments.