PaperlayerView exact revision on arXiv

Source-linked AI summary

Fast Best Subset Selection: Coordinate Descent and Local Combinatorial Optimization Algorithms

Hussein Hazimeh, Rahul Mazumder

arXiv:1803.01454v3stat.COmath.OCstat.ML

TL;DR

Exact algorithms incur a steep computational price, motivating more efficient approaches to sparse learning. The paper introduces hierarchical optimality conditions and efficient algorithms, with Algorithm 2 matching global MIO solutions while running much faster.

  • Problem

    Exact algorithms for sparse learning are accompanied by a steep computational price.

  • Method

    The paper introduces a hierarchy of necessary optimality conditions and develops efficient algorithmic counterparts for finding higher-quality local minima.

  • Results

    Algorithm 2 solutions match those of global MIO solvers for Problem (1), while running much faster; Algorithm 1 produces better solutions and is faster than IHT and random CD.

  • Takeaways & Limitations

    Higher classes in the proposed hierarchy correspond to better-quality local minima, while the algorithms improve optimization performance relative to the cited alternatives.

  • Takeaways & Limitations

    Prediction errors for the proposed methods can potentially be improved by using Algorithm 2.

Abstract

from arXiv · show

The $L_0$-regularized least squares problem (a.k.a. best subsets) is central to sparse statistical learning and has attracted significant attention across the wider statistics, machine learning, and optimization communities. Recent work has shown that modern mixed integer optimization (MIO) solvers can be used to address small to moderate instances of this problem. In spite of the usefulness of $L_0$-based estimators and generic MIO solvers, there is a steep computational price to pay when compared to popular sparse learning algorithms (e.g., based on $L_1$ regularization). In this paper, we aim to push the frontiers of computation for a family of $L_0$-regularized problems with additional convex penalties. We propose a new hierarchy of necessary optimality conditions for these problems. We develop fast algorithms, based on coordinate descent and local combinatorial optimization, that are guaranteed to converge to solutions satisfying these optimality conditions. From a statistical viewpoint, an interesting story emerges. When the signal strength is high, our combinatorial optimization algorithms have an edge in challenging statistical settings. When the signal is lower, pure $L_0$ benefits from additional convex regularization. We empirically demonstrate that our family of $L_0$-based estimators can outperform the state-of-the-art sparse learning algorithms in terms of a combination of prediction, estimation, and variable selection metrics under various regimes (e.g., different signal strengths, feature correlations, number of samples and features). Our new open-source sparse learning toolkit L0Learn (available on CRAN and Github) reaches up to a three-fold speedup (with $p$ up to $10^6$) when compared to competing toolkits such as glmnet and ncvreg.

1 Introduction

The paper develops fast algorithms for L0-based regression with optional L1 or L2 regularization, targeting computationally challenging high-dimensional problems while studying statistical performance.

  • Motivation: Low-SNR regimes can cause vanilla L0 penalization to overfit, motivating additional L1 or L2 regularization.λ0 controls selected-variable count, while λq controls shrinkage.
  • Motivation: L0-based estimators can offer superior variable selection, prediction, and estimation compared with computationally friendlier sparse-learning alternatives.The comparison is stated for suitable choices of λ0 and λq.
  • Computational challenge: Problem (1) is NP-hard, and existing exact approaches face steep computational costs relative to fast proxy methods.MIO methods address p around 1000 in minutes, whereas Lasso solvers can handle much larger problems within a second.
  • Approach: The proposed approach combines cyclic coordinate descent with local combinatorial search to obtain solutions that cannot be improved by small support changes.The paper also establishes convergence guarantees for these algorithms.
  • Computational results: Local search runs in seconds to minutes for p from 10^3 to 10^6, while L0Learn reports speedups of 25% to 300% over glmnet and ncvreg.The toolkit is open-source and implemented in R/C++.
  • Statistical results: Experiments find solutions often similar to exact MIO solutions with shorter runtimes, alongside superior combined estimation, prediction, and variable-selection performance.The reported comparisons cover real and synthetic datasets.

2 Necessary Optimality Conditions

The paper organizes necessary optimality conditions for L0-regularized problems into increasingly restrictive local-minimum classes and develops corresponding swap-based interpretations.

  • Partial swap minima: PSI(k) minima prohibit improvements from removing or adding up to k support variables followed by optimization over the added subset.As k increases, PSI(k) classes shrink toward global minimizers.
  • Full swap minima: FSI(k) minima use the same support swaps as PSI(k), but optimize over the whole new support.For sufficiently large k, FSI(k) and PSI(k) coincide with global minimizers.
  • Hierarchy: The hierarchy becomes more restrictive from stationary solutions toward FSI(k), with higher classes representing better-quality local minima.For fixed k, PSI(k) contains FSI(k).
  • Stationary solutions: Stationary solutions require nonnegative lower directional derivatives in every direction, equivalently zero gradient components on the active support.Stationary solutions are the weakest class in the hierarchy.
  • Coordinate-wise minima: Coordinate-wise minima cannot improve the objective by optimizing one coordinate while holding all others fixed.The stationary-solution class contains CW minima, generally strictly.
  • Algorithms: The proposed algorithms use cyclic coordinate descent and local combinatorial optimization to obtain PSI(k) and FSI(k) minima, with convergence guarantees.The local search problems are formulated as structured MIO problems, with efficient solvers for special cases.

3 Algorithms

The paper develops cyclic coordinate descent and local combinatorial algorithms for L0-regularized problems, with convergence guarantees to increasingly strong local optimality conditions.

  • 3.1 Cyclic Coordinate Descent: Cyclic coordinate descent with full coordinate minimization is the paper’s main algorithmic workhorse for L0-regularized problems.Spacer steps are occasionally used to stabilize behavior and ensure stationarity of limit points.
  • 3.1 Cyclic Coordinate Descent: Using a conservative step size can harm solution quality because its stationary solutions include CW minima.The resulting class of stationary solutions is broader than the target CW-minimum set.
  • 3.1.1 Convergence Analysis: The convergence guarantees rely on minor assumptions, including linear independence of every set of min{n, p} columns and an initialization condition when p > n.Under these assumptions, coordinate-descent iterates have support size at most min{n, p}.
  • 3.1.1 Convergence Analysis: Theorem 2 proves finite support stabilization and convergence of the iterates to a CW minimum.The support becomes constant after finitely many iterations, and the sequence converges to a limit with that support.
  • 3.1.1 Convergence Analysis: Algorithm 1 has an asymptotic linear convergence rate, extending a result for smooth convex minimization to this non-convex setting.The result applies to full-cycle iterates and uses the stabilized-support behavior.
  • 3.2 Local Combinatorial Optimization: Algorithm 2 terminates in finitely many iterations and returns a PSI(k) minimum for the L0 and related problems.The local combinatorial formulation searches over limited support changes around the current solution.
  • 3.2 Local Combinatorial Optimization: The combinatorial subproblem has a smaller search space than the full MIO formulation, making small-k solves usually faster.An MIO framework can also quickly produce a feasible solution to the restricted problem.

4 Efficient Computation of the Regularization Path

The implementation accelerates regularization-path computation through continuation, coordinate ordering, screening, active sets, and sparsity-aware updates, enabling high-dimensional use.

  • Toolkit and computational strategies: L0Learn is an extensible C++ toolkit with an R interface implementing most of the paper’s algorithms.It is available through CRAN and Github.
  • Adaptive Selection of Tuning Parameters: Continuation uses a decreasing grid of λ0 values, warm-starting each solution from the preceding value while avoiding duplicate solutions.The paper derives conditions for selecting λ0 values that produce different outputs.
  • Coordinate ordering: Partially greedy cyclic ordering sorts coordinates by decreasing residual correlation, while partial sorting can reduce work when p is very large.Only the top t coordinates need to be sorted; the remainder retain their initial order.
  • Correlation Screening: Correlation screening restricts updates to the warm-start support and highly residual-correlated coordinates, then checks excluded coordinates for CW violations.The screened solution typically satisfies the CW conditions after one pass over all coordinates.
  • Active Set Updates: Active-set updates exploit finite support stabilization; empirically, Algorithm 1 typically stabilizes its support in fewer than 10 full cycles.After repeated unchanged cycles, updates are restricted to the current support before an outside-support check.
  • Fast Coordinate Updates: Covariance updates cost O(∥βk∥0) per coordinate after preprocessing and are cheaper than residual updates when ∥βk∥0 < n.Residual updates instead use O(n) operations and can benefit from sparse unchanged residuals.

5 Computational Experiments

The experiments compare the proposed coordinate-descent and local combinatorial methods with sparse-learning baselines across optimization, correlation, sample-size, signal-strength, and high-dimensional settings. The results show that optimization quality, regularization, and problem difficulty determine when each method has an advantage.

  • Optimization performance: Greedy Cyclic CD improves the objective by roughly 12% over Random CD and 55% over IHT.In the same comparison, it requires roughly 28% fewer iterations than Random CD and 75% fewer than IHT.
  • Varying sample size: Algorithm 2 performs best in highly correlated settings, dominating prediction error and the L∞ norm while outperforming Algorithm 1.The reported advantage is attributed to its local combinatorial optimization procedures.
  • Varying sample size: In the mild-correlation setting, Algorithm 1 (L0) and Algorithm 2 (L0L2) have similar profiles, outperform competing methods, and fully recover support with roughly 300 samples.The swap variants do not provide significant improvements in this easier setting.
  • Local-minimum comparisons: FSI minima reduce false positives relative to PSI minima, while coordinate descent increases true positives and local combinatorial search reduces false positives.FSI minima require solving relatively more difficult optimization problems but remain within reasonable computation times.
  • Statistical performance: For high-SNR synthetic settings, L0L1 and L0L2 fully recover the true support and attain the lowest prediction error, whereas Lasso and Relaxed Lasso use many false positives.In another setting, Lasso, Relaxed Lasso, and Elastic Net never achieve full support recovery.
  • Large high-dimensional experiments: L0Learn is faster than glmnet and ncvreg on all considered datasets, exceeding twice their speed on the Amazon Reviews dataset.For L0L2 and L0L1, it also provides much sparser supports with competitive testing MSE.

6 Conclusion

The paper develops cyclic coordinate-descent and local combinatorial-search algorithms for L0-based problems, establishes convergence properties, and evaluates their optimization and statistical behavior. These methods can match global MIO solutions in difficult settings while reducing computation and revealing regime-dependent differences from popular sparse-learning methods.

  • Contributions: The proposed algorithms combine cyclic coordinate descent with local combinatorial search and are analyzed through a hierarchy of necessary optimality conditions.Solutions higher in the hierarchy have higher quality.
  • Optimization performance: Algorithm 1 produces better and faster solutions than IHT and random coordinate descent, while Algorithm 2 often improves further.The local-search procedure extends the coordinate-descent solution with additional combinatorial optimization.
  • Optimization performance: In many difficult settings, Algorithm 2 matches global MIO solutions for Problem (1) while running much faster.This demonstrates competitive optimization quality without the computational cost of exact MIO methods.
  • Statistical properties: At low signal strength or small sample sizes, L0L2, Elastic Net, and ridge typically achieve the best predictive performance, with L0L2 selecting much smaller supports.The paper characterizes these regimes as settings where full support recovery or very small estimation error may be impossible.
  • Statistical properties: Across varied settings, L0-based estimators typically outperform state-of-the-art sparse-learning algorithms on prediction, variable-selection, and estimation metrics.The paper reports no overall winner among vanilla Lasso, stepwise selection, and L0; modifications such as Problem (1) or Relaxed Lasso appear necessary.
  • Implementation: L0Learn provides an open-source implementation that achieves up to a 3x speed-up over competing sparse-learning toolkits.The comparison includes toolkits such as glmnet and ncvreg.

A.4 Proof of Lemma 5

The proof shows that every iteration produces a non-increasing objective value. Because the objective is bounded below, the sequence of objective values converges.

  • Descent property: Non-spacer steps do not increase the objective, and spacer steps preserve this descent through simultaneous loss decrease and non-increasing support size.Thus F(β^k) ≤ F(β^(k−1)) for either type of step.
  • Convergence: The objective sequence is non-increasing and bounded below by zero, so it converges to some finite F*.This establishes objective-value convergence, not necessarily convergence of the iterates themselves.

A.5 Proof of Lemma 6

The proof bounds every iterate’s support size by the sample size n when p > n. It applies the descent property to both the L0 and L0L1 problems.

  • L0 problem: For the L0 problem, the descent bound implies λ0||β^k||0 ≤ λ0n, hence every iterate satisfies ||β^k||0 ≤ n.The argument uses the initial objective bound F(β^0) ≤ λ0n.
  • L0L1 problem: For the L0L1 problem, comparison with the lasso optimum yields f(β^k) − f(β^ℓ1) ≥ 0 and therefore ||β^k||0 ≤ n.The lasso solution’s optimality supplies the nonnegative loss-difference term needed for the bound.

A.6 Proof of Theorem 2

The proof establishes boundedness and characterizes the limit behavior of Algorithm 1. Repeatedly generated supports eventually stabilize, with iterates converging to a unique stationary solution on that support.

  • The sequence generated by Algorithm 1 is bounded.
  • For any support generated infinitely often, the iterates eventually have exactly that support.
  • A subsequence associated with such a support converges to a stationary solution.
  • Every subsequence with the same support converges to the same limit point.
  • Strong convexity on a fixed support makes this stationary solution unique.
  • Every limit-point support appears infinitely often, ensuring that support-specific convergence results apply to the sequence’s limit points.

A.7 Proof of Theorem 3

The proof analyzes Algorithm 1 after its support and signs stabilize. It then applies linear convergence of cyclic coordinate descent to the resulting strongly convex smooth objective.

  • After stabilization, Algorithm 1 updates coordinates on a fixed support with fixed signs.
  • The stabilized iterates are generated by coordinate descent with step size 1/(1+2λ2).
  • The restricted objective has Lipschitz gradient parameter L = M_S + 2λ2 and strong-convexity parameter σ_S = m_S + 2λ2.
  • Known results for cyclic coordinate descent on strongly convex, continuously differentiable functions yield a linear convergence rate.

A.8 Proof of Theorem 4

The proof shows that Algorithm 2 terminates after finitely many iterations. At termination, its output is a PSI(k) minimum.

  • Each iterate produced by Algorithm 1 is a coordinatewise minimum.
  • Distinct iterates of Algorithm 2 must have distinct supports because coordinatewise minima on one support share the same objective.
  • Because only finitely many supports exist, Algorithm 2 terminates in a finite number of iterations.
  • Termination occurs when no feasible solution improves the current objective value.
  • The terminating iterate is therefore a minimizer of the local problem and a PSI(k) minimum.

B.1 Statistical Performance for Varying Number of Samples

The experiments vary the sample size and compare proposed algorithms with one another and with state-of-the-art methods. Relative Risk is defined using the covariance-weighted estimation error.

  • Sample-size experiments vary n from 100 to 1000.
  • The top comparison includes Algorithm 2 with L0, Algorithm 2 with L0L2, and other state-of-the-art algorithms.
  • The bottom comparison evaluates all proposed algorithms.
  • Additional experiments vary SNR from 0.01 to 100 and compare two proposed methods with state-of-the-art algorithms.
  • Relative Risk is the covariance-weighted squared estimation error between the estimator and β†.

C.1 Statistical Performance for Varying Number of Samples

Figures 11–14 describe performance comparisons across varying sample sizes and signal-to-noise ratios. The comparisons include proposed methods and state-of-the-art algorithms.

  • Figure 11 varies the number of samples n from 100 to 1000.
  • Figure 11 compares Algorithm 2 (L0), Algorithm 2 (L0L2), and other state-of-the-art algorithms.
  • Figures 12–14 vary the signal-to-noise ratio (SNR) from 0.01 to 100.
  • Figures 12–14 compare two proposed methods with other state-of-the-art algorithms.

D Comparisons among PSI(k) and FSI(k)

The FSI(5) variant alternates coordinate-descent refinement with MIO-based swap optimization. The reported evolution shows gains in objective values and complementary improvements in true and false positives.

  • D Comparisons among PSI(k) and FSI(k): The FSI(5) variant starts from a CW minimum and alternates CD-PSI(1) with MIO optimization of problem (25).
  • D Comparisons among PSI(k) and FSI(k): Figure 15 marks PSI(1) minima with red circles and MIO solutions with blue squares.
  • D Comparisons among PSI(k) and FSI(k): CD-PSI(1) applied on top of MIO solutions produces important gains in objective values in most cases.
  • D Comparisons among PSI(k) and FSI(k): MIO can produce solutions unavailable through PSI(1), while CD-PSI(1) increases true positives and MIO improves false positives.
  • D Comparisons among PSI(k) and FSI(k): PSI(1) minima obtain a good number of true positives but suffer from false positives.
Loading 1803.01454v3…