Source-linked AI summary
Coordinate descent algorithms for nonconvex penalized regression, with applications to biological feature selection
Patrick Breheny, Jian Huang
TL;DR
Nonconvex penalties such as SCAD and MCP offer attractive variable-selection properties, but their fitting is difficult and can produce unstable solutions. The paper develops coordinate descent algorithms, convergence diagnostics, and related tools for these models, finding strong computational performance and evidence that MCP is often preferred to lasso and SCAD.
Problem
SCAD and MCP have attractive theoretical properties, but fitting their nonconvex objectives is difficult and solutions may be unstable.
Method
The paper fits SCAD and MCP with coordinate descent, analyzes convergence, and uses convexity diagnostics and pathwise optimization for model fitting.
Results
Coordinate descent was up to 1000 times faster than LLA and at least 100 times faster with highly correlated covariates; MCP and SCAD were worthwhile alternatives to lasso, with MCP often preferred.
Takeaways & Limitations
MCP and SCAD are practical alternatives to lasso for many applications, while MCP often offers the preferred combination of predictive accuracy and sparsity.
Takeaways & Limitations
The algorithms are not generally guaranteed to reach a global minimum for nonconvex objectives, and logistic-regression convergence is not guaranteed.
Abstract
from arXiv · showhide
A number of variable selection methods have been proposed involving nonconvex penalty functions. These methods, which include the smoothly clipped absolute deviation (SCAD) penalty and the minimax concave penalty (MCP), have been demonstrated to have attractive theoretical properties, but model fitting is not a straightforward task, and the resulting solutions may be unstable. Here, we demonstrate the potential of coordinate descent algorithms for fitting these models, establishing theoretical convergence properties and demonstrating that they are significantly faster than competing approaches. In addition, we demonstrate the utility of convexity diagnostics to determine regions of the parameter space in which the objective function is locally convex, even though the penalty is not. Our simulation study and data examples indicate that nonconvex penalties like MCP and SCAD are worthwhile alternatives to the lasso in many applications. In particular, our numerical results suggest that MCP is the preferred approach among the three methods.
2.1. MCP.
MCP reduces the lasso’s bias through a nonconvex penalty whose shrinkage relaxes for larger coefficients. Its coordinate descent updates use the univariate MCP solution applied to each predictor’s partial residuals.
- MCP penalty: MCP begins with lasso-level penalization but continuously relaxes it until the penalty derivative reaches zero for sufficiently large coefficients.The relaxation occurs when θ > γλ.
- Univariate solution: For a univariate regression, MCP uses a closed-form solution based on the soft-thresholding operator.The soft-thresholding operator is the univariate lasso solution.
- Thresholding behavior: As γ increases, MCP approaches lasso soft-thresholding; as γ approaches 1, it approaches hard thresholding.MCP therefore interpolates between soft and hard thresholding in the univariate setting.
- Thresholding behavior: With an orthonormal design, MCP corresponds to multivariate firm-thresholding, between lasso soft-thresholding and subset-selection hard-thresholding.This analogy characterizes the three procedures through their thresholding behavior.
- Coordinate descent: Coordinate descent applies the univariate MCP solution to each coordinate using the regression of its partial residuals on the corresponding predictor.The partial residual is r_-j = y − X_-jβ_-j.
- Coordinate descent: Using current residuals to obtain each coordinate’s unpenalized solution is more computationally efficient than recomputing partial residuals.The residual vector is maintained with current coefficient updates.
2.2. SCAD.
SCAD uses a nonconvex penalty that, like MCP, starts with lasso-level shrinkage and reduces it for larger coefficients. Its coordinate-wise objective is convex under a parameter condition, yielding convergence to coordinate-wise and local minima, but not generally to the global minimum.
- SCAD penalty: SCAD begins with the lasso’s penalization rate and reduces that rate to zero as the coefficient moves farther from zero.The transition differs from MCP’s transition, although both penalties reduce shrinkage for larger coefficients.
- Univariate solution: The SCAD univariate solution rescales the soft-thresholding solution toward the unpenalized solution, similarly but not identically to MCP.Unlike MCP, SCAD does not approach hard thresholding at its lower γ limit.
- Convergence: For SCAD with γ > 2 and MCP with γ > 1, the objective is convex in each coordinate when the remaining coefficients are fixed.This coordinate-wise convexity supports the convergence analysis.
- Convergence: The coordinate descent coefficient sequence converges to a point that is both a local minimum and a global coordinate-wise minimum.The convergence result applies to the SCAD and MCP algorithms.
- Global convexity: Because SCAD and MCP are nonconvex, their coordinate descent algorithms are not generally guaranteed to reach a global minimum.Global convergence is obtained when the full objective is convex under the stated γ and design conditions.
2.4. Pathwise optimization and initial values.
The pathwise procedure fits penalized logistic models across decreasing λ values, using previous solutions as initial values and iteratively updating quadratic approximations. This approach has practical safeguards, but convergence is not guaranteed for logistic regression.
- Pathwise optimization: For MCP and SCAD, λmax is the smallest penalty value at which all penalized coefficients are zero.The path can begin at λmax with β(0) = 0 and proceed toward smaller λ values.
- Pathwise optimization: Solutions are computed over 100 λ values equally spaced on the log scale.The procedure follows the pathwise strategy used for lasso fitting.
- Logistic regression: In logistic regression, the model includes an unpenalized intercept alongside standardized predictors and binary responses.The intercept is represented by an unpenalized column of ones.
- Logistic regression: Each logistic-regression iteration first forms a quadratic approximation to the loss using a Taylor expansion around the current coefficients.This produces the iteratively reweighted least-squares form with a working response and diagonal weight matrix.
- Iterative fitting: The approximation and one LLA or coordinate-descent iteration are repeated until convergence for each λ value.Coordinate descent updates all covariates and the intercept, with λ set to zero for the intercept.
- Generalized linear models: The same reweighting framework can extend to penalized generalized linear models when the working response, weights, and residuals are adapted to the response distribution and link.The paper focuses on logistic regression but states this broader extension condition.
- Convergence caveat: Neither LLA nor coordinate descent is guaranteed to converge for logistic regression, although safeguards against model saturation prevented observed convergence failures.The authors report this based on their numerical experience.
3.1. Fixed scale solution.
The fixed-scale coordinate-descent updates extend linear-regression solutions to MCP and SCAD, but logistic regression introduces covariate-dependent reweighting that can make tuning difficult and solutions unstable.
- MCP updates: MCP coordinate-descent updates use piecewise rules based on the current residual regression and the penalty parameters.The displayed update is followed by a condition on γ and residual updating analogous to linear regression.
- Logistic-regression reweighting: Logistic-regression updates reweight portions of both the numerator and denominator, unlike linear regression where v_j equals 1.For linear regression, the reweighting term therefore disappears from the solution.
- Parameter interpretation: For standardized covariates in linear regression, MCP and SCAD scaling factors have constant, universal interpretations across λ and regression problems.The constants are 1−1/γ for MCP and 1−1/(γ−1) for SCAD.
- Parameter interpretation: In logistic regression, the scaling factor changes across covariates, making γ difficult to choose and removing its consistent interpretation.This issue is tied to the continually changing logistic-regression scaling factor.
- Stability: With γ = 3.7 in logistic regression, the solution becomes discontinuous and unstable because the objective function is poorly behaved.The instability is presented as a consequence of the objective function rather than of a particular algorithm.
- Rescaling: Adaptive rescaling replaces pλ,γ(|β_j|) with pλ,γ(|v_jβ_j|) to match changing covariate scales.The resulting coordinate-descent steps are simple extensions of the linear-regression solutions.
3.2. Adaptive rescaling.
Adaptive rescaling adjusts γ to changing covariate scales, but its interpretation differs across models and can be especially difficult in logistic regression.
- Rescaling mechanism: Adaptive rescaling changes the penalty argument from |β_j| to |v_jβ_j| so γ matches the changing scale of each covariate.The rescaling affects the range over which the penalty applies, not λ's magnitude.
- Logistic-regression behavior: For covariates associated with fitted probabilities near 0 or 1, adaptive rescaling extends the range over which penalization is applied.The paper describes this as reasonable because large coefficient changes may produce only small changes in model fit.
- SCAD: Adaptively rescaled SCAD does not generally equal ordinary SCAD with different shape parameters because its penalty scale is tied to the coefficient scale.The paper suggests a more flexible SCAD parameterization as a possible remedy.
- Practical consequence: Adaptive rescaling increases interpretability and makes γ easier to select, according to the paper's Section 5 results.
4.1. Convexity and stability.
Convexity provides global convergence and smooth solution paths, but in high-dimensional problems local convexity around sparse solutions can be sufficient.
- Convexity conditions: The logistic objective is convex where c*(β) > 1/γ for MCP or c*(β) > 1/(γ−1) for SCAD.Here c*(β) is the minimum eigenvalue of n−1X′WX.
- Benefits of convexity: Convexity ensures that coordinate descent converging to a critical point reaches the unique global minimum.
- Benefits of convexity: Convexity also makes β̂ continuous in λ, providing good initial values and reducing coordinate-descent iterations.
- Stability: Without convexity, small data changes may produce large estimate changes, and such estimators tend to have high variance.
- High-dimensional settings: When p > n, global convexity is impossible or irrelevant, but local convexity around sparse solutions can still yield stable estimates and smooth paths.The relevant region is the parameter-space neighborhood containing the sparse solutions of interest.
4.2. Local convexity diagnostics.
The paper introduces diagnostics for locating locally convex portions of nonconvex solution paths and illustrates that these regions have smoother, more stable estimates than nonconvex regions.
- Diagnostic purpose: The proposed diagnostic identifies coefficient-path regions that retain convexity benefits despite SCAD or MCP penalties being nonconvex.
- Diagnostic construction: The local-convexity cutoff is based on the active and imminently active covariates through c*(λ).The active-set construction uses U(λ)=A(λ)∪A(λ−) and the corresponding design matrix.
- Diagnostic construction: The objective is locally convex on (∞, λ*) and locally nonconvex on [λ*, 0], with λ* occurring when the active set changes.
- Empirical illustration: In simulated data with n = 20 and p = 50, MCP paths were smooth and well behaved in unshaded locally convex regions but discontinuous and noisy in shaded regions.Comparable SCAD behavior is reported, while logistic-regression paths are shown separately in Figure 5.
- Interpretation: Noisy paths in nonconvex regions may reflect suboptimal numerical convergence, statistical variability, or both.The paper therefore argues that practitioners should identify the region containing their solution.
- Tuning parameters: AIC and BIC may select local minima in nonconvex regions, while cross-validation is computationally intensive over γ and λ grids.The paper proposes combining information criteria, cross-validation, and convexity diagnostics when selecting tuning parameters.
- Computational assessment: The efficiency assessment compares coordinate descent and LLA for fitting entire MCP and SCAD paths in linear and logistic regression.
5.1. Computational efficiency.
Coordinate descent and LLA converged to the same coefficient paths, but coordinate descent required dramatically less computation across the investigated settings.
- The LLA and coordinate descent algorithms converged to the same coefficient path within the convergence-criterion accuracy.
- Both algorithms’ computational burden increased approximately linearly with p on the log scale.
- Coordinate descent was up to 1000 times faster than LLA when fitting entire coefficient paths.
- Highly correlated covariates made coordinate descent 2–5 times slower, but it remained at least 100 times faster than LLA.
- The statistical comparison evaluated MCP and SCAD against lasso in simple and biomedical-application-like simulations.
5.2. Comparison of MCP, SCAD and lasso.
MCP and SCAD reduce lasso’s bias for large coefficients, while lasso remains advantageous for small coefficients; in biomedical-style simulations, MCP often achieved sparse models without sacrificing prediction.
- MCP and SCAD greatly outperformed lasso in mean squared error when the underlying coefficients were large.
- Lasso’s shrinkage was beneficial for small coefficients, where MCP and SCAD were more likely to overfit noisy data.
- MCP achieved similar predictive accuracy with fewer selected variables and fewer spurious covariates in the sparse genetic-association simulation.
- In the dense microarray simulation, MCP with γ = 20 retained lasso’s predictive accuracy while improving parsimony.
- High feature correlation in the microarray simulation limited the ability of all methods to identify the generating covariates, despite good predictive accuracy.
6.1. Genetic association.
For the sparse genetic-association problem, convexity diagnostics and BIC supported a small MCP tuning parameter, and MCP produced a sparser model with modestly lower cross-validated error.
- BIC and convexity diagnostics were used to choose γ for the genetic-association analysis.
- γ = 1.5 balanced sparsity and convexity for MCP, while γ = 2.5 was chosen for SCAD.
- Ten-fold cross-validation selected λ for MCP, SCAD, and lasso.
- MCP achieved a modestly lower cross-validated prediction error while producing a much sparser model than SCAD or lasso.
- The study used logistic regression on gene-expression data from 27 ALL and 11 AML patients, with 34 additional patients reserved for testing.
6.2. Gene expression.
The gene-expression results show that MCP can retain comparable classification performance while selecting a more parsimonious model, provided its tuning is chosen appropriately. More broadly, the paper presents MCP as a frequent preferred alternative to lasso and SCAD.
- MCP achieved the same test-set classification as lasso and SCAD, correctly classifying 31/34 observations.
- MCP used 11 predictors, compared with 13 for lasso and SCAD, while maintaining that classification accuracy.This sparsity is important for screening and diagnostic applications where fewer features can reduce assay costs.
- Using γ = 5 tripled the test error to 9/34, showing that MCP performance depends critically on choosing γ appropriately.
- The numerical results suggest that MCP is often the preferred approach among lasso, SCAD, and MCP.
- The proposed algorithms, convexity diagnostics, logistic-regression rescaling, and ncvreg package address practical barriers to using nonconvex penalties.
APPENDIX
The appendix establishes coordinate-wise convexity properties for the penalized objective and uses them to justify convergence of the proposed cyclic coordinate descent algorithms to coordinate-wise minima.
- The objective is not differentiable, but it has directional first and second derivatives at every point and direction.
- For each coordinate, the objective is strictly convex when γ > 1 for MCP and γ > 2 for SCAD.
- Tseng’s sufficient conditions apply because of coordinate-wise strict convexity and continuity, establishing convergence to coordinate-wise minima.
- The full objective is convex where c*(β) > 1/γ for MCP and c*(β) > 1/(γ − 1) for SCAD.