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

Xiang Meng; Andrés Gómez; Rahul Mazumder

math.OCApr 13, 2026

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

Authors: Xiang Meng, Andrés Gómez, Rahul Mazumder

Organizations: MIT Operations Research Center · USC Daniel J. Epstein Department of Industrial and Systems Engineering · MIT Sloan School of Management, Operations Research Center

Abstract

We study computational aspects of a key problem in robust statistics -- the penalized least trimmed squares (LTS) regression problem, a robust estimator that mitigates the influence of outliers in data by capping residuals with large magnitudes. Although statistically attractive, penalized LTS is NP-hard, and existing mixed-integer optimization (MIO) formulations scale poorly due to weak relaxations and exponential worst-case complexity in the number of observations. We propose a new MIO formulation that embeds hyperplane arrangement logic into a perspective reformulation, explicitly enforcing structural properties of optimal solutions. We show that, if the number of features is fixed, the resulting branch-and-bound tree is of polynomial size in the sample size. Moreover, we develop a tailored branch-and-bound algorithm that uses first-order methods with dual bounds to solve node relaxations efficiently. Computational experiments on synthetic and real datasets demonstrate substantial improvements over existing MIO approaches: on synthetic instances with 5000 samples and 20 features, our tailored solver reaches a 1% gap in 1 minute while competing approaches fail to do so within one hour. These gains enable exact robust regression at significantly larger sample sizes in low-dimensional settings.

Figures & tables

Figure 1: Various robust loss functions.

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

Figure 3: A depth 5 integer nodes with four relatives.

Figure 4: Comparison of loss functions for a single residual

r

: the original least trimmed loss

\phi_{\text{cap}}

(before relaxation), the perspective relaxation

\bar{\phi}(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.

Disjunction	Replace $\bar{\phi}(y_{i}-\bm{x_{i}^{\top}\beta};d_{i})$ with…
$z_{i}^{-}=0,z_{i}^{+}=0$	$\bar{\phi}_{0}(\cdot;\nu_{1},\nu_{2},\rho)$
$z_{i}^{-}=1,z_{i}^{+}=0$	$\bar{\phi}_{-}(\cdot;\nu,\rho)$
$z_{i}^{-}=0,z_{i}^{+}=1$	$\bar{\phi}_{+}(\cdot;\nu,\rho)$

Table 1: Subproblem construction

Method	$n=1000$			$n=2000$			$n=5000$
Method	$p\!=\!10$	$p\!=\!20$	$p\!=\!50$	$p\!=\!10$	$p\!=\!20$	$p\!=\!50$	$p\!=\!10$	$p\!=\!20$	$p\!=\!50$
*BnB (P)	(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)
BnB (H)	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)
Gurobi (H)	358.5	414.0	766.4	(1.5%, 8)	2301.2	(6.6%, 5)	(14.0%, 0)	(16.1%, 0)	(19.6%, 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.

Method	$n=50$			$n=100$			$n=200$
Method	$p\!=\!5$	$p\!=\!10$	$p\!=\!20$	$p\!=\!5$	$p\!=\!10$	$p\!=\!20$	$p\!=\!5$	$p\!=\!10$	$p\!=\!20$
Gurobi (P)	2.5	2.2	4.2	21.9	26.3	19.2	447.0	87.5	97.0
Gurobi (H)	2.0	2.8	7.1	7.5	23.3	33.3	27.2	39.8	77.9
BnB (H)	2.0	3.9	34.5	10.5	12.5	11.9	9.5	9.0	50.0

Table 3: 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.

$\lambda_{0}$	Gurobi (P)	*Gurobi (H)	*BnB (H)
0.0005	(66.3%, 0)	(61.7%, 0)	(32.4%, 0)
0.001	(57.9%, 0)	(53.6%, 0)	(19.1%, 0)
0.002	(46.4%, 0)	(43.4%, 0)	(5.2%, 4)
0.005	(27.8%, 0)	(18.2%, 1)	34.7
0.01	(14.3%, 0)	(1.2%, 8)	10.3
0.02	(4.9%, 0)	1627.2	4.3

Table 4: 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.

SNR	Gurobi (P)	*Gurobi (H)	*BnB (H)
5.0	(52.7%, 0)	(52.2%, 0)	(45.0%, 0)
10.0	(45.7%, 0)	(44.8%, 0)	(34.8%, 0)
20.0	(32.7%, 0)	(33.6%, 0)	(13.1%, 1)
30.0	(24.1%, 0)	(18.5%, 1)	(1.3%, 8)
40.0	(18.4%, 0)	(1.0%, 9)	18.9
50.0	(14.3%, 0)	(1.2%, 8)	9.4

Table 5: 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.

Method	alcohol	education	food	milk	pulpfiber	radar	wagner
Method	$(44,6)$	$(50,4)$	$(150,3)$	$(86,7)$	$(62,7)$	$(1573,4)$	$(63,6)$
*BnB (P)	(28.6%)	(14.6%)	(7.6%)	(10.9%)	(29.0%)	(6.7%)	(14.7%)
BnB (H)	3.4	118.3	1.0	128.9	2.0	23.7	56.2
Gurobi (M)	11.9	754.6	(18.5%)	1564.5	16.6	(23.5%)	80.0
Gurobi (P)	25.4	38.5	39.2	51.8	7.6	(4.4%)	19.3
Gurobi (H)	5.4	13.6	3.0	28.8	10.6	660.2	12.9

Table 6: 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.

Method	pm10	pollen	no	fri	rmftsa	galaxy
Method	$(500,7)$	$(3848,4)$	$(500,7)$	$(1000,50)$	$(508,10)$	$(323,4)$
*BnB (P)	(23.8%)	(19.2%)	(15.0%)	(11.4%)	(26.7%)	(24.3%)
BnB (H)	(15.6%)	(16.7%)	(4.6%)	(8.8%)	(9.4%)	(14.7%)
Gurobi (M)	(31.7%)	(50.2%)	(30.2%)	(37.4%)	(44.7%)	(51.3%)
Gurobi (P)	(20.2%)	(19.6%)	(10.4%)	(10.6%)	(18.7%)	(21.7%)
Gurobi (H)	(14.8%)	(19.6%)	985.7	(10.0%)	1987.8	(4.0%)

Table 7: Solver performance on larger real benchmark datasets with 10 outliers. Methods marked with * are proposed in this work. Dataset dimensions

(n,p)

Method	alcohol	education	food	milk	pulpfiber	radar	wagner
Method	$(44,6)$	$(50,4)$	$(150,3)$	$(86,7)$	$(62,7)$	$(1573,4)$	$(63,6)$
*BnB (P)	111960	106436	88882	97331	99744	105795	100338
BnB (H)	1161	16069	616	17530	919	6253	12136
Gurobi (M)	80743	6666649	14265644	7271224	80587	1043403	481076
Gurobi (P)	84197	124485	63882	109642	5560	684574	40401
Gurobi (H)	4665	14080	578	15205	7103	9715	13031

Table 8: 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 3.5 .

Method	pm10	pollen	no	fri	rmftsa	galaxy
Method	$(500,7)$	$(3848,4)$	$(500,7)$	$(1000,50)$	$(508,10)$	$(323,4)$
*BnB (P)	92461	88564	96312	91594	83874	108458
BnB (H)	76871	70127	81740	75705	75569	76574
Gurobi (M)	4748094	22293	5370154	365364	4143545	4857546
Gurobi (P)	1525648	75180	1764024	255223	1289427	2548249
Gurobi (H)	522409	19384	170739	135830	276590	1674007

Table 9: 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 3.5 .

Figure 5: 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).

Figure 6: 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.

Robust Regression

Robust regression develops methods that predict reliably despite outliers, heavy-tailed noise, corrupted data, or shifts in the data distribution.

See all 28 papers

Robustly Reliable Learning for Regression
Oct 12, 2026Avrim Blum, Anh Tuan NguyenRobust Regression
TTIC · Carnegie Mellon University
Simpler Methods Work Better for L1 Penalized Logistic Models and Large Datasets
Sep 21, 2026Edward Raff, James HoltLogistic Regression
CrowdStrike, Inc., USA · Univ. of Maryland, Baltimore County, USA
Global Truncated Loss Minimization for Robust and Threshold-Resilient Geometric Estimation
Mar 16, 2026Tianyu Huang, Liangzu Peng, Xinyue Zhang +5Robust Loss Functions Branch and Bound
Hong Kong Centre For Logistics Robotics, The Chinese University of Hong Kong. · Center for Innovation in Data Engineering and Science (IDEAS), University of Pennsylvania. · Mobile Perception Lab of the School of Information and Technology, ShanghaiTech University. +1

Oct 12, 2026J/K move · Enter open · S save

Robustly Reliable Learning for Regression

Avrim Blum, Anh Tuan Nguyen

TTIC · Carnegie Mellon University

Robust Regression

Open paper AskPDF arXiv

We introduce the problem of robustly reliable learning for the regression setting in adversarial environments. Given a test point

x_{test}

, a training dataset

S

n

data points, a corruption budget

b \ll n

, and a real-valued hypothesis class

\mathcal H

, our goal is to output the optimal prediction range

I^*_{x_{test}}(S, b)

, such that any

L

-Lipschitz predictor consistent with all but at most

b

points in

S

must predict a value

h(x_{test}) \in I

. We then instantiate the analysis for

L

-Lipschitz hypothesis classes, a natural and widely studied class in the regression setting, focusing on two key scenarios. First, for the restricted case where the corrupted dataset is guaranteed to be

L

-Lipschitz consistent, we propose a simple and efficient

b^{th}

-best heuristic algorithm and prove its correctness in computing the tightest possible range

I^*_{x_{test}}(S, b)

. We then establish its computational complexity of

O(n)

and show that it is computationally optimal in this scenario. Building on this foundation, we address the general scenario where the dataset

S

is no longer

L

-Lipschitz consistent. We demonstrate that computing the optimal range in this setting reduces to analyzing a violation graph, which models the pairwise

L

-Lipschitz violations among data points. By leveraging the structural properties of this graph's kernel, we propose a Fixed-Parameter Tractable (FPT) algorithm that computes the exact optimal range in

O(n^2 + b^4)

time. Furthermore, under the Strong Exponential Time Hypothesis (SETH), we establish a computational lower bound of

Ω(n^{2 - o(1)})

, proving that our algorithm is computationally optimal. Finally, we discuss open directions for the problem of robustly reliable learning for regression.

Abstract

Robust Regression

Explore similar work

Robustly Reliable Learning for Regression

Simpler Methods Work Better for L1 Penalized Logistic Models and Large Datasets

Global Truncated Loss Minimization for Robust and Threshold-Resilient Geometric Estimation