Conceptio › Archive › arXiv CS

arXiv CSopen access

Computation of Least Trimmed Squares: A Branch-and-Bound framework with Hyperplane Arrangement Enhancements

Unknown · 2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

Open Source ↗Direct PDF ↓

machine learning, deep learning, neural networks

OPERATIONS RESEARCH https://pubsonline.informs.org/journal/opre Vol. 00, No. 0, Xxxxx 0000, pp. 000–000

Submitted to Operations Research

ISSN 0030-364X, EISSN 1526-5463

arXiv:2604.11584v1 [math.OC] 13 Apr 2026

Computation of Least Trimmed Squares: A Branch-and-Bound framework with Hyperplane Arrangement Enhancements Xiang Meng MIT Operations Research Center, [email protected]

Andrés Gómez USC Daniel J. Epstein Department of Industrial and Systems Engineering, [email protected]

Rahul Mazumder MIT Sloan School of Management, Operations Research Center, [email protected]

Authors are encouraged to submit new

Abstract. We study computational aspects of a key problem in robust statistics —

papers to INFORMS journals by means

the penalized least trimmed squares (LTS) regression problem, a robust estimator

of a style file template, which includes

that mitigates the influence of outliers in data by capping residuals with large magni-

the journal title. However, use of a tem-

tudes. Although statistically attractive, penalized LTS is NP-hard, and existing mixed-

plate does not certify that the paper

integer optimization (MIO) formulations scale poorly due to weak relaxations and

has been accepted for publication in

exponential worst-case complexity in the number of observations. We propose a new

the named journal. INFORMS journal

MIO formulation that embeds hyperplane arrangement logic into a perspective refor-

templates are for the exclusive purpose

mulation, explicitly enforcing structural properties of optimal solutions. We show

of submitting to an INFORMS journal

that, if the number of features is fixed, the resulting branch-and-bound tree is of poly-

and are not intended to be a true repre-

nomial size in the sample size. Moreover, we develop a tailored branch-and-bound

sentation of the article’s final published

algorithm that uses first-order methods with dual bounds to solve node relaxations

form. Use of this template to distribute

efficiently. Computational experiments on synthetic and real datasets demonstrate

papers in print or online or to submit

substantial improvements over existing MIO approaches: on synthetic instances with

papers to another non-INFORM publi-

5000 samples and 20 features, our tailored solver reaches a 1% gap in 1 minute while

cation is prohibited.

competing approaches fail to do so within one hour. These gains enable exact robust regression at significantly larger sample sizes in low-dimensional settings. Funding: This research was supported by AFOSR Grant No. FA9550-24-1-0086 and NSF Grant No. 2346058. RM acknowledges research funding from ONR N00014-25-1-2504. Key words: Robust statistics, mixed-integer nonlinear optimization, branch-and-bound, hyperplane arrangements

Meng, Gómez and Mazumder: Computation of LTS: BnB and Hyperplane Arrangements Article submitted to Operations Research

1. Introduction Outliers—informally speaking, observations that deviate substantially from the bulk of the data— are pervasive in modern datasets. They arise from sensor failures, human entry mistakes, transient process instabilities, or simply corruptions in data. For example, in regression problems, contaminated observations or outliers can severely affect the quality of statistical estimates compared to estimates obtained from clean data. Classical least squares estimation is notoriously sensitive to outliers: a single high-leverage observation can arbitrarily distort the regression coefficients, potentially hurting model prediction and/or statistical inference (Huber 1981, Maronna et al. 2006). Ridge regression (i.e., least squares penalized with a squared ℓ2 -norm penalty) — despite its regularization benefits for ill-conditioned problems, inherits this sensitivity since the squared loss grows without bounds with the residual magnitude. To address this limitation, several robust estimators have been studied in the statistics literature (Rousseeuw and Leroy 1987). Given a model matrix X ∈ Rn×p and response vector y ∈ Rn , consider the optimization problem n n X 2 λ 1X 2 ⊤ yi − xi β (1 − zi ) + ∥β∥2 + µ zi , min β∈Rp ,z∈{0,1}n 2 2 i=1 i=1

(1)

p ⊤ where x⊤ i denotes the i-th row of X, β ∈ R denotes the regression coefficients and ri := yi −xi β,

i ∈ [n] are the residuals. We are interested in the traditional setting of robust statistics where the number of observations is larger than the dimension p, that is n > p. Binary variable zi ∈ {0, 1} indicates whether observation i is an outlier—the weight (1 − zi ) serves to exclude large absolute values of the residual ri (for i ∈ [n]) and include all other residuals in the computation of the regression coefficients. The excluded points correspond to outliers, and the remaining points are inliers. The parameter λ ≥ 0 controls ridge-regularization, while µ ≥ 0 is a parameter that penalizes the number of outliers. Problem (1) is a penalized variant of the Least Trimmed Squares (LTS) estimator introduced by Rousseeuw (1984), which fits a model by minimizing the sum of the smallest squared residuals. The cardinality-constrained LTS places a constraint on the number of trimmed observations, while (1) represents its penalized form that allows the trimming level to be implicitly determined by the penalty µ. Compared to other robust estimators such as Least Median of Squares (LMS)(Rousseeuw 1984), LTS has desirable statistical properties — it enjoys √ n-consistency, and higher asymptotic efficiency (Rousseeuw and Leroy 1987, Rousseeuw and Van Driessen 2006a).

Meng, Gómez and Mazumder: Computation of LTS: BnB and Hyperplane Arrangements Article submitted to Operations Research

The estimator in problem (1) can also be interpreted through the lens of a robust loss functions. Indeed, upon minimizing Problem (1) over the binary indicators z, the problem reduces to: min

β∈Rp

n X

ϕcap yi − x⊤ i β

i=1

n1 o λ def 2 2 + ∥β∥2 , where ϕcap (r) = min r , µ . 2 2

The objective has a capped quadratic loss ϕcap (r)

Comparison of robust losses for a single residual

with an additional ridge regularization with penalty

stant level µ for large (absolute) residuals, discouraging grossly contaminated observations from dominating the fit. The nonconvex capped quadratic belongs to a broader class of robust losses that temper growth

Huber Tukey biweight Geman-McClure

2.0

parameter λ/2 ≥ 0. The loss r 7→ ϕcap (r) is quadratic

Cauchy Capped quadratic

1.5 Loss

for small (absolute) residuals but saturates at the con-

(2)

1.0 0.5 0.0 3

0 Residual r

Figure 1: Various robust loss functions.

in the tails, as illustrated in Figure 1, and has close ties to the redescending M-estimators (Maronna et al. 2006). A canonical example of a convex loss function is the Huber loss, which transitions from quadratic to linear growth and is a standard baseline in robust regression (Huber 1964, 1981). Well known bounded, nonconvex choices of the loss function include Tukey’s biweight (Holland and Welsch 1977), Geman–McClure (Chang et al. 2012), and Cauchy/Lorentzian losses (Motulsky and Brown 2006). Relative to these alternatives, the objective in (2) offers a transparent interpretation: observations with sufficiently large residuals contribute a constant effect µ and are effectively trimmed from the fit. Despite its attractive properties, solving problem (1) poses a formidable challenge. The problem is NP-hard (Bernholt 2006), and early approaches relied on heuristics such as such as the feasible solution algorithm (Hawkins 1994) and its subsequent refinement (Hawkins and Olive 1999), or the widely used FAST-LTS algorithm (Rousseeuw and Van Driessen 2006a). While efficient, these methods provide no optimality guarantees and can produce solutions that are far away from the optimum (Gómez and Neto 2025). In addition, the desirable statistical properties of the LTS estimator are predicated on obtaining a global minimizer of (1), and a suboptimal solution need not inherit any of these guarantees. Problem (1) can be formulated as a mixed-integer quadratic program, enabling globally optimal solutions via modern mixed-integer optimization (MIO) solvers (Gómez and Neto 2025, Sun et al. 2021, Zioutas and Avramidis 2005a). However, these methods have a worst case complexity of O(2n ). While MIO methods often exhibit practical runtimes that are substantially better than

Meng, Gómez and Mazumder: Computation of LTS: BnB and Hyperplane Arrangements Article submitted to Operations Research

the theoretical worst-case complexities; unfortunately, for the case of (1), convex relaxations are typically weak and branch-and-bound algorithms require a prohibitive number of nodes to demonstrate optimality. It appears that existing MIO methods for LTS can deliver optimal solutions to problem instances with n ≲ 100 and p ≲ 10 (practical performance is also highly dependent on hyperparameters λ and µ and the dataset itself). Outside of the MIO literature, some exact combinatorial algorithms have been proposed. Agulló (2001a) proposed both a probabilistic exchange algorithm and an branch-and-bound algorithm for solving LTS. Hofmann et al. (2010) designed a row adding algorithm that computes exact LTS solutions for a range of coverage values, and Klouda (2015) studied an exact algorithm for solving LTS whose computational cost scales as O(np+1 ). These exact combinatorial methods are generally limited to relatively small or low-dimensional problem instances. In a different line of work, exact methods based on topological sweeps for solving LTS have been proposed in the statistics literature. These methods inspired by computational geometry are based on an arrangement of hyperplanes, and require solving O(np ) least squares problems. While topological sweep methods can in principle scale to a large number of points provided that p is small, we are unaware of practical implementations apart from the special case of p = 2 (e.g., Edelsbrunner and Souvaine 1990, Hössjer 1995). However, as far as we can tell, there are no available practical implementations for the topological sweep method for p ≥ 3. Contributions In this paper we propose new MIO formulations and algorithms for Problem (1), designed with the common regime of n ≫ p in mind. The proposed formulations build on state-of-the-art perspective reformulations for Problem (1) with strong convex relaxations (Gómez and Neto 2025), but additionally include additional constraints enforcing hyperplane arrangements logic. The proposed formulation avoids big-M constraints altogether, and branch-and-bound methods based on the proposed formulation terminate after exploring at most O(np+1 ) nodes. Unlike topological sweep methods, the practical runtime of branch-and-bound algorithms appear to be substantially better than the worst-case performance. In fact, we propose the first practical implementation of methods based on hyperplane arrangement for Problem (1) that can handle instances with p ≥ 3. The proposed formulation leads to at least an order-of-magnitude improvement over alternative state-of-the-art MIO approaches (Gómez and Neto 2025, Insolia et al. 2022), and the approach of Bertsimas and Mazumder (2014) for the LMS problem.

Meng, Gómez and Mazumder: Computation of LTS: BnB and Hyperplane Arrangements Article submitted to Operations Research

In addition, we also develop a tailored branch-and-bound solver to tackle the proposed formulation. Our algorithm is developed in Python, uses first-order methods to solve the continuous subproblems and requires no commercial software. Our code builds upon and significantly extends the earlier work of Hazimeh et al. (2021) proposing tailored branch-and-bound approaches for a different problem sparse linear models (See Section 4 for differences with this work). The computational performance of our algorithm is competitive with state-of-the-art commercial solver Gurobi to solve the proposed formulation, and appears substantially faster than using other commercial MIO solvers. On synthetic datasets, our method scales to significantly larger sample sizes: on instances with n = 5000 samples and p = 20 features, it reaches a 1% optimality gap in 1 minute, whereas no competing approach does so within a one-hour time limit. On 13 real benchmark regression datasets, our method and Gurobi applied to the proposed formulation achieve the strongest performance, each attaining the best runtime on roughly half of the instances; the remaining baselines are substantially slower or fail to close the gap within the time limit. Outline The rest of the paper is organized as follows. In Section 2 we review existing methods for (1) (MIO and topological sweep methods as well as heuristics). In Section 3 we describe the proposed MIO formulation, study its strength and prove its worst-case complexity of O(np+1 ). In Section 4 we describe the proposed branch-and-bound algorithm based on first-order methods. Finally, in Section 5 we provide computational experiments with synthetic and real data.

2. Review of existing solution approaches We review exact solution approaches for Problem (1) from the literature. First, we discuss exact MIO approaches that have recently been proposed in the Operations Research community. We then discuss classical exact solution approaches based on hyperplane arrangements from the statistics literature. Finally, we examine approximate methods (also known as heuristics) for finding good solutions, which are currently the most widely used approach in practice. 2.1. MIO formulations At a high level, MIO approaches are based on the branch-and-bound method. The natural approach is to branch on the binary variables z, that is, branching on whether a data point should be deemed an outlier or not. Earlier approaches from the mathematical optimization community (Agulló 2001b, Giloni and Padberg 2002) studied formulation (1) directly. However, nonconvexity of the

Meng, Gómez and Mazumder: Computation of LTS: BnB and Hyperplane Arrangements Article submitted to Operations Research

continuous relaxation of Problem (1) due to the presence of cubic terms in the objective makes the design of efficient solution methods challenging. The first exact MIO approaches were proposed by Zioutas and Avramidis (2005b) and Zioutas et al. (2009), which reformulated (1) into mixed-integer quadratically constrained quadratic programs by introducing additional variables representing the loss incurred from fitting each point. More recently, Insolia et al. (2022) proposed an alternative mixed-integer quadratic optimization reformulation by introducing variables representing “corrections” associated with each datapoint, which we review in §2.1.1. Finally, stronger formulations were proposed in the literature, see Gómez (2021) for methods specific to time-series data and Gómez and Neto (2025) for approaches for the general least trimmed squares problem (1), which we review in Section 2.1.2. In general, MIO methods have a theoretical worst case complexity of O(2n ), corresponding to all possible combinations of inliers and outliers. Branch-and-bound methods typically explore substantially fewer branch-and-bound nodes (especially if strong relaxations are used). Nonetheless, both the theoretical complexity and computational experiments reported in the aforementioned papers suggest that the performance of existing MIO methods deteriorates substantially as the sample size n increases. Current MIO approaches have been shown to perform well with real data for n ≈ 100, but struggle to solve to optimality instances with larger values of n. 2.1.1. Big-M Formulation We first consider a MIO formulation of Problem (1) similar to the one proposed by Insolia et al. (2022), incorporating Big-M constraints. We introduce auxiliary variables w ∈ Rn representing corrections to the response vector, where wi captures the residual when observation i is treated as an outlier—that is, wi = yi − x⊤ i β. We define a constant M > 0 such that any optimal correction satisfies |wi∗ | ≤ M . Binary variables zi , i ∈ [n], indicate whether observation i is an outlier, enforced through the constraints −M zi ≤ wi ≤ M zi , i ∈ [n]. This results in the following Big-M formulation: n n X 2 λ 1X ⊤ 2 yi − wi − xi β + ∥β∥2 + µ zi min β,w,z 2 2 i=1 i=1

s.t. − M zi ≤ wi ≤ M zi ,

(3)

i ∈ [n],

β ∈ Rp , w ∈ Rn , zi ∈ {0, 1},

i ∈ [n].

In the above formulation, if zi = 0, then the big-M constraints force wi = 0; in that case, the loss associated with datapoint i reduces to the standard quadratic loss. On the other hand, if zi = 1, then wi is free to take any value; in particular, for a fixed β, setting wi = yi − x⊤ i β is optimal, and the associated loss with datapoint i vanishes.

Meng, Gómez and Mazumder: Computation of LTS: BnB and Hyperplane Arrangements Article submitted to Operations Research

There are two concerns in using formulation (3). Firstly, we need to specify the value of the maximum magnitude M , which can be large for arbitrary outliers. This is different from formulation (1), which doesn’t need such a specification. If we select M to be sufficiently large so that a solution to Problem (3) is also a solution to Problem (1) (e.g., using arguments similar to Bertsimas et al. 2016), the value of M can be very conservative, possibly causing numerical issues for an algorithm. Additionally, convex relaxations of (3) are weak since solution β = 0, w = y, z = ϵ1 are feasible provided ϵ ≥ ∥y∥∞ /M , with objective values close to 0 for any reasonable value of M . As a consequence, branch-and-bound algorithms struggle to prune nodes and make informed branching decisions, leading to prohibitive solution times. 2.1.2. Perspective Formulation Gómez and Neto (2025) propose stronger, big-M free formulations of Problem (1), obtained by applying the perspective reformulation (Frangioni and Gentile 2006, Günlük and Linderoth 2010). The key observation is that the objective of Problem (3) can be written as 1

λ 1 1 ∥y − w − Xβ∥22 + ∥β∥22 = ∥y∥22 − y ⊤ (Xβ − w) + β⊤ w⊤ 2 2 2 2   ⊤ ⊤ X X + λI −X . Given d ∈ Rn+ , define where Σ =  −X I   X ⊤ X + λI −X ⊤ ed =   Σ −X I − 2Diag(d)

  β Σ , w

and note that the following quadratic term can be decomposed as     n X β β 1 1 ⊤ ⊤ Σ  = ⊤ ⊤ Σ e   + di wi2 , β w β w d 2 2 w w i=1

(4)

e d is positive semi-definite. The perwhere the first term remains convex provided that matrix Σ P P spective reformulation then replaces the separable term i di wi2 with i di wi2 /zi . The latter term equals di wi2 when zi = 1 and enforces wi = 0 when zi = 0, using the convention that 0/0 = 0 and x/0 = ∞ for x > 0. We thus obtain the perspective formulation   n 2 X β 1 1 d w i i 2 ⊤ ⊤ ⊤ ed  + µzi + min ∥y∥2 − y (Xβ − w) + β w Σ β,w,z 2 2 zi w i=1 (5) s.t. − M zi ≤ wi ≤ M zi ,

i ∈ [n],

β ∈ Rp , w ∈ Rn , zi ∈ {0, 1},

i ∈ [n].

Meng, Gómez and Mazumder: Computation of LTS: BnB and Hyperplane Arrangements Article submitted to Operations Research

Observe that if d > 0, then the big-M constraints −M z ≤ w ≤ M z can be removed from the formulation. The perspective formulation (5) is equivalent to (3) but yields tighter continuous relaxations, leading to smaller branch-and-bound trees and faster solve times. To find a vector d that results in the best convex relaxation of Problem (5), Gómez and Neto (2025) propose an algorithm that solves a sequence of positive semi-definite programs (SDPs) with cones of order p + 1. We now describe an alternative and simpler approach we use in our numerical experiments. First, let λ̃ ∈ (0, λ) be a small number used to guarantee strong convexity ed ≻ 0 and numerical stability. We note that, by the Schur complement, we can guarantee that Σ holds by ensuring that I − X(X ⊤ X + (λ − λ̃)I)−1 X ⊤ − 2Diag(d) ⪰ 0. Thus, setting 1 dj = γmin I − X(X ⊤ X + (λ − λ̃)I)−1 X ⊤ 2

∀j ∈ [n],

(6)

e d ≻ 0. where γmin (·) denotes the minimum eigenvalue of the argument matrix, guarantees that Σ While choice (6) results in a weaker relaxation than using the method proposed by Gómez and Neto (2025), this can be computed faster and does not require access to (commercial) conic optimization solvers or specialized algorithms. 2.2. Solution via hyperplane arrangements Different from MIO approaches, exact algorithms for (1) were proposed in the statistics literature (Rousseeuw and Leroy 2003) based on hyperplane arrangements. The algorithms involved solving O(np ) least squares regression problems. Thus, unlike MIO methods, the complexity is polynomial in the number of datapoints n but exponential in the number of features p. As far as we can tell, these algorithms are difficult to implement in practice. We now review methods based on hyperplane arrangements. We note that the method we describe differs from classical algorithms for the least trimmed squares problems, since we describe a method for the penalized version whereas classical approaches tackle cardinality constrained problems with an upper bound on the number of outliers. Nonetheless, the ideas behind the algorithms are similar. The hyperplane arrangement algorithm relies on two key results: the first (Proposition 1) concerns structure of optimal solutions of (1), and the second (Proposition 2) is a classical result from combinatorial geometry (Cover 1965, Zaslavsky 1975). P ROPOSITION 1. Any optimal solution (β ∗ , z ∗ ) of Problem (1) satisfies: √ ∗ 1. if zi∗ = 0, then |yi − x⊤ 2µ; i β |≤ √ ∗ 2. if zi∗ = 1, then |yi − x⊤ 2µ. i β |≥

Meng, Gómez and Mazumder: Computation of LTS: BnB and Hyperplane Arrangements Article submitted to Operations Research

P ROPOSITION 2. An arrangement of n hyperplanes in Rp passing through the origin divide the P n p space into at most 2 p−1 i=0 i = O(n ) regions, where the upper bound is attained if the hyperplanes are in general position. We can use Proposition 2 to solve Problem (1) as follows. Each datapoint (xi , yi ) induces three regions in Rp : n p o Bi≤ = β ∈ Rp : x⊤ β ≤ y − 2µ i i n p o p β ≤ y + 2µ Bi= = β ∈ Rp : yi − 2µ ≤ x⊤ i i n p o β ≥ y + Bi≥ = β ∈ Rp : x⊤ 2µ . i i

Note that for any solution β in the band defined by Bi= , setting zi = 0 is a better choice than zi = 1 in Problem (1)—that is, point i is an inlier. Similarly, for any β ∈ Bi≤ ∪ Bi≥ we have that zi = 1 is a better choice in Problem (1), thus point i would be categorized as an outlier. The three regions √ √ are defined by two hyperplanes, x⊤ 2µ and x⊤ 2µ corresponding to data i β = yi − i β = yi + point i ∈ [n]. Thus, combining all 2n hyperplanes creates a division of Rp into O((2n)p ) regions, each corresponding to a selection of outliers. Hence an optimal solution to Problem (1) can be obtained by solving a ridge regularized least squares problem for each one of the regions (with a fixed selection of outliers) and choosing the solution that yields the best objective value. E XAMPLE 1. Consider three points in R2 : (x1 , y1 ) = (−1, 8), (x2 , y2 ) = (0, 0.7) and (x3 , y3 ) = (1, 1). Figure 2 (left) shows the three points in the plane as well as the optimal ridge regularized LTS regression line y = β0 + β1 x, obtained by solving (1) with λ = µ = 1 – note that the first point is flagged as an outlier in the optimal solution (β0∗ , β1∗ ) = (0.8, 0.1). Figure 2 (right) shows the associated hyperplane arrangement in the β space: each point i is associated with a band Bi= , and for any solution (β0 , β1 ) ∈ Bi= the best choice is to keep point i as an inlier; point i is discarded as an outlier for choices (β0 , β1 ) ̸∈ Bi= . Since the three bands do not intersect, we can conclude that z = 0 cannot be optimal, thus the solution with all three points as inliers should not be considered. ■ Note that for the most part, the algorithmic idea of traversing all regions while solving least squares problems is is computationally expensive and hence may be impractical. Efficient implementations do exist if p = 2 (e.g., Edelsbrunner and Souvaine 1990, Hössjer 1995): using the

Meng, Gómez and Mazumder: Computation of LTS: BnB and Hyperplane Arrangements Article submitted to Operations Research

10 2

0 -1.5

-1

-0.5

0.5

1.5

-4

β1

-2

-6

-15

-10

-5

-8 -10

Figure 2

point 1 point 2 point 3

-10

β0

Illustration of Example 1: an ridge regularized LTS regression in the plane (left), and associated hyperplane arrangement (right).

observation that the intersection of two lines occurs at a single point, one can efficiently enumerate all regions with proper data structures. Unfortunately, we are not aware of practical implementations for p > 2. In theory, it is possible to do so via linear optimization for example (Černỳ et al. 2019), but the additional computational cost can be excessive. Moreover, the computational cost of computing O(np ) least squares estimators can be prohibitive even if p = 3 for sufficiently large values of n. 2.3. Heuristics Given the perceived inefficiencies of exact methods, the preferred solution approach for Problem (1) is via heuristics. The most popular heuristic, called FAST-LTS (Rousseeuw and Van Driessen 2006b), is based on alternating minimization and works as follows. Starting from an initial point β, the heuristic alternates between: (i) flag the points with largest residuals as outliers (optimization over z); (ii) solve a least squares problem for the new inlier/outlier combination (optimization over β). Variants of this method are common in the literature for problems similar to Problem (1) (Shen and Sanghavi 2019a,b). While alternating minimization heuristics have been observed to work well in high signal-to-noise ratios, and performance guarantees can be established under appropriate conditions (Bhatia et al. 2015), they have also been shown to lead to suboptimal solutions in challenging instances (Gómez and Neto 2025).

3. Enhanced MIO formulation with hyperplane arrangements The perspective formulation (5) treats each observation independently: it strengthens the relaxation of each wi2 term, but does not exploit the relationship between the outlier indicators z and the regression coefficients β. Hyperplane arrangements, on the other hand, capture precisely this

Meng, Gómez and Mazumder: Computation of LTS: BnB and Hyperplane Arrangements Article submitted to Operations Research

relationship through the threshold structure of Proposition 1. We propose a formulation that incorporates this structure into the perspective formulation, by enforcing the logical relationships

zi = 0 =⇒ β ∈ Bi=

(7a)

zi = 1 =⇒ β ∈ Bi≤ ∪ Bi≥

(7b)

for all datapoints i ∈ [n]. A direct approach to enforce (7) is by including additional binary variables z − , z + ∈ {0, 1}n such that z − + z + = z and in the big-M constraints

p p 2µ − M (1 − zi+ ) ≤ x⊤ 2µ + M (1 − zi− ) i β ≤ yi − p p yi − 2µ − M (zi− + zi+ ) ≤ x⊤ 2µ + M (zi− + zi+ ) i β ≤ yi + yi +

(8a) (8b)

for all i ∈ [n]. The additional variables z − , z + have the interpretation that “zi+ = 1 if and only if β ∈ Bi+ ” and “zi− = 1 if and only if β ∈ Bi− ”. Note that constraints (8) are not valid for (3) or (5) in the usual sense, as they remove feasible integer points. Nonetheless, they do not remove any optimal solutions of the MIO problems and can be used as optimality cuts or constraints. A limitation of constraints (8) is the presence of the big-M terms. Indeed, the constant M is required to be at least as large as the maximum residual, which in principle can be arbitrarily large. We now propose an alternative formulation that imposes the logical considerations (7) on the perspective formulation of the LTS problem (5) without using big-M constraints:

1 ∥y∥22 − y ⊤ (Xβ − w) + β⊤ w⊤ β,w,z 2 2 min

w± ,z ±

  n X β di (wi+ )2 di (wi− )2 e   Σd + µzi + + zi+ zi− w i=1

p wi+ ≥ 2µ zi+ , i ∈ [n], p |yi − x⊤ 2µ(1 − zi− − zi+ ), i ∈ [n], i β − wi | ≤

(9b)

zi = zi+ + zi− ,

(9c)

s.t. wi− ≥

p 2µ zi− ,

wi = wi+ − wi− ,

i ∈ [n],

β ∈ Rp , w ∈ Rn , z ∈ {0, 1}n , w− , w+ ∈ Rn+ , z − , z + ∈ {0, 1}n .

(9a)

(9d)

Meng, Gómez and Mazumder: Computation of LTS: BnB and Hyperplane Arrangements Article submitted to Operations Research

Note that in Problem (9), the optimization variables z and w can be projected out and constraints (9c) removed, but we keep them for now to facilitate comparisons with Problem (5). All constraints in problem (9) are either linear or can be linearized in the case of (9b). Intuitively, variables w+ and w− represent the positive and negative parts of w, respectively; variable zi+ = 1 (zi− = 1) if point i is an outlier by overestimating (underestimating) the response variable. The rest of this section is devoted to studying formulation (9). 3.1. Formulation Strength The proposed formulation (9) is at least as strong as the perspective reformulation (5), and avoids using big-M constraints such as (8). The perspective reformulation (5) is obtained from the convex hull of the epigraph of the quadratic regularization term and indicator constraints, that is, Zpersp = (z, w, t) ∈ {0, 1} × R2 : t ≥ w2 , w(1 − z) = 0 . Our formulation is based on a richer structure, involving a residual ri corresponding to a term of the form yi − x⊤ i β for some i ∈ [n]. Namely, given b ∈ R+ , we consider the set ZHA (b) = (z, w, r, t) ∈ {0, 1} × R3 : t ≥ w2 , w(1 − z) = 0, |r|(1 − z) ≤ b(1 − z), |r|z ≥ bz, w = rz . Using decomposition (4), we can reformulate Problem (1) as the MIO

  n X β e d  + (µzi + di ti ) Σ w i=1

1 ∥y∥22 − y ⊤ (Xβ − w) + β⊤ w⊤ β,w,z,t 2 2 p s.t. (zi , wi , yi − x⊤ ∀i ∈ [n] i β, ti ) ∈ ZHA ( 2µ)

(10b)

β ∈ Rp , w ∈ Rn , z ∈ {0, 1}n , t ∈ Rn .

(10c)

min

(10a)

Note that constraints (10b) include the epigraph constraint ti ≥ wi2 . In addition: √ 2µ, encapsulating 1. if zi = 0, then the constraint in (10b) reduces to wi = 0 and |yi − x⊤ i β| ≤ precisely relationship (7a). 2. if zi = 1, then the constraint (10b) reduces to |yi − x⊤ i β| ≥

√

2µ –corresponding precisely to

(7b)– and w = yi − x⊤ i β, encapsulating the additional optimality condition that the additional variables wi are exactly the residuals to ensure that the associated term vanishes.

Meng, Gómez and Mazumder: Computation of LTS: BnB and Hyperplane Arrangements Article submitted to Operations Research

Naturally, set ZHA (b) is highly nonconvex due to the presence of binary variables and quadratic constraints. Thus, instead of relying directly on it, our formulation relaxes constraints (10b) to √ β, t ) ∈ cl conv Z ( (zi , wi , yi − x⊤ 2µ) , ∀i ∈ [n]. The description of this convex hull is given i HA i in Proposition 3, and the proof is deferred to Appendix A.1. P ROPOSITION 3. The closure of the convex hull of ZHA (b) is given by cl conv (ZHA (b)) = proj(z,w,r,t) Z̄(b) where n (w+ )2 (w− )2 Z̄(b) = (z, w, r, t, w− , w+ , z − , z + ) : t ≥ + + − , w− ≥ bz − , w+ ≥ bz + , z z o |r − w| ≤ b(1 − z − − z + ), z − + z + = z ≤ 1, w = w+ − w− , z − , z + ≥ 0 .

The constraints defining set Z̄(b) in Proposition 3 are precisely (9a)-(9c). We formalize this property in the following corollary. C OROLLARY 1. Formulation (9) is a correct formulation for Problem (1), obtained from Prob√ lem (10) by replacing constraints (10b) with (zi , wi , yi − x⊤ i β, ti ) ∈ cl conv ZHA ( 2µ) , ∀i ∈ [n]. So far we have established that formulation (9) is strong, in the sense of directly using the convex hull of a region involving a mix of feasibility and optimality conditions. The formulation is indeed the strongest in its class, that is, no stronger formulation can exist unless additional structure is included in the set (e.g., considering sets with multiple datapoints with non-trivial interactions). However, we state in Proposition 4 that the proposed relaxation does not improve the continuous relaxation of (5). This result is a consequence of the more general Proposition 7, and we defer its proof until then. P ROPOSITION 4. Constraints (9a)-(9b) are redundant for the continuous relaxation of Problem (9). C OROLLARY 2. The optimal solutions and optimal objective values of the convex relaxations of formulations (5) and (9) coincide. In light of Proposition 4, preferring formulation (9) over formulation (5) seems counterintuitive: the continuous relaxation is harder to solve due to the presence of additional constraints, and there appear to be no associated gains in terms of relaxation quality. However, we note that the result of Proposition 4 holds only for the relaxation of Problem (9), and need not hold if any additional

Meng, Gómez and Mazumder: Computation of LTS: BnB and Hyperplane Arrangements Article submitted to Operations Research

constraint is added or objective terms are modified. In particular, the result of Proposition 4 does not hold when additional branching constraints of the form zi± ≤ 0 or zi± ≥ 1 are introduced. Thus, even if the root relaxations obtained from both formulations coincide, using formulation (9) results in better relaxations at every node other than the root node of the branch-and-bound tree. In fact, as we show in the next section, the inclusion of constraints (9a)-(9b) guarantees that branch-andbound algorithms explore a number of nodes polynomial in the number of datapoints n when the dimension p is fixed. 3.2. Computational Cost Every time a branching constraint zi± ≤ 1 or zi± ≥ 1 is introduced, the feasible values for β ∈ Rp are restricted: the regression coefficients are forced to be on one side of a hyperplane in the hyperplane √ arrangement induced by {x⊤ 2µ}ni=1 . Every time a branching decision is incompatible i β = yi ± with a region of the hyperplane arrangement, then the branch-and-bound node is pruned due to being infeasible for Problem (9). In contrast, with formulation (5), no branch-and-bound nodes are infeasible, and pruning only occurs due to bounding or by finding integer solutions. In Proposition 5 we show that, thanks to the possibility of fathoming by infeasibility, branch-and-bound algorithms run in time polynomial in n. P ROPOSITION 5. If formulation (9) is solved via the branch-and-bound method where: • all convex interval relaxations are solved to optimality at each node of the branch-and-bound tree, • variable dichotomy is used for branching, that is, branching is performed only via disjunctions “zi± ≤ 0 or zi± ≥ 1”, and • only fractional variables are selected for branching, then the branch and bound tree has at most O (min{4n , np+1 }) nodes. Proof The computational cost of O (4n ) is obtained since there are 2n binary variables, and a binary tree of depth 2n has at most 2 · 22n − 1 nodes. To prove the complexity of O (np+1 ), assume the worst-case scenario that the branch-and-bound algorithm never fathoms nodes by bounding, and pruning only occurs at integer or infeasible solutions; these integer or infeasible nodes are thus the only leaves of the trees. Observe that there is a one-to-one correspondence between feasible integer solutions of Porblem (9) and regions of the √ n hyperplane arrangement induced by {x⊤ i β = yi ± 2µ}i=1 . Thus, from Proposition 2 we find that there are at most O(np ) leaves corresponding to integer feasible solutions.

Meng, Gómez and Mazumder: Computation of LTS: BnB and Hyperplane Arrangements Article submitted to Operations Research

We now count the leaves corresponding to infeasible solutions. Consider an arbitrary leaf node that was pruned by infeasibility, which we refer to as the infeasible leaf. Since we assumed the convex relaxations are solved to optimality, pruning by infeasibility is identified immediately, as soon as a subproblem is infeasible. In other words, the parent of the infeasible leaf has a feasible continuous relaxation, that is, there exist β̄ ∈ Rp satisfying all constraints (9a)-(9b) and additional branching constraints added to reach that node. Therefore, this parent has a feasible descendant; in particular, there is a (unique) descendant obtained by branching at each step according to β̄ until an integer solution is found. We say that this unique feasible descendant and the original infeasible node are relatives to each other. Finally, note that an integer node at depth q of the tree can have at most q − 1 infeasible relatives, where the upper bound is attained in the case shown in Figure 3. Since the maximum depth is 2n and there are at most O(np ) integer nodes, it follows that there are at most O(2np+1 ) infeasible nodes, at most O(np + 2np+1 ) leaves and twice as many nodes in a □

binary tree, proving the result.

inf

int

Figure 3

inf

A depth 5 integer nodes with four relatives.

R EMARK 1. Proposition 5 is stated in terms of the proposed formulation (9) for simplicity. However, we note that the same proof can be used for any correct formulation including the logic (7). In particular, using the big-M formulation (3) and big-M constraints (8) would have the same worst-case complexity (provided an adequate value of M can be found). The second and third assumption in Proposition 5 are imposed for simplicity, but the polynomial time complexity holds for most reasonable branching and variable selection strategies. The first assumption—relaxations being solved to optimality—on the other hand is crucial as the proof requires infeasible nodes to be identified immediately.

■

Meng, Gómez and Mazumder: Computation of LTS: BnB and Hyperplane Arrangements Article submitted to Operations Research

Proposition 5 suggests that formulation (9) is particularly effective in low-dimensional settings, even if the number of datapoints (and binary variables) is large. Moreover, Proposition 5 is concerned with a worst-case performance when no bounding occurs at all. In practice, branchand-bound algorithms perform much better than their worst-case complexity, especially when using strong formulations. Thus, even if computational costs of the order O(np+1 ) may appear prohibitively large for p ≈ 3 and moderately large values of n, we show in our computational experiments that the proposed method explores far fewer branch-and-bound nodes in practice and exhibits good computational performance in practice. Formulation (9) can be used with most off-the-shelf MIO solvers by reformulating the per2 spective terms wi± /zi± as rotated cone constraints (e.g., see Aktürk et al. 2009, Günlük and Linderoth 2010) and linearizing the absolute value constraint (9b). We found in our computations that this direct approach (with Gurobi as the MIO solver) works quite well provided that n is small (in the low hundreds). However, we also observe that the direct approach can struggle if n ≥ 1000, even if p is very small. Indeed, formulation (9) requires creating 4n variables and sets of linear and conic constraints. Solving a large number of convex optimization problems of this form requires a substantial computational effort. In Section 4 we describe a custom branch-and-bound framework to address this issue. Rather than dealing with additional 4n variables, we handle the node relaxation problems only in the space of the regression coefficients β ∈ Rp , which is substantially smaller in dimension. This tailored approach avoids the computational overhead incurred by the auxiliary variables introduced in formulation (9), yielding significant practical efficiency gains, as demonstrated in Section 5. R EMARK 2. Huchette and Vielma (2019) propose, in a different context, MIO formulations for settings defined by hyperplane arrangements. However, their approach is not comparable with ours, nor can it be directly applied to the setting we consider in a practical manner. At a high level, they propose to consider every region defined by the arrangement and construct an ideal formulation for the complete problem. Such an approach would require enumeration of every region a priori (which is computationally demanding for p ≥ 3), and may lead to a prohibitively large number of regions: for example, in a dataset we consider in our experiments with n = 86 and p = 7, the bound from Proposition 2 suggests we may need to enumerate 109 regions and then solve a convex problem with billions of decision variables. In contrast, formulation (9) does not lead to a strengthening of the continuous relaxation (Proposition 4), but instead ensures that enumeration of the regions of the arrangement occurs implicitly during branching. Thanks to the ability to bound, the number

Meng, Gómez and Mazumder: Computation of LTS: BnB and Hyperplane Arrangements Article submitted to Operations Research

of regions considered is substantially less than the worst case scenario: for example, in the same instance with n = 86 and p = 7, Gurobi terminates after only 15, 000 nodes. ■

4. A Custom Branch-and-Bound Algorithm based on First Order Methods In this section, we present a custom nonlinear branch-and-bound (BnB) framework for Problem (9), in which the the primal and dual bounds of the node relaxations are obtained via specialized first-order methods. Our work draws inspiration from the work of Hazimeh et al. (2021) for sparse linear models, but has important differences: Hazimeh et al. (2021) study ridge regularized least squares problem with a penalty on the number of nonzero regression coefficients. They employ standard binary branching on the binary variables (indicating whether a variable is zero or not), whereas we consider a nonconvex optimization problem expressed directly in terms of β ∈ Rp — having projected out all other variables — and adopt a three-way branching scheme tailored to the structure of our formulation. For additional details on the general BnB framework, we refer the reader to Wolsey (2020, Chapter 7); details on the specific components of our customized procedure are discussed in Section 4.2. Formally, defining function f : Rp → R as X n 1 ⊤ ⊤ e di (wi+ )2 di (wi− )2 1 β 2 ⊤ β w Σd + ∥y∥2 − y (Xβ − w) + + µzi + f (β) = min w,z w 2 2 zi+ zi− ± ± i=1 def

s.t. (9a), (9b), (9c) w, w− , w+ ∈ Rn+ , z, z − , z + ∈ {0, 1}n ,

we observe that problem (9) is equivalent to minβ∈Rp f (β). Function β 7→ f (β) is nonconvex, and our branch-and-bound procedure needs to compute, at each node, a lower bound on f restricted to the subregion defined by the node’s branching decisions. In particular, these branching decisions can be described in terms of bounds ℓ± , u± ∈ {0, 1}n with ℓ± ≤ u± on the binary variables z ± . Given such bounds, we consider lower bounding functions inspired by the interval relaxation def

X n 1 ⊤ ⊤ e di (wi+ )2 di (wi− )2 1 β β w Σd ∥y∥22 − y ⊤ (Xβ − w) + + µzi + + w 2 2 zi+ zi− ± i=1

fR (β; ℓ± , u± ) = min w,z w± ,z

s.t. (9a), (9b), (9c) ℓ− ≤ z − ≤ u− , ℓ+ ≤ z + ≤ u+ w, w− , w+ ∈ Rn+ , z, z − , z + ∈ [0, 1]n .

(11)

First, in Section 4.1, we describe explicit forms and relaxations of fR (β; ℓ± , u± ) based on augmented Lagrangians. Next, in Section 4.2, we describe the overall branching scheme, pruning conditions and other details for the BnB algorithm. In Section 4.3, we discuss how to solve

Meng, Gómez and Mazumder: Computation of LTS: BnB and Hyperplane Arrangements Article submitted to Operations Research

node relaxations and efficiently compute dual bounds. Finally, in Section 4.4 we discuss additional details and computational improvements of the BnB algorithm.

4.1. Closed form expressions and differentiable lower bounds for fR We begin by showing that fR admits a separable closed-form representation in terms of onedimensional functions, which forms the basis of our subsequent algorithmic developments. P ROPOSITION 6. For any β ∈ Rp and any bounds ℓ± , u± ∈ {0, 1}n with ℓ± ≤ u± , the relaxation fR defined in (11) can be written as n X λ − + − + 2 fR (β; ℓ , u ) = ∥β∥2 + ϕ yi − x⊤ i β; ℓi , ℓi , ui , ui , di , 2 i=1 ±

(12)

where, for each r ∈ R and d > 0, the one-dimensional function ϕ is defined as ϕ(r; ℓ− , ℓ+ , u− , u+ , d) =

min

w− ,w+ ∈R+ z − ,z + ∈R+

1 2

+ 2

− 2

(r − w+ + w− )2 + µ(z − + z + ) + d(wz+ ) + d(wz− ) − d(w+ − w− )2 p

(13a)

p w+ ≥ 2µ z + , p |r − w+ + w− | ≤ 2µ (1 − z − − z + ),

(13b)

ℓ− ≤ z − ≤ u− ,

(13d)

s.t. w− ≥

2µ z − ,

ℓ+ ≤ z + ≤ u+ ,

z − + z + ≤ 1.

(13c)

Moreover, provided that d satisfies the criteria of Section 2.1.2, the function β 7→ fR (β; ℓ± , u± ) is convex. The representation (12) follows from the observation that, for fixed β, the inner minimization in (11) is separable across the per-datapoint variable groups (wi , zi , wi− , wi+ , zi− , zi+ ): constraints (9a)–(9c) and the bounds on z ± all decouple across i; adding and subtracting w⊤ Diag(d)w to the objective and minimizing out wi in closed form yields the per-datapoint subproblems (13) after the substitution r = yi − x⊤ i β. Convexity of fR is a direct consequence of the discussion in Section 2.1.2. Proposition 7 provides a closed form expression of function ϕ when ℓ± = 0 and u± = 1.

Meng, Gómez and Mazumder: Computation of LTS: BnB and Hyperplane Arrangements Article submitted to Operations Research

P ROPOSITION 7. If 0 < d < 1/2, ℓ± = 0 and u± = 1, then ϕ defined in (13) is given by:  √ 1   r2 , |r| ≤ 2 µd (z − + z + = 0)   2    r  −dr2 + 2√µd|r| − 2µd √ µ ϕ(r; 0, 0, 1, 1, d) = , 2 µd < |r| < , (0 < z − + z + < 1)  1 − 2d d   r   µ   µ, |r| ≥ , (z − + z + = 1) d

(14)

where we indicate in parenthesis the optimal values of variables (z − , z + ). Moreover, the result holds if constraints (13c) are removed from (13). Interestingly, the closed form solution reveals new interpretations of the convex relaxations of the perspective formulation (5) and proposed formulation (9): they are equivalent to using a special family of nonconvex loss functions that underestimate the capped quadratic loss ϕcap . To simplify the notation, we will use the shorthand: def

ϕ̄(r; d) = ϕ(r; 0, 0, 1, 1, d). Observe that the second statement of the proposition, stating that constraints (13c) are redundant, imply Proposition 4: at the root node, when all lower bounds are zero and upper bounds are one, the hyperplane arrangement constraints do not improve the strength of the perspective relaxation. Figure 4 illustrates function r 7→ ϕ̄(r; d) for different values of d. We observe that function r 7→ ϕ̄(r; d) underestimates the capped loss ϕcap and becomes a better approximation of ϕcap as d → p √ 1/2, and the two functions coincide for |r| < 2 µd and |r| ≥ µ/d. The function r 7→ ϕ̄(r; d) is differentiable: we provide explicit forms of its derivative in Appendix A.2. While function r 7→ ϕ̄(r; d) is nonconvex, we point out that fR is convex as stated in Proposition 6 because the presence of the strongly convex term λ2 ∥β∥22 offsets the nonconvexities introduced by ϕ̄. We conclude this section by deriving explicit closed-form expressions for ϕ at the remaining bound configurations. We also construct differentiable lower bounds for ϕ: as the closed-form expressions involve hard indicator constraints on the residual r, which effectively turn the minimization of relaxation (12) into a constrained optimization problem, preventing direct application of gradient descent. To avoid this, we employ the augmented Lagrangian construction (Bertsekas 1982). Concretely, for a constrained problem of the form minr∈R g(r) subject to cj (r) ≤ 0 for j = 1, . . . , m, the corresponding augmented Lagrangian with multipliers ν ∈ Rm + and penalty paramPm P def m eter ρ ≥ 0 is Lρ (r; ν) = g(r) + j=1 νj cj (r) + ρ2 j=1 [cj (r)]2+ , where [ · ]+ = max{·, 0}. Whenever g and the cj ’s are differentiable, so is Lρ (·; ν); moreover, by weak duality, minr Lρ (r; ν) ≤

Meng, Gómez and Mazumder: Computation of LTS: BnB and Hyperplane Arrangements Article submitted to Operations Research

Comparison of cap with = 0.5 and its relaxations

0.5

Loss

0.4

Big-M (M=20) (d=0.003) (d=0.01) (d=0.03) Capped quadratic

0.3 0.2 0.1 0.0

Figure 4

0 r

Comparison of loss functions for a single residual r: the original least trimmed loss ϕcap (before relaxation), the perspective relaxation ϕ̄(r; d), and the analogous one-dimensional function obtained by applying the same projection argument to the standard Big-M relaxation of (9) in place of the perspective relaxation. The perspective relaxation provides a tighter convex approximation than the Big-M relaxation.

minr {g(r) : cj (r) ≤ 0, ∀j} for any ν ≥ 0 and ρ ≥ 0, so that Lρ (·; ν) serves as a differentiable lower-bounding function for the constrained problem. Applying this construction to each of the constrained forms of ϕ that arise during branching yields the following proposition. P ROPOSITION 8. Let δS (·) denote the indicator function of a set S, taking value 0 on S and +∞ elsewhere. For any d > 0 and any ν, ν1 , ν2 , ρ ≥ 0, the function ϕ defined in (13) admits the following closed-form expressions, each accompanied by a differentiable lower bound (in r) obtained from the augmented Lagrangian construction described above: • ϕ(r; 0, 0, 0, 0, d) = 12 r2 + δ[−√2µ,√2µ] (r) and a differentiable lower bound is given by p p ρ h p i2 h p i2 1 2 ϕ̄0 (r; ν1 , ν2 , ρ) = r + ν1 −r − 2µ + ν2 r − 2µ + −r − 2µ + r − 2µ . 2 2 + + • ϕ(r; 1, 0, 1, 0, d) = µ + δ[−∞,−√2µ] (r) and a differentiable lower bound is given by p ρh p i2 ϕ̄− (r; ν, ρ) = µ + ν r + 2µ + r + 2µ . 2 + • ϕ(r; 0, 1, 0, 1, d) = µ + δ[√2µ,∞] (r) and a differentiable lower bound is given by p ρh p i2 ϕ̄+ (r; ν, ρ) = µ + ν −r + 2µ + −r + 2µ . 2 + The functions ϕ̄0 , ϕ̄− , ϕ̄+ play a key role in the remainder of the paper: in Section 4.2, they are used to construct the augmented Lagrangian of the full node relaxation (12), and in Section 4.3 this Lagrangian is in turn the objective minimized by our augmented Lagrangian method (ALM) for computing dual bounds at each BnB node.

Meng, Gómez and Mazumder: Computation of LTS: BnB and Hyperplane Arrangements Article submitted to Operations Research

4.2. Description of the BnB method We now describe the main components of our custom branch-and-bound algorithm. Root node The BnB process begins by solving the interval relaxation of problem (9) at the root node. From (12) we find that the relaxation is given by n X λ 2 minp ∥β∥2 + ϕ̄(yi − x⊤ i β; di ), β∈R 2 i=1

(15)

where function ϕ̄ is explictly described in Proposition 7. All functions involved are differentiable, thus (15) can be solved via standard first-order methods. In our approach, we solve it approximately via gradient descent with line search. We describe our specific implementation and the computation of valid dual bounds from approximate solutions in Section 4.3. Branching variable selection Once a root (or node) relaxation has been solved, typical branchand-bound algorithms for mixed-integer optimization would then select a fractional integer variable to branch on. In our case, the binary variables are not defined explicitly, but appear implicitly in the definition of function ϕ in (13). In particular, function ϕ̄(·; d) underestimates the capped loss ϕcap (·), and is a strict underestimator if and only if the binary variables appearing in (13) are fractional. Thus, given an optimal solution β̄ ∈ Rp to Problem (15), the algorithm selects an index ⊤ i ∈ [n] such that ϕ̄(yi − x⊤ i β̄; di ) < ϕcap (yi − xi β̄).

Branching scheme The definition of ϕ in Problem (13) involves two binary variables (z − and z + ) and three feasible configurations (z − = z + = 0; z − = 1, z + = 0; and z − = 0, z + = 1; note that z − = z + = 1 is infeasible). Thus, instead of using the standard variable dichotomy for branching, we use three-way branching for our custom branch-and-bound. After fixing the branching variables, the child

Disjunction

Replace ϕ̄(yi − x⊤ i β; di ) with...

node’s relaxation objective is obtained from

zi− = 0, zi+ = 0

ϕ̄0 (·; ν1 , ν2 , ρ)

the parent’s by replacing ϕ̄(yi − x⊤ i β̄; di ) in

zi− = 1, zi+ = 0

ϕ̄− (·; ν, ρ)

the parent node’s objective with the appropri-

zi− = 0, zi+ = 1

ϕ̄+ (·; ν, ρ)

ate differentiable lower bounds described in

Table 1: Subproblem construction

Proposition 8, as summarized in Table 1. Recursion The algorithm proceeds recursively, applying the same branching scheme to each newly created node. To describe the relaxation solved at a generic node, let F0 , F− , F+ ⊆ [n] denote the disjoint sets of indices recording the branching decisions made along the path from the root to the current node: i ∈ F0 if zi− = zi+ = 0 has been fixed, i ∈ F− if zi− = 1 has been fixed, and i ∈ F+

Meng, Gómez and Mazumder: Computation of LTS: BnB and Hyperplane Arrangements Article submitted to Operations Research

if zi+ = 1 has been fixed. Following the construction summarized in Table 1, the term ϕ̄(·; di ) associated with each branched index i is replaced by the corresponding differentiable lower bound from Proposition 8, parameterized by dual multipliers ν and a penalty parameter ρ ≥ 0. This yields the node objective X X def λ ϕ̄− (yi − x⊤ Lρ (β, ν) = ∥β∥22 + β; (ν ) , (ν ) , ρ) + ϕ̄0 (yi − x⊤ i 1 i 2 i β; νi , ρ) i 2 i∈F− i∈F0 X X ϕ̄+ (yi − x⊤ ϕ̄(yi − x⊤ + i β; νi , ρ) + i β; di ). i∈F+

2|F |+|F− |+|F+ |

where ν ∈ R+ 0

(16)

i̸∈(F0 ∪F− ∪F+ )

collects all dual multipliers associated with the branched indices. By

Proposition 8, each function ϕ̄0 , ϕ̄− , ϕ̄+ is precisely the augmented Lagrangian of the corresponding constrained function ϕ. Consequently, Lρ (β, ν) is the augmented Lagrangian of (12) in which the bounds (ℓ± , u± ) encode the branching decisions in F0 , F− , F+ . Since the objective in (12) is convex (in β) and all constraints are linear, standard augmented Lagrangian duality (Bertsekas 1982, Rockafellar 1976) guarantees that ζ̄ = max min Lρ (β, ν), ν≥0

(17)

where ζ̄ denotes the optimal value of the node relaxation (12). Solving the node thus reduces to (approximately) solving this max-min problem; the procedure is described in detail in the next subsection, and produces a primal iterate β̄ together with a valid lower bound ζ̄ ′ on ζ̄. The node is then processed according to the standard BnB rules: if ζ̄ ′ exceeds the current incumbent, the ⊤ node is pruned by bounding; if ϕ̄(yi − x⊤ i β̄; di ) = ϕcap (yi − xi β̄) for all i ∈ [n] \ (F0 ∪ F− ∪ F+ ),

the node is pruned by integrality and the incumbent is updated (in our implementation, we adopt a best-bound strategy, so this case arises only at the final node explored, if at all); otherwise, an ⊤ index i ∈ [n] \ (F0 ∪ F− ∪ F+ ) with ϕ̄(yi − x⊤ i β̄; di ) < ϕcap (yi − xi β̄) is selected for branching,

and the recursion continues. R EMARK 3. Since our BnB implementation does not solve the node relaxation problem to optimality (due to the nature of the first order methods which can be slow to obtain high accuracy solutions), the first condition of Proposition 5 is no longer satisfied. As a consequence, the proposed algorithm is no longer guaranteed to explore at most O(np+1 ) nodes. In practice however, when optimally solving the convex interval relaxation at a node becomes expensive, our approach

Meng, Gómez and Mazumder: Computation of LTS: BnB and Hyperplane Arrangements Article submitted to Operations Research

is typically still able to prune the node by bounding, so the total number of nodes explored remains largely unaffected. This is confirmed empirically in Tables 8–9 of Section 5, which compare the node counts of our solver against those of Gurobi applied directly to Problem (9). 4.3. Node relaxations and dual bounds At each node of the branch-and-bound algorithm, we compute dual bounds for Problem (17). Observe that for any fixed ν̄

ζ̄(ν̄) = min Lρ (β; ν̄)

(18)

is a lower bound on ζ̄. Function Lρ (β; ν) is differentiable as a function of β, and the gradient is

∇β Lρ (β, ν) = λβ −

xi ϕ̄′0 (yi − x⊤ i β; (νi )1 , (νi )2 , ρ) −

i∈F0

−

xi ϕ̄′− (yi − x⊤ i β; νi , ρ)

i∈F−

xi ϕ̄′+ (yi − x⊤ i β; νi , ρ) −

i∈F+

xi ϕ̄′ (yi − x⊤ i β; di )

(19)

i̸∈(F0 ∪F− ∪F+ )

where functions ϕ̄′ ϕ̄′0 , ϕ̄′− , ϕ̄′+ denote the derivatives, which admit a closed-form expression and are continuous with respect to β (see Appendix A.2). In the proposed approach, we obtain a solution to (18) via gradient descent with line search. Note that ζ̄(ν̄) is a lower bound only if the associated problem is solved to optimality. However, first order methods can be slow to converge to highprecision solutions. Therefore, we need to consider methods to get reliable lower bounds from an approximate solution, as we discuss next. We recall that for any λ̃-strongly convex function β 7→ Φ(β) where Φ(β) is differentiable, the optimal value satisfies min Φ(β) ≥ Φ(β̄) − β

1 ∥∇Φ(β̄)∥22 2λ̃

(20)

for any point β̄. As mentioned in Section 2.1.2, the perspective reformulation and relaxations can be ensured to be λ̃-strongly convex for any λ̃ ∈ (0, λ), see (6) (note that larger values of the strong convexity parameter λ̃ come at the expense of weaker convex relaxations and thus less effective pruning of the BnB algorithm). Using this observation, we obtain the following proposition.

Meng, Gómez and Mazumder: Computation of LTS: BnB and Hyperplane Arrangements Article submitted to Operations Research

P ROPOSITION 9. Let (β̄, ν̄) be any primal-dual pair of Problem (17), then the inequality ζ̄(ν̄) ≥ Lρ (β̄, ν̄) −

∥∇β Lρ (β̄, ν̄)∥22 2λ̃

(21)

holds true. The proof of Proposition 9 is provided in Appendix A.3. Using Proposition 9, we can compute dual bounds for the convex relaxations in the BnB algorithm efficiently from any iterate of a firstorder method and its gradient, even when the relaxation is solved approximately. Building on this result, we solve Problem (17) via the augmented Lagrangian method, which alternates between minimizing Lρ (β, ν) over the primal variables β with multipliers ν fixed, and updating ν. In the primal step, we minimize Lρ (β, ν) over β using gradient descent with an Armijo backtracking line search, as detailed in Algorithm 1. The complete augmented Lagrangian procedure for solving (17) is summarized in Algorithm 2. Here, Proposition 9 is applied at each iteration of both the inner and outer loops to obtain a valid dual bound from the current iterate. As the multipliers ν are progressively updated through the outer loop, the quality of the resulting dual bounds improves. Algorithm 1 Gradient descent with Armijo backtracking for minimizing Lρ (·, ν) Input: Initial point β, dual values ν, initial step size t0 , Armijo parameter α ∈ (0, 1), backtracking factor γ ∈ (0, 1), tolerance ϵ 1: t ← t0 2: repeat 3:

Compute g = ∇β Lρ (β, ν)

while Lρ (β − tg, ν) > Lρ (β, ν) − αt∥g∥22 do

t ← γt

end while

β ← β − tg

t ← t/γ

Compute dual bound D using Proposition 9

10: until (Lρ (β, ν) − D)/Lρ (β, ν) ≤ ϵ 11: return β

▷ warm start for next iteration

Meng, Gómez and Mazumder: Computation of LTS: BnB and Hyperplane Arrangements Article submitted to Operations Research

Algorithm 2 Augmented Lagrangian Method for solving the relaxation problem Input: Initial solution β (0) , multipliers ν (0) ≥ 0, penalty ρ > 0, tolerances ϵin , ϵout 1: for k = 0, 1, 2, . . . do 2:

β (k+1) ← minimize Lρ (β, ν (k) ) using Algorithm 1 with tolerance ϵin , starting from β (k)

Update dual multipliers: νj = max{0, νj + ρ cj (β (k+1) )} for all j ∈ C n √ o (k+1) (k) ⊤ (k+1) + 2µ for all i ∈ F− νi = max 0, νi + ρ yi − xi β n √ o (k+1) (k) (k+1) 2µ for all i ∈ F+ νi = max 0, νi + ρ −yi + x⊤ β + i n √ o (k+1) (k) (k+1) β − (νi )1 = max 0, (νi )1 + ρ −yi + x⊤ 2µ for all i ∈ F0 i n √ o (k+1) (k) (k+1) (νi )2 = max 0, (νi )2 + ρ yi − x⊤ β − 2µ for all i ∈ F0 i

Compute dual bound D(k+1) using Proposition 9

if (D(k+1) − D(k) )/D(k) ≤ ϵout then

4: 5: 6:

10: 11:

(k+1)

(k)

break end if

12: end for 13: return β (k+1) , dual bound D (k+1)

4.4. Implementation Details We discuss several practical techniques to make our proposed BnB approach computationally efficient. Additionally, we discuss how our work differs from L0BnB (Hazimeh et al. 2021), who also use a custom BnB framework with fist order methods to solve the node relaxations for the sparse linear regression which is different from the robust statistics problem we consider here. BnB configuration For node selection, we employ a best-bound search strategy (Linderoth and Savelsbergh 1999), which selects the open node with the smallest lower bound as the next node to explore. For selecting the branching variable, we branch on the most fractional variable. Note that the integer variables z − , z + are not explicitly defined, but can be computed implicitly: choosing the index i whose zi− + zi+ is closest to 1/2 is equivalent to choosing r µ i⋆ ∈ arg min | |ri | − r0 | , where r0 := 12 + d , ri = yi − x⊤ i β. i d p Equivalently, since ϕ(r) is strictly increasing in |r| when |r| ∈ [0, µ/d], we find µ 34 − d − d2 ⋆ i ∈ arg min | ϕ(ri ) − ϕ0 | , ϕ0 := ϕ(r0 ) = . i 1 − 2d

Meng, Gómez and Mazumder: Computation of LTS: BnB and Hyperplane Arrangements Article submitted to Operations Research

The BnB procedure terminates when the gap between bounds falls below a predetermined tolerance (1% in our experiments), or when all nodes have been explored or pruned. Hyperparameter selection In our implementation, we set the quadratic penalty term ρ to 100 if n ≥ 300 and equal to 5 otherwise. The strong-convexity parameter is set to λ̃ = 0.1λ. We use Armijo type line search with initial step size t0 = 1, sufficient-decrease parameter α = 10−4 , and backtracking factor γ = 0.5. We set tolerances ϵin = 10−4 and ϵout = 10−4 in Algorithm 2. We found the method to be fairly robust to the choice of hyperparameters: ρ was selected from {1, 5, 10, 50, 100, 500} on 10 synthetic datasets of various sizes; the remaining hyperparameters were not tuned. Warm starting Rather than solving the node relaxation (17) from scratch, we warm start Algorithm 2 using the parent node’s solution. Specifically, β (0) is initialized to the final primal iterate (0)

of the parent node. For each dual multiplier νi , if the corresponding constraint is also present (0)

in the parent node, we inherit its final dual value; otherwise, we set νi

= 0. Since a child node

differs from its parent only by fixing a single binary variable zi , the parent’s solution typically provides a high-quality starting point, substantially reducing the number of iterations required for convergence in practice. Selective upper bound computation At any node of the BnB tree, we can use heuristics similar to those discussed in Section 2.3 based on alternating minimization to compute upper bounds for Problem (9). We provide a precise description of the heuristics in Appendix A.4. As computing incumbent solutions at every node is expensive, we compute these upper bounds only at the root node and nodes whose depth is a multiple of a constant (10 in our experiments). This reduces computational overhead and cost while maintaining effective pruning. Early pruning. Since dual bounds via Proposition 9 are already computed at every iteration of Algorithms 1 and 2, we can monitor them on the fly: as soon as the dual bound exceeds the current incumbent, we terminate the inner solver and prune the node immediately. This step is critical for the efficiency of our BnB procedure, as the original relaxation (12) at a given node may be infeasible, in which case Algorithm 2 produces a sequence of dual bounds D(k) diverging to infinity. Without early pruning, the algorithm would expend substantial effort driving D(k) toward infinity on a node that can already be safely discarded. We note that the use of heuristics at the root node guarantees that an incumbent solution is always available for this comparison. Accelerated gradient computation We use Numba (Lam et al. 2015), a just-in-time compiler for Python, to accelerate the gradient evaluation in (19).

Meng, Gómez and Mazumder: Computation of LTS: BnB and Hyperplane Arrangements Article submitted to Operations Research

Relation to L0BnB At a high level, our proposed BnB algorithm extends L0BnB (Hazimeh et al. 2021)—a BnB framework originally designed for ℓ0 -penalized least squares regression—to the outlier detection problem (1). Several key modifications are needed to address the distinct challenges posed by (1). First, we build our BnB on an enhanced MIO formulation (9) and derive a new form of the continuous relaxation (Section 4.1). Second, we develop a three-way branching scheme that efficiently encodes the status of each observation—inlier or outlier, with positive or negative residual (Section 4.2). Third, the node relaxations appearing in L0BnB are in the composite form (without constraints) and coordinate descent procedures are used to compute solutions to these problems. In this work, the relaxation subproblems are more complicated as they have constraints, and we use augmented Lagrangian approaches and also introduce a new procedure for extracting valid dual bounds (Section 4.3).

5. Experiments We evaluate our proposed BnB solver against commercial MIO solvers on both synthetic and real datasets. Our experiments focus on three aspects: computational efficiency compared to existing methods, sensitivity to problem parameters, and the value of exact optimization over heuristic approaches. 5.1. Experimental Setup Solver settings: We implement our BnB solver in Python with gradient computations accelerated using Numba (Lam et al. 2015). We compare against Gurobi (Gurobi Optimization, LLC 2022) and Mosek (ApS 2022) applied to the Big-M formulation (3) (M), the perspective formulation (5) (P), and the strengthened perspective formulation with hyperplane-arrangement (9) (H). Accordingly, BnB (P) and BnB (H) denotes our solver on perspective formulation and hyperplane-arrangement formulation, respectively. And Gurobi/Mosek (M, P, H) denotes the commercial solvers on the same formulations. For all solvers, we define the relative optimality gap as (U B −LB)/U B, where U B is the incumbent objective value and LB is the best dual bound. We terminate when the gap falls below 1% or when the runtime exceeds 1 hour. All experiments are conducted on a computing cluster, utilizing an AMD EPYC 9474F machine with 10 CPU cores and 40GB RAM. Synthetic data generation: Following Hazimeh and Mazumder (2020), we generate the data matrix X ∈ Rn×p by drawing each row from a multivariate Gaussian N (0, Ip ). The true coefficient vector β † ∈ Rp is generated with entries drawn uniformly from [0, 1]. We generate clean responses

Meng, Gómez and Mazumder: Computation of LTS: BnB and Hyperplane Arrangements Article submitted to Operations Research

as yclean = Xβ † + ε, where ε ∼ N (0, σ 2 In ) and σ is chosen to achieve a target signal-to-noise ratio SNR = Var(Xβ † )/σ 2 . To introduce outliers, we contaminate a subset of observations by adding heavy-tailed noise to their responses. Specifically, for each contaminated observation i, we set yi = yclean,i + δ · σy · Zi , where σy = std(yclean ), δ controls the outlier magnitude, and Zi is drawn from a t-distribution with 3 degrees of freedom. Unless otherwise specified, we set SNR = 50, δ = 10, and introduce 10 outliers in each dataset. Real datasets. We evaluate solver performance on 13 benchmark regression datasets: 7 smaller datasets that have been used in prior work on the LTS problem (Gómez and Neto 2025) and 6 larger datasets obtained from the OpenML repository (Vanschoren et al. 2014). The dimensions of each dataset are summarized in the headers of the corresponding results tables. Parameter selection: We set the ridge regularization parameter as λ = λ0 · Mean(Diag(X ⊤ X)), where λ0 = 0.01 for synthetic datasets and λ0 = 0.2 for real datasets, unless otherwise specified. The outlier penalty µ is chosen so that the model identifies exactly the number of outliers present in the synthetic dataset (or the specified number for real datasets). 5.2. Comparison on Synthetic Datasets 5.2.1. Varying dataset size Table 2 compares runtime performance across different numbers of samples n ∈ {1000, 2000, 5000} and features p ∈ {10, 20, 50}. The results reveal several key findings. n = 1000

Method p=10 *BnB (P) *BnB (H)

p=20

n = 2000 p=50

p=10

p=20

n = 5000 p=50

p=10

p=20

p=50

(10.0%, 0) (8.8%, 1) (7.0%, 2) (14.6%, 0) (14.3%, 0) (13.8%, 0) (17.5%, 0) (17.7%, 0) (18.4%, 0) 2.5 3.6 12.9 6.3 10.0 53.4 35.7 49.2 750.6

Gurobi (M) (33.3%, 0) (31.5%, 0) (25.6%, 0) (53.8%, 0) (49.9%, 0) (51.0%, 0) (65.0%, 0) (64.6%, 0) (64.2%, 0) Gurobi (P) (5.3%, 2) (4.9%, 3) (4.5%, 6) (13.7%, 0) (14.1%, 0) (15.3%, 0) (19.2%, 0) (18.8%, 0) (19.4%, 0) 414.0 766.4 (1.5%, 8) 2301.2 (6.6%, 5) (14.0%, 0) (16.1%, 0) (19.6%, 0) *Gurobi (H) 358.5 Mosek (M) (32.3%, 0) (30.1%, 0) (16.9%, 0) (47.4%, 0) (47.0%, 0) (45.4%, 0) (58.9%, 0) (58.6%, 0) (58.8%, 0) Mosek (P) (7.5%, 0) (7.5%, 0) (4.7%, 1) (16.1%, 0) (15.4%, 0) (15.0%, 0) (19.4%, 0) (20.0%, 0) (20.9%, 0) *Mosek (H) (8.6%, 0) (1.0%, 2) (1.8%, 0) (12.5%, 0) (11.0%, 0) (18.9%, 0) (20.1%, 0) (20.1%, 0) (21.3%, 0) Table 2

Runtime performance on synthetic datasets, averaged over 10 random seeds. Methods marked with * are proposed in

this work. Methods: (M) Big-M, (P) perspective, (H) hyperplane arrangement formulation. Numeric values indicate solve time in seconds when all 10 instances are solved within 1 hour. Otherwise, we report (gap, count) where gap is the average optimality gap and count is the number of instances solved.

Meng, Gómez and Mazumder: Computation of LTS: BnB and Hyperplane Arrangements Article submitted to Operations Research

Benchmarking methods the Big-M formulation (M) yields the weakest relaxation bounds, with optimality gaps exceeding 30% for smaller instances and over 60% for larger ones. The perspective formulation (P) significantly improves upon Big-M, reducing gaps to 5–20% depending on problem size – these results are consisting with those reported by Gómez and Neto (2025). Hyperplane arrangement with off-the-shelf solvers The hyperplane arrangement formulation (H) allows for more aggresive pruning of the branch-and-bound tree. When used with Gurobi, problems with n = 1000 can now be solved to optimality within minutes, while just using the perspective reformulation results in 5% gaps. Gurobi is also able to solve more instances and reduce optimality gaps in problems with n = 2, 000, but the benefits seem to disappear in problems with n = 5, 000. The improvements are far less notable with Mosek, and in many cases the use of hyperplane arrangement actually degrades performance. These findings suggest that Gurobi is better equipped to handling problems with heavy constraints. Hyperplane arrangement with the proposed BnB solver The proposed BnB (H) consistently outperforms all commercial solvers, achieving up to 100× speedup over Gurobi (H) on comparable instances. This efficiency stems from our reformulation of the relaxation problem in the β-space, which involves only p variables and has low complexity dependence on n. For problems with n = 1000, the method solves instances in seconds while Gurobi requires several minutes. For the largest instances (n = 5000, p = 50), BnB (H) completes in approximately 12 minutes, whereas competing methods fail to achieve the 1% optimality gap within the time limit. Notably, even without the additional hyperplane arrangement constraints, BnB (P) achieves comparable performance to Gurobi (P) for large n, confirming the method’s scalability with sample size. Additional results with smaller datasets We also examine solver performance on smaller datasets. To maintain problem difficulty and clearly distinguish solver performance, we set SNR = n/25 and δ ∈ {5, 6, 9} for n ∈ {50, 100, 200}, respectively. Table 3 presents the results. n = 50

Method

n = 100

n = 200

p=5 p=10 p=20 p=5 p=10 p=20 p=5 p=10 p=20 Gurobi (P) *Gurobi (H) *BnB (H) Table 3

2.5 2.0 2.0

2.2 2.8 3.9

4.2 7.1 34.5

21.9 7.5 10.5

26.3 23.3 12.5

19.2 33.3 11.9

447.0 27.2 9.5

87.5 39.8 9.0

97.0 77.9 50.0

Runtime performance on small datasets, averaged over 10 random seeds. Methods marked with * are proposed in this work. Numeric values indicate solve time in seconds.

From Table 3, we observe that Gurobi (H) is highly efficient when the dataset scale is very small (n = 50), solving instances in 2–7 seconds and outperforming BnB (H). However, Gurobi’s

Meng, Gómez and Mazumder: Computation of LTS: BnB and Hyperplane Arrangements Article submitted to Operations Research

runtime increases as n grows. For n ≥ 100, BnB (H) achieves better performance in most of the settings considered. 5.2.2. Ablation studies We investigate how various problem parameters affect solver performance. We fix n = 2000 and p = 20, adopt the parameter settings described in Section 5.1, and vary one parameter at a time to conduct the ablation. Table 4 examines the impact of regularization strength λ0 (the ridge is defined as λ = λ0 · Mean(Diag(X ⊤ X))) and the number of outliers, while Table 5 analyzes the effect of signal-to-noise ratio and outlier magnitude δ. λ0

Gurobi (P) *Gurobi (H) *BnB (H)

0.0005 0.001 0.002 0.005 0.01 0.02 0.05

(66.3%, 0) (57.9%, 0) (46.4%, 0) (27.8%, 0) (14.3%, 0) (4.9%, 0) 453.3

Table 4

(61.7%, 0) (53.6%, 0) (43.4%, 0) (18.2%, 1) (1.2%, 8) 1627.2 731.6

(32.4%, 0) (19.1%, 0) (5.2%, 4) 34.7 10.3 4.3 0.9

# outliers Gurobi (P) *Gurobi (H) *BnB (H) 5 10 15 20 25 30

(9.9%, 3) (14.2%, 0) (15.5%, 0) (18.5%, 0) (16.4%, 0) (17.5%, 0)

(1.1%, 9) (1.0%, 9) (5.7%, 5) (10.3%, 4) (8.8%, 5) (12.3%, 4)

7.3 9.2 (2.8%, 9) (4.7%, 8) (4.2%, 8) (4.4%, 6)

Impact of regularization strength (left) and number of outliers (right). Numeric values indicate solve time in seconds

when all 10 instances are solved within 1 hour. Otherwise, we report (gap, count) where gap is the average optimality gap and count is the number of instances solved.

From Table 4, we observe that stronger regularization (λ0 ≥ 0.01) significantly improves tractability for all methods, while weak regularization (λ0 ≤ 0.002) makes the problem harder due to weaker relaxations. The table also shows that problem difficulty increases with the number of outliers. The BnB solver handles up to 10 outliers efficiently, solving instances in under 10 seconds. As the outlier count increases beyond 15, performance degrades for all methods, though BnB (H) consistently achieves smaller optimality gaps. SNR

Gurobi (P) *Gurobi (H) *BnB (H)

5.0 10.0 20.0 30.0 40.0 50.0 100.0

(52.7%, 0) (45.7%, 0) (32.7%, 0) (24.1%, 0) (18.4%, 0) (14.3%, 0) (2.3%, 4)

Table 5

(52.2%, 0) (44.8%, 0) (33.6%, 0) (18.5%, 1) (1.0%, 9) (1.2%, 8) 1268.1

(45.0%, 0) (34.8%, 0) (13.1%, 1) (1.3%, 8) 18.9 9.4 3.6

Gurobi (P) *Gurobi (H) *BnB (H)

2.0 4.0 6.0 8.0 10.0 15.0

(40.8%, 0) (36.9%, 0) (29.2%, 0) (22.0%, 0) (14.3%, 0) (2.4%, 5)

(40.7%, 0) (36.7%, 0) (29.4%, 0) (11.3%, 2) (1.7%, 8) 951.2

(36.9%, 0) (31.0%, 0) (15.5%, 0) 474.4 9.3 3.3

Impact of signal-to-noise ratio (left) and outlier magnitude (right). Numeric values indicate solve time in seconds

when all 10 instances are solved within 1 hour. Otherwise, we report (gap, count) where gap is the average optimality gap and count is the number of instances solved.

Meng, Gómez and Mazumder: Computation of LTS: BnB and Hyperplane Arrangements Article submitted to Operations Research

Table 5 demonstrates that both higher SNR and larger outlier magnitudes make outliers easier to identify, thereby improving solver performance. When SNR ≤ 20 or δ ≤ 6, distinguishing outliers from inliers becomes difficult, resulting in large optimality gaps for all methods. The proposed solver achieves the best performance across all parameter settings, solving most instances to optimality when SNR ≥ 30 or δ ≥ 8. 5.3. Evaluation on Real Datasets We report solver performance on 13 benchmark regression datasets. For each dataset, we set the outlier penalty µ so that the model identifies exactly 10 outliers. Table 6-7 presents the runtime and optimality gap results, and Table 8-9 reports the number of BnB nodes explored by each solver. The results are summarized in Figure 5. Method

alcohol education food (44, 6) (50, 4) (150, 3)

*BnB (P) (28.6%) *BnB (H) 3.4 Gurobi (M) 11.9 Gurobi (P) 25.4 *Gurobi (H) 5.4 Mosek (M) 127.7 Mosek (P) 21.4 *Mosek (H) (24.8%) Table 6

(14.6%) 118.3 754.6 38.5 13.6 (8.8%) 235.4 (–%)

milk (86, 7)

pulpfiber radar wagner (62, 7) (1573, 4) (63, 6)

(7.6%) (10.9%) (29.0%) 1.0 128.9 2.0 (18.5%) 1564.5 16.6 39.2 51.8 7.6 3.0 28.8 10.6 (26.0%) (11.2%) 129.8 746.2 563.2 8.4 (–%) (12.9%) (26.2%)

(6.7%) 23.7 (23.5%) (4.4%) 660.2 (23.6%) (7.7%) (–%)

(14.7%) 56.2 80.0 19.3 12.9 1234.2 101.8 (12.9%)

Solver performance on smaller real benchmark datasets with 10 outliers. Methods marked with * are proposed in this

work. Dataset dimensions (n, p) are shown below each name. Values in parentheses indicate the optimality gap when the solver fails to reach the 1% tolerance within the time limit; numeric values indicate solve time in seconds. Bold entries highlight the best result for each dataset. (–%) indicates solver failure.

We observe that the hyperplane arrangement formulations (H) with Gurobi or the proposed BnB algorithm consistently outperform the Big-M and perspective formulations across all datasets, demonstrating the value of incorporating threshold constraints (similarly to results with synthetic data, hyperplane arrangement formulations with Mosek perform poorly). Neither BnB (H) or Gurobi (H) consistently outperforms the other: from the performance profile in Figure 5, BnB (H) seems to be able to prove optimality faster in easier instances (less than 10 seconds), and Gurobi (H) seems to be more effective in instances requiring more than 1,000 seconds. Tables 8–9 reveal that the hyperplane arrangement formulation (H) requires fewer nodes to solve to optimality compared to the perspective formulation (P), but each node takes longer to solve. This is evident from instances where all solvers reach the time limit: the perspective formulation

Meng, Gómez and Mazumder: Computation of LTS: BnB and Hyperplane Arrangements Article submitted to Operations Research

Table 7

Method

pm10 pollen no fri rmftsa galaxy (500, 7) (3848, 4) (500, 7) (1000, 50) (508, 10) (323, 4)

*BnB (P) *BnB (H) Gurobi (M) Gurobi (P) *Gurobi (H) Mosek (M) Mosek (P) *Mosek (H)

(23.8%) (19.2%) (15.0%) (15.6%) (16.7%) (4.6%) (31.7%) (50.2%) (30.2%) (20.2%) (19.6%) (10.4%) (14.8%) (19.6%) 985.7 (33.2%) (49.1%) (32.3%) (27.4%) (19.9%) (13.5%) (–%) (19.9%) (15.0%)

(11.4%) (8.8%) (37.4%) (10.6%) (10.0%) (36.2%) (11.3%) (10.9%)

(26.7%) (9.4%) (44.7%) (18.7%) 1987.8 (46.6%) (21.9%) (–%)

(24.3%) (14.7%) (51.3%) (21.7%) (4.0%) (52.1%) (22.6%) (23.2%)

Solver performance on larger real benchmark datasets with 10 outliers. Methods marked with * are proposed in this

Method

alcohol (44, 6)

education (50, 4)

food (150, 3)

milk (86, 7)

*BnB (P) 111960 106436 88882 97331 *BnB (H) 1161 16069 616 17530 Gurobi (M) 80743 6666649 14265644 7271224 Gurobi (P) 84197 124485 63882 109642 *Gurobi (H) 4665 14080 578 15205 Mosek (M) 131507 4816500 854332 2019163 Mosek (P) 12767 105370 59726 90755 *Mosek (H) 1129695 – – 398930 Prop. 5 2471988 41752 22652 ∼ 109 Table 8

pulpfiber radar wagner (62, 7) (1573, 4) (63, 6) 99744 919 80587 5560 7103 105449 3022 582820 ∼ 108

105795 6253 1043403 684574 9715 80660 1191 – ∼ 109

100338 12136 481076 40401 13031 757731 31376 695613 ∼ 107

Number of BnB nodes explored by each solver on smaller real benchmark datasets. “–” indicates solver failure. The last row shows the theoretical upper bound from Proposition 5.

Method

pm10 pollen no fri rmftsa galaxy (500, 7) (3848, 4) (500, 7) (1000, 50) (508, 10) (323, 4)

*BnB (P) 92461 *BnB (H) 76871 Gurobi (M) 4748094 Gurobi (P) 1525648 *Gurobi (H) 522409 Mosek (M) 224224 Mosek (P) 9475 *Mosek (H) – Prop. 5 ∼ 1014 Table 9

88564 70127 22293 75180 19384 25079 47 47 ∼ 1010

96312 81740 5370154 1764024 170739 260771 20571 21270 ∼ 1014

91594 75705 365364 255223 135830 94819 4377 3913 ∼ 10104

83874 108458 75569 76574 4143545 4857546 1289427 2548249 276590 1674007 254029 372766 21400 51221 – 44744 ∼ 1019 ∼ 107

Number of BnB nodes explored by each solver on larger real benchmark datasets. “–” indicates solver failure. The last row shows the theoretical upper bound from Proposition 5.

Meng, Gómez and Mazumder: Computation of LTS: BnB and Hyperplane Arrangements Article submitted to Operations Research

% of instances solved

60 50 40 30 20

BnB (H) Gurobi (M) Gurobi (P) Gurobi (H)

10 0 10 1 Figure 5

100

101 Time (s)

102

103

Performance profile on the real datasets. For each time t, the curve reports the fraction of instances solved within time t (to the 1% optimality-gap tolerance).

explores more nodes than the strengthened formulation. Comparing Gurobi (H) to Gurobi (P), the hyperplane arrangement formulation explores only 30–60% as many nodes. In contrast, BnB (H) explores approximately 80% of the nodes that BnB (P) does. In other words, while the time required to solve the continuous relaxations increases with the hyperplane arrangement considerations, the proposed BnB method is comparatively less impacted, showcasing the benefit of using first order methods to solve the continuous relaxations. We observe that, in general, Gurobi (H) requires fewer branch-and-bound nodes than BnB (H) to prove optimality, although: (i) there are notable exceptions, see for example results with the “pulpfiber” dataset; (ii) the differences in number of nodes are typically small. In other words, we observe that the proposed BnB method does not seem to be substantially hampered from not solving the constrained problems to optimality. The last row of each table shows the theoretical upper bound on the number of nodes required for the strengthened formulation from Proposition 5. In practice, both BnB (H) and Gurobi (H) require only a small fraction of these nodes, showing that the methods perform much better in practice than what the worst-case theoretical bound would suggest. Overall, as summarized in Figure 5, BnB (H) attains the highest solve rate at short-to-moderate runtimes, while Gurobi (H) remains competitive and becomes advantageous on a subset of harder

Meng, Gómez and Mazumder: Computation of LTS: BnB and Hyperplane Arrangements Article submitted to Operations Research

instances. Both substantially outperform the remaining formulations, confirming that using hyperplane arrangement (H) is the dominant driver of performance on these benchmarks. 5.4. Value of Optimal Solutions We assess the practical benefit of computing optimal solutions compared to heuristic approaches. For each real dataset, we vary the number of outliers from 5 to 20 and compare the objective values obtained by the popularly used alternating minimization heuristic (Algorithm 3) against the solutions computed by the BnB solver. Figure 6 presents the relative suboptimality gap (fheuristic − foptimal )/foptimal as a box plot.

Heuristic v.s. MIP solution quality Gap 0 Figure 6

0.1%

10%

Box plot of relative suboptimality of heuristic solutions (Algorithm 3) compared to optimal solutions from BnB (H) on 13 real datasets, with outlier counts from 5 to 20.

The figure demonstrates the value of exact optimization. By construction, the optimal method always achieves an objective value at least as good as the heuristic (since the heuristic is called at the root node). In practice, the optimal method achieves a strictly better objective in more than 50% of cases. Furthermore, the gap exceeds 1% in more than 10% of cases, with the largest gap reaching approximately 30%. These findings indicate that while heuristics provide reasonable solutions in many cases, exact methods can yield substantially better solutions for a significant fraction of problem instances.

References Agulló J (2001a) New algorithms for computing the least trimmed squares regression estimator. Computational Statistics & Data Analysis 36(4):425–439. Agulló J (2001b) New algorithms for computing the least trimmed squares regression estimator. Computational statistics & data analysis 36(4):425–439. Aktürk MS, Atamtürk A, Gürel S (2009) A strong conic quadratic reformulation for machine-job assignment with controllable processing times. Operations Research Letters 37(3):187–191.

Meng, Gómez and Mazumder: Computation of LTS: BnB and Hyperplane Arrangements Article submitted to Operations Research

ApS M (2022) The MOSEK optimization toolbox for Python manual. Version 9.3. URL https://docs.mosek. com/latest/pythonfusion/index.html. Bernholt T (2006) Robust estimators are hard to compute. Technical Report Technical Report No. 2005,52, Universität Dortmund, SFB 475. Bertsekas DP (1982) Constrained Optimization and Lagrange Multiplier Methods (Academic Press). Bertsimas D, King A, Mazumder R (2016) Best subset selection via a modern optimization lens. The annals of statistics 44(2):813–852. Bertsimas D, Mazumder R (2014) Least quantile regression via modern optimization. The Annals of Statistics 42(6):2494–2525. Bhatia K, Jain P, Kar P (2015) Robust regression via hard thresholding. Advances in neural information processing systems 28. Ceria S, Soares J (1999) Convex programming for disjunctive convex optimization. Mathematical Programming 86(3):595–614. Černỳ M, Hladı́k M, Rada M (2019) Walks on hyperplane arrangements and optimization of piecewise linear functions. arXiv preprint arXiv:1912.12750 . Chang LC, Jones DK, Pierpaoli C (2012) Restore: robust estimation of tensors by outlier rejection. Magnetic Resonance in Medicine 68(2):538–552, URL http://dx.doi.org/10.1002/mrm.23268. Cover TM (1965) Geometrical and statistical properties of systems of linear inequalities with applications in pattern recognition. IEEE transactions on electronic computers (3):326–334. Edelsbrunner H, Souvaine DL (1990) Computing least median of squares regression lines and guided topological sweep. Journal of the American Statistical Association 85(409):115–119. Frangioni A, Gentile C (2006) Perspective cuts for a class of convex 0–1 mixed integer programs. Mathematical Programming 106:225–236. Giloni A, Padberg M (2002) Least trimmed squares regression, least median squares regression, and mathematical programming. Mathematical and Computer Modelling 35(9-10):1043–1060. Gómez A (2021) Outlier detection in time series via mixed-integer conic quadratic optimization. SIAM Journal on Optimization 31(3):1897–1925. Gómez A, Neto J (2025) Outlier detection in regression: conic quadratic formulations. INFORMS Journal on Computing . Günlük O, Linderoth J (2010) Perspective reformulations of mixed integer nonlinear programs with indicator variables. Mathematical programming 124:183–205. Gurobi Optimization, LLC (2022) Gurobi Optimizer Reference Manual. URL https://www.gurobi.com. Hawkins DM (1994) The feasible solution algorithm for least trimmed squares regression. Computational Statistics & Data Analysis 17(2):185–196.

Meng, Gómez and Mazumder: Computation of LTS: BnB and Hyperplane Arrangements Article submitted to Operations Research

Hawkins DM, Olive DJ (1999) Improved feasible solution algorithms for high breakdown estimation. Computational Statistics & Data Analysis 30(1):1–11. Hazimeh H, Mazumder R (2020) Fast best subset selection: Coordinate descent and local combinatorial optimization algorithms. Operations Research 68(5):1517–1537. Hazimeh H, Mazumder R, Saab A (2021) Sparse regression at scale: Branch-and-bound rooted in first-order optimization. Mathematical Programming 1–42. Hofmann M, Gatu C, Kontoghiorghes EJ (2010) An exact least trimmed squares algorithm for a range of coverage values. Journal of Computational and Graphical Statistics 19(1):191–204. Holland PW, Welsch RE (1977) Robust regression using iteratively reweighted least-squares. Communications in Statistics - Theory and Methods 6(9):813–827, URL http://dx.doi.org/10.1080/ 03610927708827533. Hössjer O (1995) Exact computation of the least trimmed squares estimate in simple linear regression. Computational Statistics & Data Analysis 19(3):265–282. Huber PJ (1964) Robust estimation of a location parameter. The Annals of Mathematical Statistics 35(1):73–101, URL http://dx.doi.org/10.1214/aoms/1177703732. Huber PJ (1981) Robust statistics. Wiley Series in Probability and Mathematical Statistics . Huchette J, Vielma JP (2019) A geometric way to build strong mixed-integer programming formulations. Operations Research Letters 47(6):601–606. Insolia L, Kenney A, Chiaromonte F, Felici G (2022) Simultaneous feature selection and outlier detection with optimality guarantees. Biometrics 78(4):1592–1603. Klouda K (2015) An exact polynomial time algorithm for computing the least trimmed squares estimate. Computational Statistics & Data Analysis 84:27–40. Lam SK, Pitrou A, Seibert S (2015) Numba: A llvm-based python jit compiler. Proceedings of the Second Workshop on the LLVM Compiler Infrastructure in HPC, 1–6. Linderoth JT, Savelsbergh MW (1999) A computational study of search strategies for mixed integer programming. INFORMS Journal on Computing 11(2):173–187. Maronna RA, Martin RD, Yohai VJ (2006) Robust Statistics: Theory and Methods (Chichester: John Wiley & Sons). Motulsky HJ, Brown RE (2006) Detecting outliers when fitting data with nonlinear regression: a new method based on robust nonlinear regression and the false discovery rate. BMC Bioinformatics 7:123, URL http://dx.doi. org/10.1186/1471-2105-7-123. Rockafellar RT (1976) Augmented lagrangians and applications of the proximal point algorithm in convex programming. Mathematics of Operations Research 1(2):97–116. Rousseeuw PJ (1984) Least median of squares regression. Journal of the American Statistical Association 79(388):871–880.

Meng, Gómez and Mazumder: Computation of LTS: BnB and Hyperplane Arrangements Article submitted to Operations Research

Rousseeuw PJ, Leroy AM (1987) Robust Regression and Outlier Detection (New York: John Wiley & Sons). Rousseeuw PJ, Leroy AM (2003) Robust regression and outlier detection (John wiley & sons). Rousseeuw PJ, Van Driessen K (2006a) Computing lts regression for large data sets. Data Mining and Knowledge Discovery 12(1):29–45. Rousseeuw PJ, Van Driessen K (2006b) Computing lts regression for large data sets. Data mining and knowledge discovery 12(1):29–45. Shen Y, Sanghavi S (2019a) Iterative least trimmed squares for mixed linear regression. Advances in Neural Information Processing Systems 32. Shen Y, Sanghavi S (2019b) Learning with bad training data via iterative trimmed loss minimization. International conference on machine learning, 5739–5748 (PMLR). Sun Q, Mao R, Zhou WX (2021) Adaptive capped least squares. arXiv preprint URL http://dx.doi.org/10. 48550/arXiv.2107.00109. Vanschoren J, Van Rijn JN, Bischl B, Torgo L (2014) Openml: networked science in machine learning. ACM SIGKDD Explorations Newsletter 15(2):49–60. Wolsey LA (2020) Branch and bound. Integer Programming, chapter 7 (John Wiley & Sons), 2 edition. Zaslavsky T (1975) Facing up to arrangements: Face-count formulas for partitions of space by hyperplanes, volume 154 (American Mathematical Soc.). Zioutas G, Avramidis A (2005a) Deleting outliers in robust regression with mixed integer programming. Acta Mathematicae Applicatae Sinica, English Series 21(2):323–334. Zioutas G, Avramidis A (2005b) Deleting outliers in robust regression with mixed integer programming. Acta Mathematicae Applicatae Sinica 21(2):323–334. Zioutas G, Pitsoulis L, Avramidis A (2009) Quadratic mixed integer programming and support vectors for deleting outliers in robust regression. Annals of Operations Research 166(1):339–353.

Meng, Gómez and Mazumder: Computation of LTS: BnB and Hyperplane Arrangements Article submitted to Operations Research

Appendix A: Additional Technical Details A.1. Proof of Proposition 3 In this section we derive the convex hull of set ZHA (b), which we repeat for convenience.

ZHA = (z, w, r, t) ∈ {0, 1} × R3 : t ≥ w2 , w(1 − z) = 0, |r|(1 − z) ≤ b(1 − z), |r|z ≥ bz, w = rz ;

we omit the explicit dependence on parameter b ∈ R+ for convenience, as it is fixed for the purposes of this section. To prove the result, we use disjunctive programming (Ceria and Soares 1999). Indeed, set ZHA = Z− ∪ Z+ ∪ Z= where

def Z− = (z, w, r, t) ∈ R4 : t ≥ w2 , z = 1, r ≤ −b, w = r def Z+ = (z, w, r, t) ∈ R4 : t ≥ w2 , z = 1, r ≥ b, w = r def Z= = (z, w, r, t) ∈ R4 : t ≥ w2 , z = 0, |r| ≤ b, w = 0 .

Using standard disjunctive programming reformulations, we find that (z̄, w̄, r̄, t̄) ∈ cl conv(ZHA ) if and only if there exists additional variables (z− , w− , r− , t− , α− ), (z+ , w+ , r+ , t+ , α+ ) and (z= , w= , r= , t= , α= ) such that the system

α− + α+ + α= = 1, α− ≥ 0, α+ ≥ 0, α= ≥ 0 z̄ = z− + z+ + z= , w̄ = w− + w+ + w= , r̄ = r− + r+ + r= , t̄ = t− + t+ + t= 2 t− α− ≥ w− , z− = α− , r− ≤ −bα− , w− = r− 2 t+ α+ ≥ w+ , z+ = α+ , r+ ≥ bα+ , w+ = r+ 2 t= α= ≥ w= , z= = 0, |r= | ≤ bα= , w= = 0.

Meng, Gómez and Mazumder: Computation of LTS: BnB and Hyperplane Arrangements Article submitted to Operations Research

To obtain the result, we project out the additional variables. First, we remove auxiliary variables (z, w) by using the equality constraints, resulting in the simplified system

α− + α+ + α= = 1, α− ≥ 0, α+ ≥ 0, α= ≥ 0 z̄ = α− + α+ , w̄ = r− + r+ , r̄ = r− + r+ + r= , t̄ = t− + t+ + t= 2 t− α− ≥ r− , r− ≤ −bα− 2 , r+ ≥ bα+ t+ α+ ≥ r+

t= ≥ 0, |r= | ≤ bα= .

2 2 Next we can project out auxiliary variables t, leading inequality t̄ ≥ r− /α− + r+ /α+ , and auxiliary

variable α= using the first equality constraint. The simplified system reads

α− ≥ 0, α+ ≥ 0, α− + α+ ≤ 1 z̄ = α− + α+ , w̄ = r− + r+ , r̄ = r− + r+ + r= 2 2 t̄ ≥ r− /α− + r+ /α+

r− ≤ −bα− , r+ ≥ bα+ , |r= | ≤ b(1 − α− − α+ ).

Finally, we project out r= since r̄ = r− + r+ + r= ⇔ r̄ = w̄ + r= , and thus r= satisfying constraints exists if and only |r̄ − w̄| ≤ b(1 − α− − α+ ). The result of Proposition 3 then follows by renaming variables as α− ↔ z− , α+ ↔ z+ , r− ↔ −w− and r+ ↔ w+ . A.2. Proof of Proposition 7 and Proposition 8 We first establish a key lemma about the structure of optimal solutions of the relaxation problem. L EMMA 1. Any optimal solution of (11) with any given ℓ− , ℓ+ , u− , u+ satisfies wi+ wi− = 0 and zi+ zi− = 0 for all i ∈ [n]. Proof. Given an index i, we discuss the following three cases. + + − Case 1: u− i = ui = 0. This implies zi = zi = 0. The perspective terms

di (wi+ )2 di (wi− )2 and then + zi zi−

force wi+ = wi− = 0. Thus wi+ wi− = zi+ zi− = 0. + + − Case 2: ℓ− i = 1 or ℓi = 1. This implies zi + zi = 1. Constraint (9b) becomes |ri − wi | ≤ 0, which

forces wi = ri . For a fixed residual ri , we determine the optimal values of the auxiliary variables.

Meng, Gómez and Mazumder: Computation of LTS: BnB and Hyperplane Arrangements Article submitted to Operations Research

If ℓ+ i = 1 (i.e., wi ≥ 0), then ri = wi ≥ 0. If the problem is feasible, the objective is minimized by setting zi+ = 1, zi− = 0, wi+ = ri , and wi− = 0, which satisfies wi = wi+ − wi− = ri and minimizes the perspective term. If ℓ− i = 1 (i.e., wi ≤ 0), then ri = wi ≤ 0. If the problem is feasible, the objective is minimized by setting zi+ = 0, zi− = 1, wi+ = 0, and wi− = −ri , which satisfies wi = wi+ − wi− = ri and minimizes the perspective term. In both cases, wi+ wi− = zi+ zi− = 0. + − + Case 3: ℓ− i = ℓi = 0, ui = ui = 1. We prove by contradiction. Suppose at an optimal solution, √ both zi+ > 0 and zi− > 0 for some i. By constraint (9a), we have wi+ ≥ 2µzi+ > 0 and wi− ≥ √ 2µzi− > 0.

Consider a perturbation: replace (zi+ , zi− , wi+ , wi− ) by (αzi+ , βzi− , αwi+ , βwi− ) for some 0 < α, β < 1 chosen such that αwi+ − βwi− = wi+ − wi− . This equation can be rewritten as wi+ (1 − α) = wi− (1 − β). Since wi+ , wi− > 0, we can choose α, β ∈ (0, 1) satisfying this constraint. We verify that all constraints remain satisfied: √ √ • Constraint (9a): αwi+ ≥ 2µαzi+ and βwi− ≥ 2µβzi− hold since the original constraints hold. • Constraint (9b): The left-hand side |ri − wi | = |ri − (wi+ − wi− )| is unchanged since αwi+ − √ βwi− = wi+ − wi− . The right-hand side 2µ(1 − zi ) increases since the new zi = αzi+ + βzi− < zi+ + zi− . Thus the constraint remains satisfied. • The new zi = αzi+ + βzi− < zi+ + zi− , so the constraint zi ∈ [0, 1] remains satisfied. The objective changes only in the perspective terms. The new contribution is di (αwi+ )2 di (βwi− )2 di (wi+ )2 di (wi− )2 di (wi+ )2 di (wi− )2 + =α +β < + , αzi+ βzi− zi+ zi− zi+ zi− since α, β < 1. This contradicts the optimality of the original solution. Therefore, at any optimal solution, zi+ zi− = 0. By constraint (9a), this implies wi+ wi− = 0.

□

Lemma 1 allows us to project out the auxiliary variables zi+ , zi− , wi+ , wi− by writing di (wi+ )2 di (wi− )2 di wi2 + = , zi zi+ zi−

Meng, Gómez and Mazumder: Computation of LTS: BnB and Hyperplane Arrangements Article submitted to Operations Research

− + where zi = zi+ + zi− and wi = wi+ − wi− . Note that ℓ− i = 1 implies zi = 1, zi = 0 and therefore wi+ = 0, w ≤ 0. Similarly, ℓ+ i = 1 implies w ≥ 0. We can rewrite (13) as

2 1 w 2 2 ϕ(r; ℓ , ℓ , u , u , d) = min (w − r) + µz + d −w w∈R,z∈[0,1] 2 z p p s.t. |w| ≥ 2µz, |r − w| ≤ 2µ(1 − z) −

−

z ≤ u+ + u− , w ≥ 0 if ℓ+ = 1, w ≤ 0 if ℓ− = 1.

(22a) (22b) (22c)

We now discuss the following four Cases, which gives us the closed form expression of ϕ in Proposition 7 and Proposition 8. Case 1: u− = u+ = 0. We have z = 0. The constraint |r − w| ≤ √ (forced by the perspective term) yields |r| ≤ 2µ, which gives

√ 2µ combined with w = 0

1 ϕ(r; 0, 0, 0, 0, d) = r2 + δ[−√2µ,√2µ] (r). 2 Case 2: l+ = 1. we have z = 1 and w ≥ 0. The constraint |r − w| ≤ 0 forces w = r. Combined √ √ with |w| ≥ 2µ, we obtain r ≥ 2µ, , which gives ϕ(r; 0, 1, 0, 1, d) = µ + δ[√2µ,∞] (r). √ √ Case 3: l− = 1. We have z = 1 and w ≤ 0. Similarly, w = r and w ≤ − 2µ yield r ≤ − 2µ, which gives ϕ(r; 1, 0, 1, 0, d) = µ + δ[−∞,−√2µ] (r). Case 4: ℓ− = ℓ+ = 0, u− = u+ = 1. We first ignore the constraints in (22b) and derive the closed-form solution for the optimization problem 2 1 w 2 2 ϕ(r) = min (w − r) + µz + d −w , w∈R,z∈[0,1] 2 z

(23)

where µ > 0 and 0 < d < 1/2. Step 1: Minimization over w for fixed z ∈ (0, 1]. For any z > 0, the objective is strictly convex in w. The coefficient of w2 is 1 1 1 d a(z) := + d − 1 = − d + > 0. 2 z 2 z The unique minimizer is w∗ (z) =

r 1 + 2d

1 −1 z

r . 1 − 2d + 2d z

(24)

Meng, Gómez and Mazumder: Computation of LTS: BnB and Hyperplane Arrangements Article submitted to Operations Research

Substituting back and completing the square, the minimal value over w is 1 r2 . g(z) := min ϕ(w; z) = µz + r2 − w 2 2 1 − 2d + 2d z

(25)

When z = 0, feasibility forces w = 0, yielding g(0) = 12 r2 . At z = 1, we have w∗ (1) = r and g(1) = µ. The original problem thus reduces to ϕ(r) = min g(z). z∈[0,1]

Step 2: Stationarity in z and the interior solution. Define D(z) := 1 − 2d + 2d = (1−2d)z+2d for z z z ∈ (0, 1]. Differentiating (25) gives g ′ (z) = µ −

r2 d . z 2 D(z)2

Setting g ′ (z) = 0 and noting that z 2 D(z)2 = N (z)2 where N (z) := (1 − 2d)z + 2d, we obtain s d N (z) = |r| . µ Solving for z yields the candidate interior minimizer q |r| µd − 2d ∗ z (r) = . 1 − 2d p √ This is feasible (i.e., 0 ≤ z ∗ ≤ 1) precisely when 2 µd ≤ |r| ≤ µ/d.

(26)

p √ Step 3: Evaluation at the optimizer and boundary regimes. When 2 µd ≤ |r| ≤ µ/d, we p have z ∗ (r) ∈ (0, 1) and g ′ (z ∗ ) = 0. Using N (z ∗ ) = |r| d/µ and D(z ∗ ) = N (z ∗ )/z ∗ in (25): 1 r2 z ∗ ϕ(r) = g(z ∗ ) = µz ∗ + r2 − . 2 2N (z ∗ ) Substituting z ∗ from (26) and simplifying yields √ −dr2 + 2 µd|r| − 2µd ϕ(r) = . 1 − 2d The optimal w∗ in this regime is obtained from (24): rz ∗ r w (r) = = = sign(r) D(z ∗ ) N (z ∗ ) ∗

µ ∗ z (r). d

√ When |r| < 2 µd, the stationary z ∗ in (26) is negative and infeasible. Since g is convex on [0, 1], the minimum occurs at z = 0, which forces w = 0. Thus ϕ(r) = 12 r2 .

Meng, Gómez and Mazumder: Computation of LTS: BnB and Hyperplane Arrangements Article submitted to Operations Research

When |r| >

p µ/d, we have z ∗ > 1, which is infeasible. The minimum occurs at z = 1, where

w∗ (1) = r and ϕ(r) = µ. Summary. Combining the three cases:  √ 1 2   r , |r| ≤ 2 µd,   2    r  −dr2 + 2√µd|r| − 2µd √ µ ϕ(r) = , 2 µd ≤ |r| ≤ ,  1 − 2d d   r   µ   µ, |r| ≥ , d with optimal (w∗ , z ∗ ) given by

∗

z (r) =

   0,    q     |r| d − 2d µ  1 − 2d        1,

√ |r| ≤ 2 µd, ,

√ 2 µd ≤ |r| ≤ r |r| ≥

   0,    r   µ ∗ ∗ z (r), w (r) = sign(r) d       r, Derivative of ϕ(r). Differentiating each piece yields:    r,       −2dr + 2√µd sign(r) , ϕ′ (r) = 1 − 2d        0,

µ , d

√ |r| ≤ 2 µd, √

2 µd ≤ |r| ≤ r µ |r| ≥ . d

µ , d

√ |r| < 2 µd, √ 2 µd < |r| < r µ |r| > . d

µ , d

√ At the transition points, continuity can be verified directly. At |r| = 2 µd: the left derivative is √ √ √ √ −2d·2 µd sign(r)+2 µd sign(r) sign(r) · 2 µd, and the right derivative is = sign(r) · 2 µd. At |r| = 1−2d √ √ p −2d µ/d sign(r)+2 µd sign(r) µ/d: the left derivative is = 0. Thus ϕ′ (r) is continuous everywhere. 1−2d Checking feasibility in (22b). √ √ • When |r| ≤ 2 µd: We have z ∗ = w∗ = 0. Both constraints reduce to 0 ≥ 0 and |r| ≤ 2µ. √ √ The latter holds since |r| ≤ 2 µd ≤ 2µ (using d ≤ 1/2).

Meng, Gómez and Mazumder: Computation of LTS: BnB and Hyperplane Arrangements Article submitted to Operations Research

p p √ • When 2 µd ≤ |r| ≤ µ/d: We have w∗ = sign(r) µ/dz ∗ . The first constraint becomes p √ µ/dz ∗ ≥ 2µz ∗ , which holds since d ≤ 1/2. For the second constraint, we compute |r − p w∗ | = |r| − |w∗ | = |r| − µ/dz ∗ . Substituting z ∗ and simplifying shows that this equals √ 2µ(1 − z ∗ ). p √ • When |r| ≥ µ/d: We have z ∗ = 1 and w∗ = r. The constraints become |r| ≥ 2µ and |r − r| = 0 ≤ 0, both of which hold. Since the optimal (w∗ , z ∗ ) satisfies constraints (22b), they are also the optimal solutions to ϕ(r; 0, 0, 1, 1, d). Moreover, the constraints (13c) can be removed. A.3. Proof of Proposition 9 We first establish a preliminary lemma. L EMMA 2 (Gradient-suboptimality inequality). If Φ : Rp → R is λ̃-strongly convex and differentiable, then for any x ∈ Rp , min Φ(u) ≥ Φ(x) −

u∈Rp

1 ∥∇Φ(x)∥22 . 2λ̃

Proof. Strong convexity yields for all u: Φ(u) ≥ Φ(x) + ⟨∇Φ(x), u − x⟩ + λ̃2 ∥u − x∥22 . The right side is minimized in u at u = x − λ̃−1 ∇Φ(x), giving the result.

□

As mentioned in Section 2.1.2, the perspective relaxation (12) is designed to be λ̃-strongly convex. Since the augmented Lagrangian (16) is obtained by removing some constraints in (12) and adding some piecewise-linear/quadratic functions, it is also λ̃-strongly convex. Therefore, ζ̄(ν̄) = min Lρ (β; ν̄) ≥ Lρ (β̄, ν̄) − β

∥∇β Lρ (β̄, ν̄)∥22 2λ̃

where the right-hand side provides a lower bound of the relaxation problem.

□

A.4. Computing Incumbent Solutions We obtain incumbent solutions in the BnB procedure by solving an ℓ2 -regularized estimation problem with the outlier set S ⊂ [n] fixed. Specifically, for a given support S, we solve the original problem (1) with zi = 1 for i ∈ S and zi = 0 otherwise: F ∗ (S) := minp β∈R

1X λ 2 2 (yi − x⊤ i β) + ∥β∥2 + µ|S|. 2 2

(27)

i∈S /

This is a weighted ridge regression problem with a closed-form solution. Let W = Diag(1 − z) be the diagonal weight matrix. The optimal coefficient vector is β ∗ (S) = (X ⊤ W X + λI)−1 X ⊤ W y.

(28)

Meng, Gómez and Mazumder: Computation of LTS: BnB and Hyperplane Arrangements Article submitted to Operations Research

Before constructing the BnB tree, following (Rousseeuw and Van Driessen 2006b), we compute an initial incumbent solution using a heuristic based on alternating minimization. The method alternates between optimizing over β for fixed z and optimizing over z for fixed β. Starting with z = 0 (i.e., no outliers), we iterate as follows: 1. Solve for β: Given the current z, compute β using (28). 2. Update z: Given the current β, compute residuals ri = yi − x⊤ i β and update each zi by comparing the squared loss to the outlier penalty:   1, if 1 r2 > µ, 2 i zi =  0, otherwise.

(29)

The update rule (29) follows from the optimality condition: for fixed β, observation i should be an outlier if and only if the penalty µ is less than the squared loss 12 ri2 . This alternating procedure continues until the outlier set S = {i : zi = 1} stabilizes, which is guaranteed since the objective decreases monotonically. Algorithm 3 summarizes the method. Algorithm 3 Computing initial incumbent via alternating minimization Input: Data (X, y), parameters λ, µ 1: Initialize z = 0 2: repeat 3:

Compute β = (X ⊤ W X + λI)−1 X ⊤ W y where W = Diag(1 − z)

Compute residuals ri = yi − x⊤ i β for all i ∈ [n]

Update zi = 1{ 21 ri2 > µ} for all i ∈ [n]

6: until the outlier set S = {i : zi = 1} does not change 7: return β, outlier set S

At each node in the BnB tree, we further refine the upper bound by solving problem (27) with the support set determined by rounding the relaxed integer variables: S = {i : zi− + zi+ ≥ 0.5}. Note that variables zi− + zi+ are only implicitly computed, and this is equivalent to choosing ( ) 3 r µ 2 µ − d − d 4 1 S = i : |yi − x⊤ or S = i : ϕ(yi − x⊤ . i β| ≥ 2 + d i β) ≥ d 1 − 2d

Record · ID 10350 · SHA-256 8edb2701b0802a94
Conceptio Open Knowledge Archive — every document is proof-bundled with source, license, and retrieval metadata.