Fast and Exact Least Absolute Deviations Line Fitting via Piecewise Affine Lower-Bounding

TL;DR

The paper tackles robust LAD line fitting by introducing PALB, an exact solver that constructs piecewise-affine lower bounds from subgradients of the marginal LAD objective $J(m)=\min_t\sum_i|y_i-mx_i-t|$ and refines the slope interval until the subdifferential contains zero. PALB achieves finite termination with provable complexity bounds linked to the $k$-level/$k$-set problems, and exhibits log-linear practical performance, substantially faster than public LP-based and IRLS solvers on synthetic and real data. The implementation in Rust with a Python API, along with a comprehensive benchmark against L1pack, CPLEX, Clarabel, HiGHS, and Statsmodels, demonstrates PALB’s speed, exactness, and low memory footprint, making it suitable for embedded and real-time applications. The work also provides rigorous proofs for subdivision and expansion steps and discusses performance profiles and iteration counts, reinforcing PALB as a scalable, exact solution for LAD regression.

Abstract

Least-absolute-deviations (LAD) line fitting is robust to outliers but computationally more involved than least squares regression. Although the literature includes linear and near-linear time algorithms for the LAD line fitting problem, these methods are difficult to implement and, to our knowledge, lack maintained public implementations. As a result, practitioners often resort to linear programming (LP) based methods such as the simplex-based Barrodale-Roberts method and interior-point methods, or on iteratively reweighted least squares (IRLS) approximation which does not guarantee exact solutions. To close this gap, we propose the Piecewise Affine Lower-Bounding (PALB) method, an exact algorithm for LAD line fitting. PALB uses supporting lines derived from subgradients to build piecewise-affine lower bounds, and employs a subdivision scheme involving minima of these lower bounds. We prove correctness and provide bounds on the number of iterations. On synthetic datasets with varied signal types and noise including heavy-tailed outliers as well as a real dataset from the NOAA's Integrated Surface Database, PALB exhibits empirical log-linear scaling. It is consistently faster than publicly available implementations of LP based and IRLS based solvers. We provide a reference implementation written in Rust with a Python API.

Figures (5)

• Figure 1: An example showing how the interval boundaries $a_k, b_k$ evolve from one step to the next. Starting from an interval $[a_0,b_0]$, the approximation $J^{(1)}$ is formed by the two supporting lines at the endpoints: lines 2 and 5. The minimizer of $J^{(1)}$, found at the intersection of these two lines, becomes the new endpoint $a^{(1)}$. This yields a new interval $[a_1, b_1]$ with $b_1 = b_0$. Similarly intervals $[a_2, b_0]$ and $[a_2,b_3]$ are computed. At this point the approximation $J^{(3)}$ has the same minimizer $m^*$ as the objective and the algorithm terminates with that minimizer. Note that the lines 4 and 6 were "skipped" along the way. Finally note that the lines in this example do not necessarily correspond to those of an actual instance of the $L_1$ regression problem; they are instead manually chosen to qualitatively visualize the general algorithm.
• Figure 2: Median runtime results for the experiments on synthetic data.
• Figure 3: The performance profiles for the synthetic data.
• Figure 4: The performance profiles for the real data.
• Figure 5: The number of steps our method took to produce a solution on synthetic and real data. (a) The indicated regions show the distribution of the steps across all realizations. Each region is bounded by the indicated $p$-quantiles. (b) Each point indicates the number of steps taken for one timeseries. (The synthetic and real data are visualized differently because plotting distinct datapoints for every sample in the synthetic data looks rather cluttered, while there aren't enough samples per input size to produce a meaningful quantile plot in the case of real data. The dashed and dotted lines show $L_{1}$ and $L_{2}$ lines fitted to the measurements. The $L_{2}$ fit is accomplished using numpy's np.linalg.lstsq method, while the $L_{1}$ fit is due to our method. There are some points in the samples that can be considered outliers --- these can be seen to throw the $L_{2}$ lines off from what one might intuitively expect.)

Theorems & Definitions (10)

• Theorem 2.1
• Theorem 2.2
• Lemma 5.1
• Lemma 5.2
• Theorem 5.3
• Proposition 5.4
• Proposition 5.5
• Proposition 6.1
• Lemma 7.1
• Proposition 7.2