Coordinate Descent Algorithm for Least Absolute Deviations Regression

1. Problem Statement

The paper addresses the computational limitations of Least Absolute Deviations (LAD) regression, a robust alternative to Ordinary Least Squares (OLS) that minimizes the sum of absolute residuals rather than squared residuals.

• The Challenge: While LAD is robust to outliers and heavy-tailed noise, traditional solvers rely on Linear Programming (LP) formulations (e.g., Simplex or Interior-Point methods). These methods become computationally expensive, numerically unstable, or fail entirely in high-dimensional settings where the number of predictors ($p$) exceeds or equals the number of observations ($n$).
• The Gap: There is a lack of scalable, matrix-inversion-free algorithms that can handle the non-differentiable nature of the LAD objective function efficiently, particularly when $p \ge n$.

2. Methodology

The authors propose a Coordinate Descent (CD) algorithm tailored specifically for LAD regression. The method leverages the convexity of the LAD objective and the specific structure of coordinate-wise updates.

A. Core Mathematical Insight

The LAD objective function is:
$$L(\beta) = \sum_{i=1}^n |y_i - x_i^T \beta|$$
When all coefficients except $\beta_j$ are fixed, the problem reduces to a univariate minimization:
$$\beta_j = \arg\min_{t} \sum_{i=1}^n |r_i^{(j)} - x_{ij}t|$$
where $r_i^{(j)}$ is the partial residual.

• Closed-Form Solution: This subproblem admits an exact solution given by the weighted median of the ratios $r_i^{(j)}/x_{ij}$ (with weights $|x_{ij}|$). If $x_{ij}=1$ (intercept), it is the standard median.
• Convergence: Since the objective is convex and each coordinate update is an exact minimization, the algorithm is guaranteed to converge to a global minimizer.

B. Algorithmic Optimization

To overcome the inefficiency of naive coordinate descent (which requires $O(np^2)$ complexity per iteration due to repeated residual recomputation), the authors introduce an Incremental Residual Update strategy:

1. Maintain Global Residual: Keep a single residual vector $res = y - X\beta$.
2. Update Steps:
   • Restore: Temporarily add the contribution of the current $\beta_j$ to get the partial residual ($O(n)$).
   • Solve: Compute the new $\beta_j$ via weighted median ($O(n \log n)$).
   • Update: Adjust the global residual by the change in $\beta_j$ ($O(n)$).

• Complexity: This reduces the per-iteration complexity from $O(np^2)$ to $O(np \log n)$, making it highly scalable for large $p$.

C. Initialization Strategies

Recognizing that CD performance depends on initialization, especially in underdetermined regimes ($p \ge n$), the paper evaluates:

• Ridge Regression Warm-start: Solving a penalized least squares problem first. This provides a stable, deterministic starting point that shrinks coefficients, effectively handling rank deficiency.
• Genetic Algorithm (GA): Used as a global search initializer for highly non-convex or complex landscapes, though it is computationally heavier.
• Random Multi-start: Running the algorithm from multiple random seeds to ensure robustness.

3. Key Contributions

1. Algorithm Design: A novel coordinate descent algorithm for LAD that avoids matrix inversion and handles the non-differentiable objective via weighted medians.
2. Scalability: An optimized implementation with $O(np \log n)$ complexity per iteration, enabling LAD regression in high-dimensional settings ($p \ge n$) where LP solvers fail.
3. Theoretical Guarantees: Proof of global convergence based on the convexity of the LAD objective and the exactness of coordinate updates.
4. Hybrid Framework: Demonstration that combining Ridge regression (for stability) with LAD-CD (for robustness) yields superior predictive performance in contaminated and high-dimensional data.

4. Experimental Results

The authors evaluated the method on synthetic data and real-world benchmarks (Boston Housing, Air Quality, Concrete Compressive Strength).

• Convergence & Stability:
  • The algorithm converges monotonically within 20–50 iterations.
  • In 1,000 replications on synthetic data, the final Mean Absolute Error (MAE) was highly concentrated (SD $\approx 6 \times 10^{-4}$), proving stability regardless of random initialization.
• Robustness to Outliers:
  • In datasets with 20% outlier contamination, LAD-CD significantly outperformed OLS and matched the accuracy of LP-based solvers (QuantReg).
  • Ridge-initialized LAD-CD reduced prediction MAE by 32.5% compared to the initial Ridge fit, and 44% compared to GA-initialized fits in contaminated settings.
• High-Dimensional Performance ($p \ge n$):
  • LP Solvers: Failed to return solutions when $p > n$ or produced trivial overfitting when $p=n$.
  • LAD-CD: Remained stable and well-defined.
  • Ridge Initialization: Consistently outperformed unregularized methods and GA-based starts in high dimensions, achieving the lowest prediction error.
• Real-World Benchmarks:
  • On the Boston Housing, Air Quality, and Concrete datasets, LAD-CD achieved predictive accuracy (MAE) comparable to the established quantreg (LP-based) package.
  • While LP solvers were faster on small datasets, LAD-CD offered a matrix-inversion-free alternative that scales effectively.

5. Significance and Implications

• Practical Utility: The method provides a "solver-free" tool for robust regression that is easy to implement without specialized LP software.
• High-Dimensional Robustness: It fills a critical gap in statistical methodology by enabling robust median regression in regimes where $p \ge n$, a common scenario in modern genomics and finance.
• Computational Efficiency: By reducing complexity from quadratic to near-linear in $n$ (per sweep), it makes LAD regression feasible for large-scale datasets.
• Future Directions: The paper suggests future work on parallelizing the coordinate updates, incorporating explicit regularization (Lasso/Ridge) directly into the CD steps, and establishing finite-sample guarantees under minimal distributional assumptions.

In conclusion, the proposed LAD-CD algorithm offers a scalable, stable, and theoretically sound alternative to traditional LP-based solvers, particularly excelling in high-dimensional and outlier-contaminated environments.