LLADA.CPP: Robust LAD-LASSO Solver

• LLADA.CPP is a compact, efficient framework for solving the LAD-LASSO regression problem using coordinate descent and ternary search.
• It leverages convexity and piecewise-linear properties to guarantee global optimization and robust variable selection.
• The framework is designed for scalability, supporting parallel GPU execution and browser-based modeling for numerous small-scale regressions.

LLADA.CPP is a compact, implementation-oriented algorithmic framework for efficiently solving the Least Absolute Deviation LASSO regression problem (LAD-LASSO) via coordinate descent augmented with a ternary-search “quadrature” approach. Originally designed to support interactive GPU- and browser-based mathematical modeling with robust variable selection, LLADA.CPP delivers a practical and general-purpose solver structure particularly well-suited for simultaneous solution of hundreds of small-scale regression tasks—either on single CPU threads or in parallelized browser GPU environments (Wright, 17 Oct 2025).

1. Mathematical Formulation and Problem Context

LLADA.CPP solves the LAD-LASSO regression objective, a convex optimization problem combining robust regression with L1-penalized variable selection. The task is, given data $X\in \mathbb{R}^{n\times p}$ and response vector $y\in \mathbb{R}^n$, to compute

$$F(\beta) = \sum_{i=1}^n |y_i - x_i^T\beta| + \lambda\sum_{j=1}^p |\beta_j|$$

for regression coefficients $\beta \in \mathbb{R}^p$ and penalty parameter $\lambda\geq 0$. The LAD loss offers superior robustness to outliers compared to OLS, while the L1 penalty provides sparsity for variable selection. This combination is relevant in applications where both statistical interpretability and robustness are required. The convexity of $F$ ensures a unique global minimizer, but the non-smooth, piecewise-linear nature complicates optimization using standard second-order or smooth-gradient methods.

2. Axis-wise Minima, Convex Locus, and Algorithmic Principle

A foundational observation is that—despite multiple axis-wise local minima—the locus of these minima forms a monotonic, convex one-dimensional set. Specifically, if all coordinates but $j$ are fixed, $F$ as a function of $\beta_j$ is convex and piecewise linear, with a unique minimizer. By varying one distinguished coordinate (typically $\beta_1$) and, for each trial value, running cyclic coordinate descent (CCD) over the other $y\in \mathbb{R}^n$0 coordinates (fixing $y\in \mathbb{R}^n$1), one traces a curve of axis-wise minima with the crucial property of monotonicity in every coordinate. Restricting $y\in \mathbb{R}^n$2 to this curve yields a univariate convex function, making global optimization on this set tractable by single-variable convex search (e.g., ternary search).

This structure justifies the two-level approach of LLADA.CPP: an outer ternary search along the distinguished coordinate, and an inner coordinate descent to minimize $y\in \mathbb{R}^n$3 over the remaining axes for each fixed value.

3. Cyclic Coordinate Descent and Coordinate-wise Update Mechanics

For a fixed coordinate $y\in \mathbb{R}^n$4, the update in LLADA.CPP minimizes a univariate penalized LAD subproblem: $y\in \mathbb{R}^n$5 where $y\in \mathbb{R}^n$6 is the partial residual. This function is convex and piecewise linear in $y\in \mathbb{R}^n$7, with breakpoints at $y\in \mathbb{R}^n$8 for all $y\in \mathbb{R}^n$9 where $$F(\beta) = \sum_{i=1}^n |y_i - x_i^T\beta| + \lambda\sum_{j=1}^p |\beta_j|$$0, and $$F(\beta) = \sum_{i=1}^n |y_i - x_i^T\beta| + \lambda\sum_{j=1}^p |\beta_j|$$1 (due to the L1 penalty). The update is performed by evaluating the objective at these candidate breakpoints:

1. Construct $$F(\beta) = \sum_{i=1}^n |y_i - x_i^T\beta| + \lambda\sum_{j=1}^p |\beta_j|$$2.
2. For each $$F(\beta) = \sum_{i=1}^n |y_i - x_i^T\beta| + \lambda\sum_{j=1}^p |\beta_j|$$3, compute $$F(\beta) = \sum_{i=1}^n |y_i - x_i^T\beta| + \lambda\sum_{j=1}^p |\beta_j|$$4.
3. Set $$F(\beta) = \sum_{i=1}^n |y_i - x_i^T\beta| + \lambda\sum_{j=1}^p |\beta_j|$$5 to the $$F(\beta) = \sum_{i=1}^n |y_i - x_i^T\beta| + \lambda\sum_{j=1}^p |\beta_j|$$6 which minimizes $$F(\beta) = \sum_{i=1}^n |y_i - x_i^T\beta| + \lambda\sum_{j=1}^p |\beta_j|$$7.

The coordinate-wise update is $$F(\beta) = \sum_{i=1}^n |y_i - x_i^T\beta| + \lambda\sum_{j=1}^p |\beta_j|$$8 in the worst case but practical for moderate $$F(\beta) = \sum_{i=1}^n |y_i - x_i^T\beta| + \lambda\sum_{j=1}^p |\beta_j|$$9 ($\beta \in \mathbb{R}^p$0) and highly suitable for parallelization in GPU settings.

4. Ternary-Search Quadrature and Global Minimization

The outer optimization in LLADA.CPP proceeds via a ternary search over a selected coordinate (commonly $\beta \in \mathbb{R}^p$1 or $\beta \in \mathbb{R}^p$2). An initial bracket $\beta \in \mathbb{R}^p$3 is established (e.g., via a doubling search) to ensure the global minimizer lies within. The ternary search iteratively refines this bracket:

• Evaluate $\beta \in \mathbb{R}^p$4 at $\beta \in \mathbb{R}^p$5 and $\beta \in \mathbb{R}^p$6 by fixing $\beta \in \mathbb{R}^p$7 at $\beta \in \mathbb{R}^p$8 or $\beta \in \mathbb{R}^p$9 and performing inner CCD on the other coordinates.
• Retain the bracket half that contains the lower value of $\lambda\geq 0$0.
• Terminate when the bracket width is below the preset tolerance.

The solution is finalized via one complete CCD pass at the midpoint of the last bracket. All steps maintain convexity guarantees: the minimizer found is global for $\lambda\geq 0$1 (Wright, 17 Oct 2025).

5. Initialization, Convergence, and Hyperparameters

Initialization uses $\lambda\geq 0$2 as the starting point (always feasible). Convergence in the inner CCD loop is determined by a coordinate-wise tolerance or a maximum number of cycles. Outer ternary search accuracy is set by the bracket reduction tolerance (e.g., $\lambda\geq 0$3 to $\lambda\geq 0$4). The penalty parameter $\lambda\geq 0$5 can be selected for robustness ($\lambda\geq 0$6) or for variable selection (via $\lambda\geq 0$7-fold cross-validation from $\lambda\geq 0$8 to $\lambda\geq 0$9). It is essential to enforce $F$0 to address non-uniqueness in flat stretches of the objective.

6. Implementation Structure and Computational Considerations

LLADA.CPP is designed for ease of translation to C++ and by extension to WebGL-style GPU kernels. The data matrix $F$1 is stored as a 1D row-major array, with auxiliary vectors for $F$2, $F$3, and residuals. All primary routines are simple for-loops, and all memory management is explicit—no dynamic allocations inside core loops. Candidate sets are stored in preallocated buffers of size $F$4. No external dependency is required beyond basic standard libraries. The resulting implementation is efficient on single CPU threads and trivially parallelizable for GPU execution over multiple independent regression tasks.

A minimal C++ structure for LLADA.CPP, covering all main mechanics and control flows, can be constructed directly from the algorithmic summary above. It is fully functional while leaving room for further enhancement, such as cross-validation loops, multi-threaded parallelization, or additional logging (Wright, 17 Oct 2025).

7. Practical Scope, Extensions, and Applications

LLADA.CPP enables simultaneous robust regression and feature selection in contexts where extreme scale (hundreds of tasks) or deployment constraints (browser, WebGL, or single-threaded platforms) prohibit heavyweight optimization toolkits. The algorithm’s structure is suitable for GPU parallelization, as each coordinate-update subproblem is computationally simple and each regression instance is independent. Extensions may include more efficient one-dimensional minimization for coordinate updates, adaptive convergence tolerances, or custom initialization schemes. The outlined approach is central to interactive mathematical modeling, especially in decision support tools requiring live model fitting and interpretation with outlier robustness.

Potential avenues for further research or implementation enhancement include exploring advanced memory layouts, alternative regularization schemes, and integration into distributed computation frameworks targeting low-latency environments.

---

For detailed algorithmic explanation and pseudocode suitable for direct implementation, see "A simple LAD-LASSO coordinate descent algorithm for interactive browser-based GPU applications" (Wright, 17 Oct 2025).