Best Subset Selection via a Modern Optimization Lens

Dimitris Bertsimas; Angela King; Rahul Mazumder

Source-linked AI summary

Best Subset Selection via a Modern Optimization Lens

Dimitris Bertsimas, Angela King, Rahul Mazumder

arXiv:1507.03133v1Jul 11, 2015stat.MEmath.OCstat.COstat.ML

TL;DR

The paper addresses whether best subset selection can be solved at useful problem sizes despite its combinatorial difficulty. It combines discrete first-order methods with MIO, finding provably optimal or near-optimal solutions and reporting better sparsity and predictive performance than Lasso and other methods. The scope includes least-squares and selected least-absolute-deviation experiments.

Problem
Best subset selection is NP-hard, while existing sparse-regression approaches do not provide globally optimal solutions for the discrete problem.
Method
The approach uses discrete first-order methods to generate warm starts for an MIO formulation, with extensions for coefficient side constraints and least absolute deviation loss.
Results
The method solves problems with n in the 1000s and p in the 100s to provable optimality in minutes, and finds near-optimal solutions for n in the 100s and p in the 1000s in minutes.
Takeaways & Limitations
MIO methods provide sparse models with good predictive power and outperform Lasso and other state-of-the-art methods in the reported experiments.
Takeaways & Limitations
The LAD evaluation uses only a few representative examples, while the paper’s main focus is least-squares loss.

Abstract

from arXiv · show

In the last twenty-five years (1990-2014), algorithmic advances in integer optimization combined with hardware improvements have resulted in an astonishing 200 billion factor speedup in solving Mixed Integer Optimization (MIO) problems. We present a MIO approach for solving the classical best subset selection problem of choosing $k$ out of $p$ features in linear regression given $n$ observations. We develop a discrete extension of modern first order continuous optimization methods to find high quality feasible solutions that we use as warm starts to a MIO solver that finds provably optimal solutions. The resulting algorithm (a) provides a solution with a guarantee on its suboptimality even if we terminate the algorithm early, (b) can accommodate side constraints on the coefficients of the linear regression and (c) extends to finding best subset solutions for the least absolute deviation loss function. Using a wide variety of synthetic and real datasets, we demonstrate that our approach solves problems with $n$ in the 1000s and $p$ in the 100s in minutes to provable optimality, and finds near optimal solutions for $n$ in the 100s and $p$ in the 1000s in minutes. We also establish via numerical experiments that the MIO approach performs better than {\texttt {Lasso}} and other popularly used sparse learning procedures, in terms of achieving sparse solutions with good predictive power.

1 Introduction

Best subset selection seeks a parsimonious k-feature linear-regression fit, but its cardinality constraint makes exact optimization NP-hard and historically difficult to scale. The paper proposes combining discrete first-order methods with mixed integer optimization to obtain optimal or near-optimal solutions with practical computational guarantees and extensions.

Best subset selection seeks a k-feature linear-regression fit that preserves explanatory power while imposing sparsity.
The cardinality constraint makes best subset selection NP-hard, and popular implementations such as leaps do not scale beyond p = 30.
Lasso is computationally feasible but can select noise variables, shrink coefficients, and omit true predictors when sufficient data regularity conditions fail.
Existing continuous non-convex approaches can produce good local minimizers but do not provide guarantees on global optimality.
The proposed framework combines discrete first-order methods and MIO to solve best subset regression to optimality or near optimality within reasonable time.It also accommodates polyhedral coefficient constraints and both least-squares and least-absolute-deviation losses.
Experiments report favorable predictive performance and sparsity for best subset selection relative to competing sparse-learning procedures.

2 Mixed Integer Optimization Formulations

The paper reformulates best subset selection as mixed integer optimization, using modern solver advances, relaxations, and warm starts to obtain optimal or near-optimal solutions with certificates. The formulation supports data-driven parameter bounds and tracks convergence through improving upper and lower bounds.

2.1 Brief Background on MIO: MIO combines discrete and continuous variables in quadratic optimization, providing a framework for best subset selection.The paper introduces MIQO as the relevant subclass of mixed integer optimization.
2.1 Brief Background on MIO: 200 billion: combined hardware and software advances produced this approximate MIO speedup from 1993 to 2013.The cited speedup concerns mixed integer linear optimization, with similar factors reported for MIQO.
2.1 Brief Background on MIO: MIO solvers provide feasible solutions and lower bounds, yielding increasingly tight suboptimality certificates when optimization stops early.This distinguishes MIO from heuristic methods, which do not provide such certificates.
2.2 MIO Formulations for the Best Subset Selection Problem: Best subset selection is reformulated as an MIQO with binary variables controlling which regression coefficients are nonzero.The basic formulation uses coefficient bounds, while SOS-1 constraints provide an alternative without requiring MU.
2.2 MIO Formulations for the Best Subset Selection Problem: Adding bounds on coefficient and norm quantities typically strengthens MIO performance, especially the lower-bound certificates.MU bounds coefficient magnitudes, while Mℓ bounds the ℓ1-norm of β.
2.2 MIO Formulations for the Best Subset Selection Problem: SOS-1 formulation (8) can obtain global solutions without specifying the coefficient bound MU in advance.The formulation models the condition that at most one of βi and 1 − zi is nonzero.
2.2 MIO Formulations for the Best Subset Selection Problem: Figure 2 tracks upper bounds, lower bounds, and the MIO gap for diabetes instances with n = 350, p = 64, and k = 6 or k = 7.Optimal solutions appear within seconds, whereas certification takes 10-20 minutes and becomes slower as k increases.
2.3.1 Specification of Parameters in terms of Coherence and Restricted Strong Convexity: Theorem 2.1 derives parameter bounds from cumulative coherence and restricted eigenvalues, while discrete first-order warm starts provide practical bounds.The hybrid estimate sets MU to τ∥β̂hyb∥∞, with τ greater than one, and derives other bounds from it.

3 Discrete First Order Algorithms

The paper extends first-order convex optimization ideas to cardinality-constrained problems, using hard thresholding and related updates to obtain stationary or near-optimal solutions. Under stated conditions, the algorithm has convergence guarantees and can warm-start the MIO solver.

The discrete first-order solutions are used as warm starts that enhance the performance of the MIO formulation.
The method extends projected gradient descent to smooth convex optimization with cardinality constraints and also adapts to nonsmooth LAD loss.
The cardinality-constrained quadratic subproblem is solved by retaining the k largest-magnitude entries and setting the rest to zero.This defines the hard-thresholding operation H_k(c).
When L > ℓ and α_k > 0, the support stabilizes after finitely many iterations, while the bounded sequence converges to a first-order stationary point.
Theorem 3.1 provides an iteration-dependent guarantee for reaching an ε-approximate first-order stationary point over convex functions with Lipschitz-continuous gradients.
After support stabilization, strong convexity on the stabilized support yields a linear convergence rate without modifying the algorithm.

4 A Brief Tour of the Statistical Properties of Problem (1)

The paper reviews theoretical and empirical properties of best subset solutions across predictive performance, coefficient estimation, and variable selection. These properties are summarized elsewhere in the paper’s review section.

Best subset solutions have been studied using predictive performance, regression coefficient estimation, and variable selection metrics.

5 Computational Experiments for Subset Selection with Least Squares Loss

The experiments evaluate MIO with discrete first-order warm starts on synthetic and real least-squares problems, including both n > p and p ≫ n settings. Across these studies, MIO obtains strong sparse-model and predictive results, often improving warm-start solutions and certifying or approaching optimality.

Experimental settings: The experiments cover synthetic and real datasets in both overdetermined and high-dimensional settings for least-squares subset selection.The study includes n > p and p ≫ n regimes, with examples varying covariance structure, sparsity, signal-to-noise ratio, and coefficient patterns.
Upper bounds and warm starts: MIO with discrete first-order warm starts found the best possible relative accuracy in a few minutes.The first-order solutions were used to initialize MIO, combining fast feasible solutions with MIO improvement and certification.
Statistical performance: MIO generally selected sparser models than competing procedures while matching or improving predictive performance across the synthetic examples.For Example 1, MIO performed best overall; for Examples 2 and 3 it produced sparser models with similar predictive performance, while Example 4 favored MIO on both criteria.
Upper bounds and warm starts: MIO improved first-order solutions on Leukemia data, with smaller k values generally reaching improved upper bounds faster than larger k values.The time required to stabilize the best upper bound depended on k, and larger subsets were typically harder.
Bounding-box formulation: Bounding boxes around β certified optimality and found better solutions faster than bounding boxes around Xβ in the reported experiments.The Xβ bounding boxes made the problem harder, plausibly because of strong correlations among X's columns; SNR had little apparent runtime effect in those comparisons.

6 Computational Results for Subset Selection with Least Absolute Deviation Loss

The LAD experiments compare MIO, discrete first-order methods, and LAD-Lasso across overdetermined and high-dimensional settings. MIO and discrete first-order methods generally achieve better sparsity and predictive accuracy than ℓ1-regularization schemes, while least-squares subset selection is less suitable under heavy-tailed errors.

Experimental setup: The experiments evaluate LAD subset selection for n = 500, p = 100, n = 50, p = 1000, and n = 500, p = 1000.The setup uses k0 = 5 and ρ = 0.9, covering both overdetermined and high-dimensional cases.
Compared methods: The compared procedures include a discrete first-order method, warm-started MIO, and LAD-Lasso.Training, validation, and testing follow the least-squares experiments, with nonzero counts and prediction accuracy reported for each method.
Results: MIO and discrete first-order methods perform much better than ℓ1-regularization schemes in the reported LAD examples.The LAD-Lasso selects overly dense models and has higher prediction error when weakly correlated signal features compete with highly correlated noise features.
Loss-function comparison: Least-squares subset selection produces inferior models in these examples because the errors follow a heavy-tailed distribution.The LAD experiments use double exponential errors, for which LAD loss is described as appropriate.
Results: MIO provides an edge over Lasso in both sparsity and predictive accuracy for overdetermined and high-dimensional cases.Figure 10 reports the overdetermined comparison, while Figure 11 covers the high-dimensional settings.

7 Conclusions

The paper concludes that modern discrete optimization makes best subset selection practically tractable at useful problem sizes. Its framework combines MIO with discrete first-order methods, provides optimality guarantees, and supports coefficient side constraints and LAD loss.

Conclusions: The MIO approach solves problems with n in the 1000s and p in the 100s in minutes to provable optimality.For n in the 100s and p in the 1000s, it finds near-optimal solutions in minutes.
Conclusions: MIO solutions significantly outperform Lasso in achieving sparse models with good predictive power.The conclusion presents this as the main empirical advantage over other state-of-the-art methods.
Conclusions: Early termination still yields a solution with a guarantee on its suboptimality.This guarantee distinguishes the MIO approach from the other methods discussed in the conclusion.
Conclusions: The framework accommodates side constraints on regression coefficients and extends to best subset selection with least absolute deviation loss.These extensions broaden the formulation beyond the standard least-squares, cardinality-constrained problem.
Conclusions: The evidence suggests that MIO methods are tractable and can improve accuracy and sparsity at higher but still reasonable computational times.The paper identifies increased computational time as the trade-off for these properties.
Computational support: The paper also shows that the convex quadratic subproblems arising in its framework can be solved with small computational cost.This supports the computational practicality of the associated formulations.

A.3 Proof of Proposition 2

This proof derives bounds and optimal solutions for the quadratic subproblems used in the formulation. It uses projection and singular-value arguments to establish the needed inequalities and computational steps.

Proof of Proposition 2: The resulting bounds support the computational treatment of the quadratic formulations in the paper.The proof’s algebraic steps connect projection, singular values, and norm inequalities to the formulation’s bounds.
Proof of Proposition 2: The proof combines intermediate bounds to obtain the proposition’s stated inequalities.Equations (43) and (44) are combined to obtain (15), while (42) is substituted into (40) for another bound.
Proof of Proposition 2: The proof uses the singular value decomposition of X_I to characterize the singular values of related inverse matrices.If d1, ..., dk are the nonzero singular values of X_I, the corresponding singular values of (X_I′X_I)^−1 are expressed through 1/d_i.
Proof of Proposition 2: Projection onto the columns of X_I gives the bound ∥P_I y∥_2 ≤ ∥y∥_2.The proof combines this projection bound with √m∥a∥_2 ≥ ∥a∥_1 to derive further inequalities.
Proof of Proposition 2: The derivation applies to vectors supported on an index set I and uses the rows of A in the matrix inequalities.The proof explicitly considers vectors β_I with zero entries outside I.

B.1 Proof of Proposition 6

The proof establishes convergence properties for the discrete first-order algorithm. It shows decreasing objective values, eventual support stabilization under a separation condition, and convergence to a first-order stationary point.

Proof of Proposition 6: The objective values g(β^m) decrease and converge because the sequence is bounded below by zero.This follows by applying the descent inequality to consecutive algorithm iterates.
Proof of Proposition 6: When α_k > 0, the support of β^m eventually stabilizes after finitely many iterations.The proof argues that changing supports infinitely often would contradict the iterates’ convergence behavior.
Proof of Proposition 6: After support stabilization, Algorithm 1 becomes gradient descent restricted to a fixed support space.The proof then invokes convergence of gradient descent for a convex function over a closed convex set.
Proof of Proposition 6: The sequence β^m converges to β*, a first-order stationary point.This is the stated convergence conclusion after the support sequence becomes constant.
Proof of Proposition 6: If the smallest selected coefficient magnitude tends to zero, the gradient norm converges to zero along the iterates.The proof derives this through the limiting behavior of the k-th largest absolute coefficient.
Proof of Proposition 6: Any limit point β∞ solves the unconstrained optimization problem min g(β).Continuity of the gradient and convergence of the objective values are used to establish this result.

B.2 Proof of Proposition 3

The proof shows that the optimal support selects the k largest absolute values of the coefficients c_i. Setting selected coefficients equal to c_i reduces the objective to the squared magnitudes of excluded coefficients.

B.2 Proof of Proposition 3: The optimal support S contains the indices of the k largest |c_i| values.This follows because setting β_i = c_i on S leaves the objective as the sum of squared excluded |c_i| values.

B.3 Proof of Proposition 7

The proof uses Proposition 6(a) and the definition of H_k to show that η and an arbitrary reference point ˆη must coincide. Therefore, η is the unique element of H_k.

B.3 Proof of Proposition 7: Proposition 6(a), together with g(η) = g(ˆη) and L > ℓ, forces ∥ˆη − η∥_2 = 0.Since ˆη is arbitrary, η is the only element in H_k.

B.4 Proof of Proposition 8

The proof establishes that the k-sparse vector bβ is a first-order stationary point by combining its sparsity with Proposition 6(a). The argument uses equality of objective values and L > ℓ to make the relevant bound vanish.

B.4 Proof of Proposition 8: The k-sparse vector bβ is a first-order stationary point.Its objective equality g(bβ) = g(ˆη), together with L > ℓ, makes the right-hand side of the bound zero.
B.4 Proof of Proposition 8: The sequence g(β^{m+1}) decreases and converges to g(β*).This convergence is attributed to Proposition 6.

B.6 Proof of Proposition 5

This section reviews predictive-risk guarantees and compares subset selection with Lasso-based sparse solutions under restricted eigenvalue-type conditions. It concludes that the predictive-performance gap can be substantial when feature correlations make γ(X) small, while statistical properties of near-global solutions remain an open question.

B.6 Proof of Proposition 5: Under the linear model with iid Gaussian noise, prior work gives high-probability worst-case predictive-performance bounds for subset-selection solutions.The bound holds with probability greater than 1 − exp(−c_1 k log(p/k)) and is independent of X.
B.6 Proof of Proposition 5: Lasso-based predictive-performance bounds depend on a restricted eigenvalue-type condition parameterized by γ(X).The condition constrains β in the union of sets C(S) for supports S of size k.
B.6 Proof of Proposition 5: A significant predictive-performance gap separates subset selection from Lasso-based k-sparse solutions when γ(X) is small.High pairwise feature correlations can make γ(X) small, and the authors state that these bounds complement their Section 5 experiments.
B.6 Proof of Proposition 5: Prior analyses report superior variable-selection properties for subset selection compared with Lasso-based methods.The cited work also studies regression-coefficient errors and statistical properties of the Lagrangian formulation.
B.6 Proof of Proposition 5: The statistical properties of suboptimal solutions to the best-subset problem are left for future work.The paper notes that MIO certificates provide global lower bounds for near-global solutions, but a precise theory for their statistical properties remains unavailable.

D.2 Lasso, Debiased Lasso and MIO

The section compares Lasso, debiased Lasso, and MIO on model sparsity and predictive performance across overdetermined and high-dimensional experiments. MIO generally achieves better predictive models and greater sparsity, with its advantage most pronounced under higher correlation and lower signal-to-noise conditions.

Debiasing: Debiasing typically performs worse than Lasso when unrestricted least squares is applied to the Lasso-selected support.The unrestricted least-squares solution has the same sparsity pattern as the corresponding Lasso model.
Debiasing: The debiased-Lasso variant computes an unrestricted least-squares fit for each Lasso-selected support, yielding feasible solutions with k determined by the selected model.The tuning parameter is selected using predictive performance on a held-out validation set.
Predictive Performance: MIO performs significantly better than debiased Lasso, especially for larger ρ and smaller SNR values.The comparison concerns predictive performance across the numerical experiments.
Sparsity: Debiased Lasso produces less dense models than Lasso but more dense models than MIO in the n = 50, p = 1000 setting.The models are selected using tuning parameters chosen by held-out validation performance.
Predictive Performance: MIO consistently delivers better predictive models than debiased Lasso and ordinary Lasso in the n = 50, p = 1000 experiments.Debiased Lasso performs better than ordinary Lasso in these experiments.