Computing Leapfrog Regularization Paths with Applications to Large-Scale K-mer Logistic Regression

Abstract

High-dimensional statistics deals with statistical inference when the number of parameters or features p exceeds the number of observations n (i.e., $p\gg n$ ). In this case, the parameter space must be constrained either by regularization or by selecting a small subset of $m\le n$ features. Feature selection through $l_1$ -regularization combines the benefits of both approaches and has proven to yield good results in practice. However, the functional relation between the regularization strength $\lambda$ and the number of selected features m is difficult to determine. Hence, parameters are typically estimated for all possible regularization strengths $\lambda$ . These so-called regularization paths can be expensive to compute and most solutions may not even be of interest to the problem at hand. As an alternative, an algorithm is proposed that determines the $l_1$ -regularization strength $\lambda$ iteratively for a fixed m . The algorithm can be used to compute leapfrog regularization paths by subsequently increasing m .

1. Introduction

In statistics and machine learning, a common problem is to learn a model of the form $${\mu }_{\theta }\left(z\right)=\sigma \left({\theta }_1{f}_1\left(z\right)+\dots +{\theta }_m{f}_m\left(z\right)\right),$$

with coefficients $\theta =\left({\theta }_1,\dots ,{\theta }_m\right)\in {ℛ}^m$ and where f_j denotes the j th feature computed from the input data z. In statistics, $\sigma$ is called mean function, whereas in machine learning it is commonly referred to as activation function. The model ${\mu }_{\theta }$ is typically estimated by minimizing a given loss function on a set of n observations.

For instance, in linear regression, the mean function is $\sigma \left(x\right)=x$ and the parameters $\theta$ would be estimated by minimizing the sum of squared residuals, also called ordinary least squares (OLS). In the case of artificial neural networks, ${\mu }_{\theta }$ represents a single output neuron and f_j denotes the jth output of the preceding layers of the network. The activation function $\sigma$ is typically either a sigmoid function or some variant. Neural networks are estimated by minimizing, for instance, the mean squared or cross-entropy error.

Two fundamentally different ways can be found in the literature of how features f_j are estimated. Some machine learning methods, such as neural networks, treat each f_j as a parametric function, which is optimized jointly with the model parameters $\theta$ . In contrast, in statistics a fixed set of $p>m$ features $\left\{x_j=f_j\left(z\right)\right\}$ is considered, from which a subset of size m is selected. Although the former is computationally more tractable, it often leads to nonconvex optimization problems. The focus of this study lies in the latter. In particular, it is assumed that the feature space is very large and that only a small subset of size m must be selected ( $p\gg m$ ). Ideally, the parameters of a statistical model are estimated for a given loss function $\omega :{ℛ}^p\to ℛ$ by minimizing $\omega \left(\theta \right)$ subject to the constraint ${∥\theta ∥}_0=m$ . However, this approach is generally not feasible since it requires to test all subsets of size m.

Several strategies exist to compute approximate solutions. For instance, matching pursuit (MP) (Tropp et al., 2007) is a greedy strategy that selects one feature at a time until the desired number of m features is reached. This algorithm was found to be overly greedy in practice (Efron et al., 2004) and may provide poor subset selections. A more powerful strategy is to consider the $l_1$ constraint ${∥\theta ∥}_1=\Lambda$ (Tibshirani, 1996). In practice, this leads to the equivalent unconstrained optimization problem $$\widehat{\theta }={argmin}_{\theta }\omega \left(\theta \right)+\lambda {∥\theta ∥}_1,$$(P0)

where the penalty strength $\lambda \ge 0$ is a parameter that controls the sparseness of $\theta$ and is often chosen so that ${∥\theta ∥}_0=m$ . Although this approach often selects nearly optimal subsets, it comes with the difficulty of determining the penalty or regularization strength $\lambda$ .

Unfortunately, there exists no known functional relation between m and $\lambda$ . However, grid search can be used to compute solutions of the optimization problem P0 for a manually specified set of regularization strengths. More efficient algorithms, such as least angle regression (LARS) (Efron et al., 2004) and the homotopy (Osborne et al., 2000) algorithm, initially start with $\lambda ={\lambda }_{max}$ for which all coefficients are zero. These methods decrease $\lambda$ and determine breakpoints ${\lambda }_1>{\lambda }_2>\dots {\lambda }_q$ at which features either become active (nonzero coefficient) or inactive (zero coefficient). By interpolating between breakpoints a regularization path $\left\{\widehat{\theta }\left(\lambda \right)|0\le \lambda \le {\lambda }_{max}\right\}$ is obtained, which contains solutions for all feasible regularization strengths.

Computing full regularization paths is costly for large feature spaces and complex models. In such cases, the regularization path is only computed up to a maximal $m_q\ll p$ . Furthermore, a solution of the optimization problem P0 with m features will often behave very similarly to a solution with $m+1$ features. Therefore, the focus of this study is to present and test an algorithm that computes leapfrog regularization paths containing solutions $\widehat{\theta }\left({\lambda }_k\right)$ such that ${∥\widehat{\theta }\left({\lambda }_k\right)∥}_0=m_k$ , where $m_1<m_2<\dots <m_q$ is a given sequence of cardinalities. The main advantage of this algorithm is that m_q can be chosen much larger than with existing regularization path algorithms, because it does not compute all solutions of the regularization path. Furthermore, the sequence of cardinalities $m_1<m_2<\dots <m_q$ can be selected so that a clear difference in performance for m_k and $m_{k+1}$ is observed.

2. Algorithmic Background

Many feature selection methods were first developed for linear regression and later extended to more complex models, such as generalized linear models (GLMs). Hence, linear regression may serve as a role model for introducing feature selection methods. Let $X\in {ℛ}^{n\times p}$ denote the covariates or data matrix with n observations $\left\{x_i\right\}\left(rows\right)$ and p features $\left\{f_j\right\}\left(columns\right)$ . For linear regression, $\sigma$ is the identity function so that ${\mu }_{\theta }\left(x_i\right)=x_i\theta$ and the objective function $\omega ={\omega }_L$ is the sum of squared residuals defined as

Furthermore, $y\in {ℛ}^n$ is a vector of n dependent variables, also called response, and $\epsilon =y-X\theta$ is referred to as residuals. If $p\le n$ and X has full rank p then the coefficients $\theta$ can be estimated by OLS, which minimizes the sum of squared residuals ${\omega }_L\left(\theta \right)={∥\epsilon ∥}_2^2={∥y-X\theta ∥}_2^2$ . When $p>n$ , parameters can be estimated by minimizing the $l_1$ -penalized loss $L_{\lambda }\left(\theta \right)={\omega }_L\left(\theta \right)+\lambda {∥\theta ∥}_1$ , that is $$\widehat{\theta }={arg\; min}_{\theta }{L}_{\lambda }\left(\theta \right),$$(P1)

which is known as the Lasso (Tibshirani, 1996).

The following interpretation of OLS (Hastie et al., 2009) is essential for understanding many feature selection methods. In ${ℛ}^n$ the dot product $X\theta$ is a point in the hyperplane H spanned by the columns $\left[f_1,\dots ,f_p\right]$ of X and y is a vector typically not contained in H. The point $X\widehat{\theta }$ in H that minimizes the length of the residual vector $\epsilon$ is located directly below y, that is, it is the projection of y onto H. It follows that

where is a projection matrix. Therefore, the OLS solution can be interpreted as the scalar projection of y onto H. If the covariates are appropriately normalized, that is, Most importantly, X^Ty can be interpreted as the correlation of y with X. The OLS solution satisfies , which shows that the residuals are uncorrelated with X. Many feature selection methods decrease the correlation iteratively starting with $c\left(0\right)$ until the residuals are orthogonal to X.

2.1. Matching pursuit

MP is a greedy feature selection algorithm (Tropp et al., 2007) that selects one feature at a time. Let the columns of X have unit length, that is, for all $j=1,\dots ,p$ . The first feature $j_1\in \left\{1,\dots p\right\}$ is computed as . Hence, the algorithm selects the feature j with the largest scalar projection of y onto f_j, or the feature that is most correlated with y.

All remaining features are selected accordingly. Let $\Omega =\left\{j_1,\dots ,j_t\right\}$ denote the set of active features and the covariate matrix restricted to $\Omega$ . Furthermore, let

denote the residuals at iteration t with ${\widehat{\theta }}_{\Omega }={argmin}_{\theta }{∥y-X_{\Omega }\theta ∥}_2^2$ . At iteration $t+1$ the algorithm computes $$j_{t+1}={arg\; min}_j{∥{\epsilon }_t-f_j{\widehat{\theta }}_j∥}_2^2={arg\; max}_{j\notin \Omega }\left|f_j^T{\epsilon }_t\right|,$$

that is, it selects that feature j for which the scalar projection of ${\epsilon }_t$ onto f_j is largest. Alternatively, can be interpreted as the correlation of the jth feature with the residuals ${\epsilon }_t$ .

The MP algorithm has been applied to several other models, including logistic regression (Lozano et al., 2011), although the correlation might be less informative.

2.2. LARS and homotopy algorithm

LARS is an iterative feature selection method for linear regression (Efron et al., 2004). The algorithm is related to MP, but it is less greedy and gives only as much weight to each feature as it deserves (Hastie et al., 2009). Another method is the homotopy algorithm (Osborne et al., 2000), which is identical to LARS except that it also allows features to be removed within an iteration. It is referred to as LARS with Lasso modification in the literature (Efron et al., 2004; Donoho and Tsaig, 2008). Both methods are introduced hereunder. X is typically assumed to be standardized so that the correlations $c\left(\theta \right)$ are not tainted by heterogeneous scaling.

The LARS and homotopy algorithms maintain a set of active features $\Omega \subset \left\{1,\dots ,p\right\}$ all equally correlated with the residuals $y-X\widehat{\theta }$ for the current estimate $\widehat{\theta }$ . In contrast, features not in $\Omega$ are less correlated with the residuals. Both algorithms initially set $\theta =0$ . In each subsequent iteration, the coefficients $\theta$ are updated so that the correlations of all features in $\Omega$ uniformly shrink toward zero until some other feature $j\prime \in {\Omega }^c$ is equally correlated with the residuals. For the next iteration, $j\prime$ is added to the active set $\Omega$ . The homotopy algorithm also removes a feature j from the active set $\Omega$ when the coefficient ${\theta }_j$ becomes zero. In this way, $\lambda$ is incrementally reduced from its maximum value, where all features are inactive ( $\Omega =\oslash$ ), to zero, where all coefficients are nonzero ( $\Omega =\left\{1,2,\dots ,p\right\}$ ). The values of $\lambda$ at which a feature enters or leaves the active set $\Omega$ are referred to as breakpoints.

More specifically, the decrement in correlation at which a feature enters the active set is given by

where min⁺ is the minimum over positive elements and $c_j\left(\theta \right)$ denotes the jth element of $c\left(\theta \right).$ For LARS without Lasso modification, the step size ${\gamma }^{\ast }$ equals ${\gamma }^{+}$ ; however, the homotopy algorithm also removes a feature j from the active set when for some $\gamma$ $${\theta }_j+\gamma v_j=0,$$

so that ${\gamma }^{-}={min}_{j\in \Omega }\left\{-{\theta }_j∕v_j\right\}$ . The subsequent breakpoint is given by ${\gamma }^{\ast }=min\left\{{\gamma }^{+},{\gamma }^{-}\right\}$ .

LARS with Lasso modification is equivalent to solving the Lasso problem (P1) for all regularization strengths $\lambda \ge 0$ . A direct application of LARS or the homotopy algorithm to GLMs is not possible, because the path along which the correlations $c_{\Omega }\left(\theta \right)=\lambda 1$ uniformly decline is nonlinear. As a consequence, there is no analytical solution to the position of breakpoints. Nevertheless, a linear approximation for GLMs exists (Park and Hastie, 2007). The algorithm, however, requires several iterations until a breakpoint is reached.

3. Computing Leapfrog Regularization Paths

LARS and the homotopy algorithm compute full regularization paths $\left\{\widehat{\theta }\left(\lambda \right)|0\le \lambda \le {\lambda }_{max}\right\}$ for all feasible regularization strengths $\lambda$ . These algorithms start with $\lambda ={\lambda }_{max}$ and decrease the regularization strength until the unconstrained model is obtained. The full and exact regularization path is computed by linear interpolation between breakpoints. Extensions of LARS to GLMs (LARS-GLM) instead compute approximate regularization paths, because $\widehat{\theta }\left(\lambda \right)$ is nonlinear between breakpoints and there exists no analytical expression to compute at what values of $\lambda$ breakpoints occur.

For certain applications the full regularization path can be very costly to compute, especially when the number of features p is large. In such cases, the algorithms can be stopped when a maximum number of features m_q is reached. Nevertheless, because LARS-GLM successively decreases the regularization strength $\lambda$ , it might take very long until a desired number of features m is reached. Furthermore, estimated models with m and $m+1$ selected features might behave very similarly, for instance, in terms of classification performance. For such situations, an algorithm is introduced that computes for an a priori defined sequence $m_1<m_2<\dots <m_q$ the estimates $\widehat{\theta }\left({\lambda }_k\right)$ such that ${∥\widehat{\theta }\left({\lambda }_k\right)∥}_0=m_k$ . The set of estimates ${\left\{\widehat{\theta }\left({\lambda }_k\right)|{∥\widehat{\theta }\left({\lambda }_k\right)∥}_0=m_k\right\}}_{k=1}^q$ is called a leapfrog regularization path. The algorithm can be derived from proximal gradient descent (Boyd et al., 2004) or from a modification of LARS-GLM.

Let $\omega$ denote the objective function and ${\omega }_j^{\prime }\left(\widehat{\theta }\right)=\frac{\partial }{\partial {\theta }_j}\omega \left(\theta \right){|}_{\theta =\widehat{\theta }}$ the jth partial derivative at the current estimate $\widehat{\theta }$ . Given a fixed $m\in \left\{m_1,m_2,\dots ,m_q\right\}$ the following update rules are iterated until convergence:

The algorithm has converged when $\Omega$ is identical between two successive iterations. The first update step computes a new active set $\Omega$ . All features with nonzero coefficients remain in the active set. Additional features are added based on the absolute value of the partial derivatives. An estimate of the optimal $\lambda$ parameter is obtained during this step. In the second step, the coefficients of the active set are re-estimated, all other coefficients remain zero. The re-estimation is computationally cheap if $m\ll p$ , because the gradient descent method operates only on m selected features. Once the algorithm has converged for a fixed $m=m_k$ , the procedure is repeated for $m=m_{k+1}$ (see Algorithm 1).

Algorithm 1 Calculate leapfrog regularization path ${\left\{\widehat{\theta }\left({\lambda }_k\right)|{∥\widehat{\theta }\left({\lambda }_k\right)∥}_0=m_k\right\}}_{k=1}^q$
1: procedure LeapfrogPath $\left(m_1<m_2<\dots <m_q\right)$
2:  $r\leftarrow \oslash$
3:  $\widehat{\theta }\leftarrow 0\in {ℛ}^p$
4: for $k=1,\dots ,q$ do
5:  repeat
${\widehat{\theta }}_{\Omega }\leftarrow {arg\; min}_{\theta \in {ℛ}^{\left|\Omega \right|}}{\omega }_{\Omega }\left(\theta \right)+\lambda {∥\theta ∥}_1$
6:  until $\Omega$ remains unchanged
7:   $r\leftarrow r\cup \left\{{\widehat{\theta }}_{\Omega }\right\}$
8: end for
9: return r
10: end procedure

Algorithm 1 relies on $l_1$ -regularization for feature selection. Compared with MP it is much less greedy, because it does not simply select one feature per iteration (Hastie et al., 2009). Instead, the active set is refined until a solution to optimization problem (P0) is found. The algorithm is more similar to LARS-GLM, but it uses a different method for determining the regularization strength $\lambda$ . Whereas LARS estimates the position of subsequent breakpoints, Algorithm 1 iteratively updates $\lambda$ until a given number of m features is selected.

Despite extremely high-dimensional feature space, the estimation of parameters is relatively cheap when $m\ll p$ , because the gradient descent method operates only on a small subset of m features and not the full feature space of size p. In other words, estimating the coefficients of the logistic regression with a gradient method takes $O\left(nmT\right)$ steps, where T is the number of iterations of the gradient descent method. This is much faster than estimating the parameters on the full feature space, which would take $O\left(npT\right)$ steps.

3.1. Derivation from proximal gradient descent

The optimization problem P0 is commonly solved using proximal gradient descent (Boyd et al., 2004). In general, proximal gradient descent can be used to compute minimizers of functions $f\left(x\right)=g\left(x\right)+h\left(x\right)$ , where g is convex differentiable and h is convex but not differentiable. The update equation is given by $$x\leftarrow pro{x}_h\left(x-\gamma \nabla g\left(x\right)\right),$$

$$wherepro{x}_h\left(x\right)={arg\; min}_{x\prime }\frac{1}{2\gamma }{∥x\prime -x∥}_2^2+h\left(x\prime \right)$$

and $\gamma$ is the step size. The proximal operator $pro{x}_h$ can be solved analytically for $h\left(x\right)={∥x∥}_1$ . Note that the equation inside the proximal operator is a simple gradient descent step applied to $g.$ .

For the optimization problem $P0$ the following update equations are obtained: $$\vartheta \leftarrow \theta -\gamma \nabla \omega \left(\theta \right)$$ $$\theta \leftarrow {arg\; min}_{\theta \prime }\frac{1}{2\gamma }{∥\theta \prime -\vartheta ∥}_2^2+{∥\theta \prime ∥}_1=sign\left(\vartheta \right)max\left\{\left|\vartheta \right|-\gamma \lambda ,0\right\},$$

where $\vartheta$ is the result after the gradient descent step and before applying the proximal operator. Assume that for a feature j the coefficient ${\theta }_j$ is zero. The coefficient remains zero during subsequent update steps unless the gradient update is large enough, because of the application of the proximal operator. More specifically, for any step size $\gamma >0$ , a coefficient ${\theta }_j$ remains zero unless $\left|{\omega }_j^{\prime }\left(\widehat{\theta }\right)\right|>\lambda$ .

3.2. Derivation from LARS

In each iteration, LARS either removes a feature from the active set or selects a new one. The focus here is on the selection of new features, that is, Equation (1). By assuming that $f_j^{T}Xv=0$ , a much simplified step size

is obtained, which shows that $\lambda -{\gamma }^{+}={max}_{j\in {\Omega }^c}\left|c_j\left(\theta \right)\right|$ and remember that for linear and logistic regression $c_j\left(\theta \right)=\omega {\prime }_j\left(\theta \right)$ . Hence, m features can be selected by setting $\lambda$ equal to the mth largest absolute correlation or partial derivative.

4. Applications to Large-Scale Logistic Regression

Algorithm 1 can be easily applied to GLMs, such as logistic regression, as demonstrated in the following. The closest method in the literature is LARS-GLM and the GNU-R library glmpath is used for comparing it with Algorithm 1. Decisive for the computational cost is the number of steps required for computing a regularization path. A step here refers to the estimation of parameters for a fixed regularization strength $\lambda$ , that is, solving optimization problem $P0$ . For Algorithm 1 this corresponds to one evaluation of the two update equations. In addition, the number of passes of the gradient descent method through the entire data set (epochs) is informative. However, because different gradient descent methods are used for Algorithm 1 and LARS-GLM, a direct comparison in terms of epochs is not possible.

Enhancers are regulatory elements in the genome that control the cell type-specific expression of genes. The DNA sequence of an enhancer is predictive of the cell types in which it actively controls gene expression (Lee et al., 2011; Hashimoto et al., 2016; Kelley et al., 2016; Yang et al., 2017). Training and test data sets were compiled from DNA sequences of $56,880$ enhancers active in embryonic mice at day 15.5 in brain (positive set) and nonbrain tissues (negative set). Enhancers were detected using ModHMM (Benner and Vingron, 2020) and the data are available online (Benner, 2020). For each enhancer sequence, the K-mer content is determined. K-mers are short DNA subsequences of length K. In this example, K-mers of variable length are considered. The resulting data matrix $X=\left(x_{ij}\right)$ has n rows (observations) and p columns (features), where $x_{ij}$ is the number of occurrences of the jth K-mer in the ith enhancer sequence. A logistic regression model is used to discriminate between brain and nonbrain enhancers. The average negative log-likelihood (average loss) is given by

where $\sigma$ is the sigmoid function and $y\in {\left\{0,1\right\}}^n$ is a vector of n class labels. Analogously to linear regression, when $p>n$ , parameters are estimated by minimizing the $l_1$ -penalized average logistic loss $G_{\lambda }\left(\theta \right)={\omega }_G\left(\theta \right)+\lambda ∥\theta ∥$ (penalized average loss), that is $$\widehat{\theta }={argmin}_{\theta }{G}_{\lambda }\left(\theta \right).$$

For logistic regression the objective function $G_{\lambda }$ is convex and the optimum $\widehat{\theta }$ must satisfy . SAGA (Defazio et al., 2014) is used as gradient descent method to compute solutions of this optimization problem. To guarantee stable convergence of the gradient descent method, the data are normalized to unit variance.

In a first example, K-mers of length $2-7$ are used, that is, $K\in \left\{2,3,\dots ,7\right\}$ . This leads to a feature space of size $p=21,840$ , where only K-mers that actually occur in the data set are considered. Figure 1 shows the leapfrog regularization path computed with Algorithm 1 for $m\in \left\{1,10,500\right\}$ . In addition, a regularization path computed with LARS-GLM is shown and the algorithm was stopped when it reached 500 features. The example shows that the computation of the leapfrog path requires much fewer steps, however, at the expense of not evaluating the full path up to 500 features. Nevertheless, for larger feature spaces and larger m the LARS-GLM algorithm quickly becomes computationally too expensive.

FIG. 1.

Comparison of regularization paths. Leapfrog regularization path (left) with $m\in \left\{1,10,500\right\}$ and LARS-GLM regularization path (right). Vertical dashed gray lines show steps of the path algorithms, that is, $\lambda$ values for which parameters $\widehat{\theta }$ were estimated

In a second example, the set of K-mers is extended to lengths of $2-12$ , that is, $K\in \left\{2,3,\dots ,12\right\}$ , which results in a feature space of size $p=16,876,344$ . Figure 2 shows the estimation of parameters for a leapfrog regularization path with Algorithm 1 and targets $m\in \left\{10,100,1000,5000\right\}$ . The algorithm successively decreases the penalty strength $\lambda$ until a target is reached. After only 1500 epochs, corresponding to ∼50 steps of Algorithm 1, the leapfrog regularization path is computed. In contrast, extracting the 5000 most important K-mers with the LARS-GLM algorithm is computationally not feasible.

FIG. 2.

Estimation of a leapfrog regularization path with $m\in \left\{10,100,1000,5000\right\}$ . The vertical lines show the epoch at which m is increased. An epoch is one pass of the gradient descent method through the entire data. For a fixed penalty strength $\lambda$ , the gradient descent method is iterated until convergence. If the number of nonzero parameters does not match the target m, $\lambda$ is updated and the parameters are re-estimated.

4.1. Bias-variance trade-off

In statistics and machine learning one is often confronted with the task of choosing a model family $ℍ$ that captures the underlying structure of the observed data, but also generalizes well to new observations. The bias-variance trade-off suggests that both cannot be accomplished with arbitrary precision. When the capacity of $ℍ$ is too small, the estimated models has poor performance on the observed data as well as on new observations. The underfitting of the estimated models is caused by a strong (inductive) bias of $ℍ$ . As the capacity of $ℍ$ increases, the model is capable of extracting more complex patterns from the training data, but the variance of parameter estimates across samples increases. If the capacity of $ℍ$ is too large it causes overfitting, that is, the estimated models will perform very well on the training data but generalize poorly. A model family that generalizes well must minimize both sources of error, the bias and the variance.

For instance, in the case of linear or logistic regression, the capacity of $ℍ$ can be increased by using more features (i.e., increasing p). Regularization, in contrast, constrains the parameter space and, therefore, reduces the capacity of $ℍ$ (Hastie et al., 2009). In practice, it is often observed that models generalize best if regularization is used in combination with very large p (Bühlmann and Van De Geer, 2011).

In this section, the effect of regularization and the size of the feature space p on the capacity of $ℍ$ is studied. A balanced set of n = 10,000 brain and nonbrain enhancers is used. The size of the feature space p is controlled by the maximal length of K-mers. Figure 3 shows the training and test error for a fixed penalty $\lambda$ and increasing p. Whereas the training error decreases monotonically, the test error starts to decrease after it reaches a maximum. This double descent has been observed before (Vallet et al., 1989; Opper et al., 1990; Duin, 2000; Belkin et al., 2019) and is theoretically expected (Cheema and Sugiyama, 2020). It shows that the capacity of $ℍ$ is the result of a complex interplay between the size of the feature space p and the regularization strength $\lambda$ .

FIG. 3.

Average loss as a function of p and fixed penalty $\lambda$ .

The dependency between $\lambda$ and p makes it difficult in practice to find optimal values for both parameters, especially when double descents are possible, because an extensive search is required to eliminate the possibility of being in a local optimum. Instead of adding more parameters by increasing p, Algorithm 1 is utilized to control the number of active features m for a fixed p. Figure 4 shows that in this case no double descent is observed, suggesting a monotonic relation between m and the capacity of $ℍ$ , as expected from the bias-variance trade-off. Furthermore, the test error increases rapidly as the interpolation threshold is crossed, that is, the point where the model is able to perfectly discriminate between the positive and negative class.

FIG. 4.

Average loss as a function of m and fixed p.

Figure 5 shows the converse case, the average loss for fixed m but variable p. Also in this case no double descent of the test error is observed. However, the training error increases with p after reaching a minimum, because the regularization strength must increase to keep m fixed. The generalization gap, that is, the difference between test and training loss, steadily decreases, showing that the capacity of $ℍ$ seems to decline.

FIG. 5.

Average loss as a function of p and fixed m.

5. Conclusion

Several iterative feature selection methods were discussed in this study. MP is among the simplest methods, but can lead to poor subset selection because of its overly greedy strategy (Hastie et al., 2009). More advanced methods, such as LARS and the homotopy algorithm, often yield better results in practice. The estimates of LARS with Lasso modification are equivalent to solving the Lasso problem (P1) for all regularization strengths $\lambda$ . Extensions to GLMs are more expensive to compute, because there exists no analytical solution to the $\lambda$ values at which the active set of features changes. As a result, many more iterations are required to compute a full regularization path. In practice, one is often not interested in the full solution, because estimates with m and $m+1$ features provide very similar predictions.

Algorithm 1 instead computes for an a priori defined sequence $m_1<m_2<\dots <m_q$ a leapfrog regularization path ${\left\{\widehat{\theta }\left({\lambda }_k\right)|{∥\widehat{\theta }\left({\lambda }_k\right)∥}_0=m_k\right\}}_{k=1}^q$ . This algorithm can be derived from proximal gradient descent or a simple modification of LARS for GLMs. It is very effective for high-dimensional problems where only a small subset of features is required for making accurate predictions. It is less greedy than MP, because it does not simply select one feature per iteration. As opposed to LARS, which computes subsequent breakpoints, Algorithm 1 iteratively refines the regularization strength $\lambda$ until a given number of m features is selected.

On a data set of DNA sequences the algorithm very efficiently computes the leapfrog regularization path despite the high-dimensional feature space, where LARS-GLM is computationally too expensive. In addition, the number of selected features m allows a good control of model complexity. An implementation of the algorithm is available at: https://github.com/pbenner/kmerLr