{'fraction_non_alphanumeric': 0.02812040086901675, 'fraction_numerical': 0.0061847361412853035, 'mean_word_length': 4.770599534930745, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
ADMM Algorithmic Regularization Paths for Sparse Statistical Machine Learning

Optimization approaches based on operator splitting are becoming popular for solving sparsity regularized statistical machine learning models. While many have proposed fast algorithms to solve these problems for a single regularization parameter, conspicuously less attention has been given to computing regularization paths, or solving the optimization problems over the full range of regularization parameters to obtain a sequence of sparse models. In this chapter, we aim to quickly approximate the sequence of sparse models associated with regularization paths for the purposes of statistical model selection by using the building blocks from a classical operator splitting method, the Alternating Direction Method of Multipliers (ADMM). We begin by proposing an ADMM algorithm that uses warm-starts to quickly compute the regularization path. Then, by employing approximations along this warm-starting ADMM algorithm, we propose a novel concept that we term the ADMM Algorithmic Regularization Path. Our method can quickly outline the sequence of sparse models associated with the regularization path in computational time that is often less than that of using the ADMM algorithm to solve the problem at a single regularization parameter. We demonstrate the applicability and substantial computational savings of our approach through three popular examples, sparse linear regression, reduced-rank multi-task learning, and convex clustering.


Introduction
With the rise of Big Data and the subsequent explosion of statistical machine learning methods to analyze it, statisticians have become avid consumers of large-scale optimization procedures to estimate sparse models. The estimation problem is often cast as an optimization problem of the form: where β β β is a parameter which specifies a statistical model, L(β β β ; W) is a smooth loss function or data-fidelity term that quantifies the discrepancy between the data, W, and the model specified by β β β , and P(β β β ) is a nonsmooth penalty that encourages sparsity in model parameter β β β [3,4,15]. A regularization parameter, λ ≥ 0, explicitly trades off the model fit and the model complexity.
Directly solving the optimization problem (1) is often challenging. Operator splitting methods, such as the Alternating Direction Method of Multipliers (ADMM), have become popular because they convert solving the problem into solving a sequence of simpler optimization problems that involve only the smooth loss or nonsmooth penalty. By breaking up the problem into smaller ones, ADMM may end up taking more iterations than directly solving (1), but it often runs in less total time since the subproblems are typically easy to solve. Clearly in the context of Big Data, faster algorithms are indispensable, and the numerical optimization community has devoted a great deal of effort to solving (1) rapidly for a fixed value of λ . This goal, however, is not necessarily aligned with the application of statistical machine learning problems to real data.
In practice, statisticians are interested in finding the best sparse model that represents the data. Achieving this typically entails a two-step procedure: (i) model selection, or selecting the best sparse model or equivalently the best subset of parameters, and (ii) model fitting, or fitting the model by minimizing the loss function over the selected parameters [15]. The first step is often the most challenging computationally as this entails searching the combinatorial space of all possible sparse models. As this combinatorial search is infeasible for large-scale problems, many consider convex relaxations through constraints or penalties as computationally feasible surrogates to help search through the space of sparse models. Consider for example, sparse linear regression, where the goal is to find the subset of variables or inputs that best predicts the response or output. Searching over all possible subsets of variables, however, is an NP hard problem. Instead, many have employed the penalty or constraint, P(β β β ) = β β β 1 , which is the tightest convex relaxation to performing best subset selection and whose solution can be computed in polynomial time. The nonsmooth penalty term, P(β β β ), then serves to translate an infeasible computational problem into a tractable one for model selection purposes.
Suppose now that we focus on selecting the best sparse model by means of penalized statistical estimation as in (1). As λ varies, we trace out a continuous parametric curveβ β β (λ ) ∈ R p . Since this curve cannot be determined analytically in general, the curve is estimated for a finite sequence of regularization parameters. To choose the best model, statisticians inspect the sequence of sparse solutions to (1) over the full range of regularization parameters: {β β β (λ n ) : 0 ≤ λ 1 ≤ · · · ≤ λ max }, where λ max is the value of λ at whichβ β β (λ max ) = 0, the maximally sparse solution. This sequence of sparse solutions is often called the regularization path [9,10,14]. For model selection purposes, however, the actual parameter values,β β β (λ ), as λ varies in the regularization paths are less important than identifying the non-zero components ofβ β β (λ ). (Note that the parameter values for the optimal model are typically re-fit anyways in the second model fitting stage.) Instead, the support ofβ β β (λ ) or the sequence of active sets defined as A (λ ) = { j :β j (λ ) = 0}, are the important items; these yield a good sequence of sparse models to consider that limit computationally intensive exploration of a combinatorial model space. Out of this regularization path or sequence of active sets, the optimal model can be chosen via a number of popular techniques such as minimizing the trade-off in model complexity as with the AIC and BIC, the prediction error as with cross-validation [15] or the model stability as with stability selection [23].
To apply sparse statistical learning methods to large-scale problems, we need fast algorithms not only to fit (1) for one value of λ , but to estimate the entire sequence of sparse models in the model selection stage. Our objective in this chapter is to study the latter, which has received relatively little attention from the optimization community. Specifically, we seek to develop a new method to approximate the sequence of active sets associated with regularization paths that is (i) computationally fast and (ii) comprehensively explores the space of sparse models at a sufficiently fine resolution. In doing so, we will not try to closely approximate the parameter values,β β β (λ ), but instead try to closely approximate the sparsity of the parameters, A (λ ), for the statistical learning problem (1).
To rapidly approximate the sequence of active sets associated with regularization paths, we turn to the ADMM optimization framework. We first introduce a procedure to estimate the regularization path by using the ADMM algorithm with warm starts over a range of regularization parameters to yield a path-like sequence of solutions. Extending this, we preform a one-step approximation along each point on this path, yielding the novel method that we term ADMM Algorithmic Regularization Paths. Our procedure can closely approximate active sets given by regularization paths at a fine resolution, but dramatically reduces computational time. This new approach to estimating a sequence of sparse models opens many interesting questions from both statistical and optimization perspectives. In this chapter, however, we focus on motivating our approach and demonstrating its computational advantages on several sparse statistical machine learning examples.
This chapter is organized as follows. We first review how ADMM algorithms have been used in the statistical machine learning literature, Section 1.1. Then, to motivate our approach, we consider application of ADMM to the familiar example of sparse linear regression, Section 1.2. In Section 2, we introduce our novel Algorithmic Regularization Paths for general sparse statistical machine learning procedures. We then demonstrate how to apply our methods through some popu-lar machine learning problems in Section 3; specifically, we consider three examples -sparse linear regression (Section 3.1), reduced-rank multi-task learning (Section 3.2), and convex clustering (Section 3.3) -where our Algorithm Paths yield substantial computational benefits. We conclude with a discussion of our work and the many open questions it raises in Section 4.

ADMM in Statistical Machine Learning
The ADMM algorithm has become popular in statistical machine learning in recent years because the resulting algorithms are typically simple to code and can scale efficiently to large problems. Although ADMM has been successfully applied over a diverse spectrum of problems, there are essentially two thematic challenges among the problems that ADMM has proven adept at addressing: (i) decoupling constraints and regularizers, that are straightforward to handle individually, but not in conjunction; and (ii) simplifying fusion type penalties. We make note of these two types of problems because the ADMM Algorithmic Regularization Path we introduce in this chapter can be applied to either type of problem.
An illustrative example of the first thematic challenge arises in sparse principal component analysis (PCA). In [35] Vu et al. propose estimating sparse principal subspace estimatorB of a symmetric input matrix S with the solution to the following semidefinite program: where B 1 is 1-norm of the vectorization of B, the set F d = {B : 0 B I, tr(B) = d} is a closed and convex set called the Fantope, and λ ≥ 0 is a regularization parameter. The main algorithmic challenge is the interaction between the Fantope constraint and the 1 -norm penalty. If only either the penalty or constraint were present the problem would be straightforward to solve. Consider the following equivalent problem to which ADMM can be readily applied: where δ C (Σ Σ Σ ) denotes the indicator function of the closed convex set C, namely the function that is 0 on C and ∞ otherwise. By minimizing an augmented Lagrangian over B, the copy variable Z, and the scaled dual variable U as outlined in [3], we arrive at the following ADMM updates: Thus, the penalty and constraint are effectively decoupled resulting in simple updates: the update for Z requires the soft-thresholding operator, S µ (x) = sign(x)(|x|− µ) + , and the update for B involves the projection onto the Fantope, denoted by P F d , which has a closed form solution given in [35]. The literature abounds with many more examples of using the ADMM splitting strategy to decouple an otherwise challenging optimization problem into simpler subproblems. Boyd et al. [3] review many such applications. Other example applications include decoupling trace or nuclear norm penalties as in robust PCA [42], latent variable graphical models [21], and tensor completion [20]; decoupling different types of hierarchical constraints [2], decoupling a series of loss functions [18], decoupling joint graphical models [6], and decoupling large linear programming problems [1], among many others.
The second thematic challenge that ADMM algorithms have been used to solve involve fusion or non-separable penalties. A good illustrative example of this challenge arises in total variation (TV) denoising [31]. Consider the simple version of this problem, specifically finding a smooth estimate of a noisy one-dimensional signal y ∈ ℜ n : where the tuning parameter λ ≥ 0 trades off the smoothness of the approximation with the goodness of fit with the data y. What makes this problem challenging is that the fusion penalty couples the non-smooth terms so that they are non-separable. Note that this penalty can be written more compactly as A x 1 where A is the discrete first order differences operator matrix. More generally, this second class of problems consist of problems of the form, L(β β β ; W) + λ P(A β β β ). In the machine learning context these penalties arise because we often wish to impose structure, not on a latent variable of interest directly, but rather on a linear transformation of it. In the TV denoising example we seek sparsity in differences of adjacent time points of the latent signal. Previously, we could break the objective into a sum of simpler objectives. The issue here is different; specifically the composition of the regularizer with a linear mapping complicates matters. ADMM can again greatly simplify this problem if we let the ADMM copy variable copy the linearly transformed parameters: The ADMM subproblems for iteratively solving this problem then have the following simple form: Note that we have eliminated having to minimize any functions containing the troublesome composition penalty. In the context of the TV denoising example, the β β β update requires solving a linear system of equations, and the z update involves a straightforward soft-threshold.
The ADMM algorithm has been used to decouple fusion or non-separable types of penalties in many statistical learning problems. These include more general instances of total variation [36,12], a convex formulation of clustering [5], joint graphical model selection [24,25], overlapping group lasso penalties [40], and more generally for structured sparsity [22].
Overall, while the ADMM algorithm is gaining more widespread application in statistics and machine learning, the algorithm is applied in the traditional sense to solve a composite optimization problem for one value of the regularization parameter. In this chapter, we seek to investigate the ADMM algorithm for a different purpose, namely to find a sequence of sparse models associated with regularization paths.

Developing an Algorithmic Regularization Path: Sparse Regression
Our goal is to quickly generate a sequence of candidate sparse solutions for model selection purposes. To achieve this, we will propose a method of approximating the sequence of active sets given by regularization paths, or the path-like sequence solutions of penalized statistical models over the full range of regularization parameters. To motivate our approach, we study the familiar example of sparse linear regression. Suppose we observe a covariate matrix X ∈ ℜ n×p consisting of n iid observations of p variables and an associated response variable y ∈ ℜ n . We are interested in fitting the linear model y = X β β β +ε where ε is independent noise, but assume that the linear coefficient vector β β β is sparse, β β β 0 p where · 0 is the 0 norm or the number of non-zero elements of β β β . Minimizing a criterion subject to a constraint of the form β β β 0 ≤ k for some k, becomes a combinatorially hard task. To estimate a sparse model in reasonable time, many have proposed to use the tightest convex relaxation, the 1 -norm penalty, commonly called the LASSO [33] in the statistical literature: where λ ≥ 0 is the regularization parameter controlling the sparsity of β β β . The full regularization path of solutions for the LASSO is the set of regression coefficients {β β β (λ ) : ∀ 0 ≤ λ ≤ λ max } where λ max = 1 n X T y ∞ is the smallest amount of regularization that yields the sparse solutionβ β β = 0. The regularization paths for the LASSO have been well-studied and in particular, are continuous and piece-wise linear [28,8,30]. These paths also outline a sequence of active sets or sparse models that smoothly increase in sparsity levels as λ decreases from the fully sparse solution at λ = λ max to the fully dense solution at λ = 0. Hence for model selection, one can limit exploration of the combinatorial space of sparse models to that of the sequence of active sets outlined by the LASSO regularization paths.
Computing the full regularization paths, however, can be a computational challenge. Several path following algorithms for the LASSO [28,30] and closely related algorithms such as such as Least Angle Regression (LAR) [8] and Orthogonal Matching Pursuit (OMP) [7] have been proposed; their computational complexity, however, is O(p 3 ) which is prohibitive for large-scale problems. Consequently, many have suggested to closely approximate these paths by solving a series of optimization problems over a grid of regularization parameter values. Specifically, this is typically done for a sequence of 100 log-spaced values from λ max to λ 1 = 0. Statisticians often employ homotopy, or warm-starts, to speed computation along the regularization path [9]. Warm-starts use the solution from the previous value of λ j ,β β β (λ j ), as the initialization for the optimization algorithm to solve the problem at λ j+1 . As the coefficients, β β β , change continuously in λ , warm-starts can dramatically reduce the number of iterations needed for convergence asβ β β (λ j ) is expected to be close toβ β β (λ j+1 ) for small changes from λ j to λ j+1 . Many consider shooting methods, or coordinate descent procedures [9,37], that use warm-starts and iterate over the active set for 100 log-spaced values of λ [10] to be the fastest approximate solvers of the LASSO regularization path.
We seek to further speed the computation of the sequence of active sets given by the regularization path by using a single path approximating algorithm instead of solving separate optimization problems over a grid of regularization parameter values. Our approach is motivated by two separate observations: (i) the evolution of the sparsity level of the iterates of the ADMM algorithm used to fit (2) for one value of λ , and (ii) the behavior of a new version of the ADMM algorithm that incorporates warm-starts to expedite computation of regularization paths. We study each of these motivations separately, beginning with the first.
Consider using ADMM to solve the LASSO problem. First, we split the differentiable loss function from the non-differentiable penalty term by introducing a copy z of the variable β β β in the penalty function, and adding an equality constraint forcing them to be equal. The LASSO problem (2) can then be re-written as: with its associated augmented Lagrangian: Here, u is the scaled dual variable of the same dimension as β β β and ρ is the algorithm tuning parameter. The ADMM algorithm then follows three steps (subproblems) to solve the LASSO: The benefit of solving this reformulation is simpler iterative updates. These three steps are iterated until convergence, typically measured by the primal and dual residuals [3]. The β β β -subproblem solves a linear regression with an additional quadratic ridge penalty. Solving the z-subproblem introduces sparsity. Notice that this is the proximal operator of the 1 -norm applied to β β β k − u k which is solved analytically via soft-thresholding.Finally, the dual update ensures that β β β is squeezed towards z and primal feasibility as the algorithm progresses. Consider the sparsity of the z iterates, z k 0 , for the LASSO problem. Notice that as the algorithm proceeds, z k becomes increasingly sparse; this is illustrated for a small simulated example in the left panel of Figure 1. Let us study why this occurs and its implications. Regardless of λ , the ADMM algorithm begins with a fully dense β β β 1 as this is the solution to a ridge problem with parameter ρ. Softthresholding in the z-subproblem then sets coefficients of small magnitude to zero. The first dual update, u 1 , has magnitude at most |λ |, meaning that the second β β β 2 update is essentially shrunken towards β β β 1 . Smaller coefficients decrease further in magnitude and soft-thresholding in the z-subproblem sets even more coefficients to zero. The algorithm thus proceeds until the sparsity of the z k stabilizes to that of the solution,β β β (λ ). Hence, the support of the z k has approximated the active set of the solution long before the iterates of the β β β -subproblem; the latter typically does not reach the sparsity of the solution until convergence when primal feasibility is achieved. While Figure 1 only illustrates that z k quickly converges to the correct sparsity level, we have observed empirically in all our examples that the active set outlined by z k 0 also quickly identifies the true non-zero elements of the solution, β β β (λ ).
Interestingly then, the z k quickly explore a sequence of sparse models going from dense to sparse, similar in nature to the sequence of sparse models outlined by the regularization path. While from Figure 1 we can see that this sequence of sparse models is not desirable as it does not smoothly transition in sparsity and does not fully explore the sparse model space, we nonetheless learn two important items from this: (i) We are motivated to consider using the algorithm iterates of the z-subproblem, as a possible means of quickly exploring the sparse model space; and (ii) we are motivated to consider a sequence of solutions going from dense to sparse as this naturally aligns with the sparsity levels observed in the ADMM algorithm iterates. Given these, we ask: Is it possible to use or modify the iterates of the ADMM algorithm to achieve a path-like smooth transition in sparsity levels similar in nature to the sparsity levels and active sets corresponding to regularization paths? One possible solution would be to employ warm-starts in the ADMM algorithm along a grid of regularization parameters similar to other popular algorithms for approximating regularization paths. Recall that warm-starts use the solution obtained at the previous value of λ as initialization for the next value of λ . We first introduce this new extension of the ADMM algorithm for approximating regularization paths in Algorithm 1 and then return to our motivation of studying the sequence of active sets outlined by this algorithm.
Our ADMM algorithm with warm starts is an alternative algorithm for fitting regularization paths. It begins with λ small corresponding to a dense model, fits the ADMM algorithm to obtain the solution, and then uses the previous solution, β β β (λ j−1 ), and dual variable, u(λ j−1 ), to initialize the ADMM algorithm for λ j .
Before considering the sequence of active sets outlined by this algorithm, we pause to discuss some noteworthy features. First, notice that the ADMM tuning parameter, ρ, does not appear in this algorithm. We have omitted this as a parameter by fixing ρ = 1 throughout the algorithm. Fixing ρ stands in contrast with the burgeoning literature on how to dynamically update ρ for ADMM algorithms [3]. For example, adaptive procedures that change ρ to speed up convergence are introduced in [16]. Others have proposed accelerated versions of the ADMM algorithm that achieve a similar phenomenon [11]. Changing the algorithm tuning parameter, however, is not conducive to achieving a path-like algorithm using warm-starts. Consider the z-subproblem which is solving by soft-thresholding at the level λ j /ρ. Thus, if ρ is changed in the algorithm, the sparsity levels of z dramatically change, Algorithm 1 ADMM with Warm Starts: Sparse Regression eliminating the advantages of using warm-starts to achieve smooth transitions in sparsity levels. Second, notice that we have switched the order of the sub-problems by beginning with the z-subproblem. While technically the order of the subproblems does not matter [38], we begin with the z-subproblem as this is where the sparsity is achieved through soft-thresholding at the value, λ ; hence, the solution for z is what changes when λ is increased.
Next, notice that our regularization paths go from dense to sparse, or λ small to large, which is the opposite of other path-wise algorithms and algorithms that approximate regularization paths over a grid of λ values [10]. Recall that our objective is to obtain a smooth path-like transition in sparsity levels corresponding to a sequence of active sets that fully explores the space of sparse models. Our warm-start procedure naturally aligns with the sparsity levels of the iterates of the ADMM algorithm which go from dense to sparse, thus ensuring a smooth transition in the sparsity level of z as λ is increased. Our warm-start procedure could certainly be employed going in the reverse direction from sparse to dense, but we have observed that this introduces discontinuities in the z k iterates and consequently their active sets as well, thus requiring more iterations for convergence. This behavior occurs as the solution of the β β β -subproblem is always more dense than that of the z-subproblem, even when employing warm-starts. Now, let us return to our motivation and consider the sparsity levels and corresponding sequence of active sets achieved by the iterates of our new pathapproximating ADMM Algorithm. The sparsity of the iterates of the z-subproblems, z k 0 , are plotted for 30 log-spaced values of λ for the same simulated example in the middle panel of Figure 1. The iterates over all values of λ are plotted on the xaxis with vertical lines denoting the increase to the next λ value. Carefully consider the sparsity levels of the z iterates for each fixed value of λ in our ADMM algorithm with warm starts. Notice that the sparsity levels of z typically stabilize to that of the solution within the first few iterations after λ is increased. The remaining iterations and a large proportion of the computational time are spent on squeezing β β β towards z to satisfy primal feasibility. This means that the z-subproblem has achieved the sparsity associated with the active set ofβ β β (λ ) within a few iterations of increasing λ . One could surmise that if the increase in λ were small enough, then the z-subproblem could correctly approximate the active set corresponding to λ within one iteration when using this warm-start procedure. The right panel of Figure 1 illustrates the sparsity levels achieved by the z-subproblem for this sequence of one-step approximations to our ADMM algorithm with warm-starts. Notice that this procedure achieves a smooth transition in sparsity levels corresponding to a sequence of active sets that fully explore the range of possible sparse models, but requires only a fraction of the total number of iterations and compute time. This, then is the motivation for our new ADMM Algorithmic Regularization Paths introduced in the next section.

The Algorithmic Regularization Path
Our objective is to use the ADMM splitting method as the foundation upon which to develop a new approximation to the sequence of sparse solutions outlined by regularization paths. In doing so, we are not interested in estimating parameter values by solving a statistical learning optimization problem with high precision. Instead, we are interested in quickly exploring the space of sparse model at a fine resolution for model selection purposes by approximating the sequence of active sets given by the regularization path.
Again, consider the general sparse statistical machine learning problem of the following form: minimize where W denotes the "data" (for the sparse linear regression example, W = {X, y}), the loss function, L(β β β ; W) is a differentiable, convex function of β β β , and the regularization term, P : ℜ p → ℜ + is a convex and non-differentiable penalty function. As before, λ ≥ 0 is the regularization parameter controlling the trade-off between the penalty and loss function. Following the setup of the ADMM algorithm, consider splitting the smooth loss from the nonsmooth penalty through the copy variable, z: With scaled dual variable u, the augmented Lagrangian of general problem (3) is Now following from the motivations discussed in the previous section, there are three key ingredients that we employ in our Algorithmic Regularization Paths: (i) warm-starts to go from a dense to a sparse solution, (ii) the sparsity patterns of the z-subproblem iterates, and (iii) one-step approximations at each regularization level. We put these ingredients together in Algorithm 2 to give our Algorithmic Regularization Paths: Algorithm 2 Algorithmic Regularization Path for Sparse Statistical Learning Output {z k : k = 1, · · · , K} as algorithmic regularization path .
Our Algorithmic Regularization Path, Algorithm Path for short, outlines a sequence of sparse models going from fully dense to fully sparse. This can be used as an approximation to the sequence of active sets given by regularization paths for the purpose of model selection. Consider that the algorithm begins will the fully dense ridge solution. It then gradually increases the amount of regularization, γ, performing one full iterate of the ADMM algorithm (β β β -subproblem, z-subproblem, and dual update) for each new level of regularization. The regularization level is increased until the z-subproblem becomes fully sparse.
One may ask why we would expect our Algorithm Path to yield a sequence of active sets that well approximate those of the regularization path? While a mathematical proof of this is beyond the scope of this chapter, we outline the intuition stemming from our three key ingredients. (Note that we also demonstrate this through specific examples in the next section).
(i) Warm-starts from dense to sparse. Beginning with a dense solution and gradually increasing the amount of regularization ensures a smooth decrease in the sparsity levels corresponding to a smooth pruning of the active set as this naturally aligns with sparsity levels of the ADMM algorithm iterates. (ii) z-subproblem iterates. The iterates of the z-subproblem encode the sparsity of the active set,β β β (λ ), quickly as compared to the β β β -subproblem which achieves sparsity only in the limit upon algorithm convergence. (iii) One-step approximations. For a small increase in regularization when using warm-starts, the iterates of the z-subproblem often achieve the sparsity level of the active set within one-step.
Notice that if we iterated the three subproblems of our Algorithm Path fully until convergence, then our algorithm would be equivalent to our ADMM Algorithm with warm starts; thus, the one-step approximation is the major difference between Algorithms 1 and 2. Because of this one-step approximation, we are not fully solving (1) and thus the parameter values, β β β , will never stabilize to that of the regularization path. Instead, our Algorithm Path quickly approximates the sequence of active sets corresponding to the regularization path, as encoded in the z-subproblem iterates.
The astute reader will notice that we have denoted the regularization parameters in Algorithm 2 as γ instead of λ as in (1). This was intentional since due to the one-step approximation, we are not solving (1) and thus the level of regularization achieved, γ, does not correspond to λ from (1). Also notice that we have introduced a step size, t, that increases the regularization level, γ, at each iteration. The sequence of γ's can either be linearly spaced, as with additive t, or geometrically spaced, as with multiplicative t. Again, if t is very small, then we expect the sparsity patterns of our Algorithm Paths to well approximate the active sets of regularization paths.
We will explore the behavior and benefits of our Algorithm Paths through demonstrations on popular sparse statistical learning problems in the next section. Before presenting specific examples, however, we pause to outline three important advantages that are general to sparse statistical learning problems of the form (1).

Easy to implement. Our Algorithm Path is much simpler than other algorithms
to approximate regularization paths. The hardest parts, the β β β and z subproblems, often have analytical solutions for many popular statistical learning methods. Then, with only one loop, our algorithm can often be implemented in a few lines of code. This is in contrast to other algorithm paths which require multiple loops and much overhead to track algorithm convergence or the coordinates of active sets [10]. 2. Finer resolution exploration of the model space. Our Algorithm Path has the potential to explore the space of sparse models at a much finer resolution than other fast methods for approximating regularization paths over a grid of λ values. Consider that as the later are computed over M, typically M = 100, λ values, these can yield an upper bound of M distinct active sets; usually these yield much less than M distinct models. In contrast, our Algorithm Path yields an upper bound of K distinct models where K is the number of iterations needed, depending on the step-size t, to fully explore the sequence of sparse models. As K will often be much greater than M, our Algorithm Path will often explore a sequence of many more active sets and at a finer resolution than comparable methods. 3. Computationally fast. Our Algorithm Path has the potential to yield a sequence of sparse solutions much faster than other methods for computing regularization paths. Consider that our algorithm takes at most K iterations. In contrast, regularization paths of a grid of M λ values require M times the number of iterations needed to fully estimateβ β β (λ j ); often this will be much larger than K. In each iteration of our algorithm, the β β β and z subproblems require the most computational time. The β β β subproblem consists of the loss function with a quadratic penalty which can be solved via an analytical form for many loss functions. The z subproblem has the form of the proximal operator of P [29]: prox λ P (x) = arg min u x − u 2 2 + λ P(u). For many popular convex penalty-types such as the 1 -norm, group lasso, and nuclear norm, the proximal operator has an analytical solution. Thus, for a large number of statistical ma-chine learning problems, the iterations of our Algorithm Path are inexpensive to compute.
Overall, our Algorithmic Regularization Paths give a novel method for finding a sequence of sparse solutions by approximating the active sets of regularization paths. Our methods can be used in place of regularization paths for model selection purposes with many sparse statistical learning problems. In this chapter, instead of studying the mathematical and statistical properties of our new Algorithm Paths, which we leave for future work, we study our method through applications to several statistical learning problems in the next section.

Examples
To demonstrate the versatility and advantages of our ADMM Algorithmic Regularization Paths, we present several example applications to sparse statistical learning methods: sparse linear regression, reduced-rank multi-task learning and convex clustering.
Let us first discuss computational aspects of our Algorithm Path for sparse linear regression. Notice that the β β β -subproblem consists of solving a ridge-like regression problem. Much of the computations involved, however, can be pre-computed, specifically the matrix inversion, (X T X /n + I) −1 , and matrix-vector multiplication, X T y. In cases where p n, inverting a p × p matrix is is highly computationally intensive, requiring O(p 3 ) operations. We can reduce the computational cost to O(n 3 ), however, by invoking the Woodbury Matrix Identity [13]: X T X /n + I p −1 = I p − X T (n I n + X X T ) −1 X and caching the Cholesky decomposition of the smaller n-by-n matrix n I n + X X T . Thus, the iterative updates for β β β k are reduced to O(n 2 ), the cost of solving two n-by-n triangular linear systems of equations. The z-subproblem is solved via soft-thresholdingwhich requires only O(p) operations. We study our Algorithmic Regularization Path for sparse linear regression through a real data example. We use the publicly available 14-cancer microarray data from [15] to form our covariate matrix. This consists of gene expression measurements for n = 198 subjects and 16063 genes; we randomly sample p = 2000 genes to use as our data matrix X. We simulate sparse true signal β β β * with s = 16 non-zero features of absolute magnitude 5-10, and with the signs of the non-zero signals assigned randomly; the 16 non-zero variables were randomly chosen from the 2000 genes. The response variable y is generated as y = X β β β * +ε, where ε i.i.d.
∼ N(0, 1). A visualization of regularization paths, stability paths, and our Algorithmic Regularization Paths is given in Figure 2 for this example.
First, we verify empirically that that our ADMM algorithm with warm starts is equivalent to the regularization path (top left and top middle). Additionally, notice that, as expected, our Algorithm Path with a tiny step size (bottom right) also well approximates the sequence of active sets given by the regularization paths. With a larger step-size, however, our algorithm path (bottom left) yields a sequence of sparse models that differ markedly from the sparsity patterns of the regularization paths. This occurs as the change in regularization levels of each step are large enough so that the sparsity levels of the z-subproblem after the one-step approximation are not equivalent to that of the solution to (2).
Despite this, Figure 2 suggests that our Algorithm Paths with larger step sizes may have some advantages in terms of variable selection. Notice that regularization paths select many false positives (blue and gray dashed lines) before the true positives (red lines). This is expected as we used a real microarray data set for X consisting of strongly correlated variables that directly violate the irrepresentable conditions under which variable selection for sparse regression is achievable [4]. Our method, however, selects several true variables before the first false positive enters the model. To understand this further, we compare our approach to the Stability Paths used for stability selection [23], a re-sampling scheme with some of the strongest theoretical guarantees for variable selection. The stability paths, however, also select several false positives. This as well as other empirical results that are omitted for space reasons suggest that our Algorithm Path with moderate or larger step sizes may perform better than convex optimization techniques in terms of vari- able selection. While a theoretical investigation of this is beyond the scope of this book chapter, the intuition for this is readily apparent. Our Algorithm Path starts from a dense solution and uses a ridge-like penalty. Thus, coefficients of highly correlated variables are likely to be grouped and have similar magnitude coefficient values. When soft-thresholding is performed in the z-subproblem, variables which are strongly correlated are likely to remain together in the model for at least the first several algorithm iterations. By keeping correlated variables together in the model longer and otherwise eliminating irrelevant variables, this gives our algorithm a better chance of selecting the truly non-zero variables out of a correlated set. Hence, the fact that we start with a dense solution seems to help us; this is in contrast to the LASSO, LAR and OMP paths which are initialized with an empty active set and greedily add variables most correlated with the response [28,8]. We plan on investigating our methods in terms of variable selection in future work. Finally, we compare our Algorithm Paths to state-of-the-art methods for computing the sparse regression regularization paths in terms of computational time in Table 1. The regularization paths were computed using the glmnet R package [10] which is based on shooting (coordinate descent) routines [10]. This approach and software is widely regarded as one of the fastest solvers for sparse regression. Notice that our Algorithm Paths, coded entirely in Matlab, run in about a fifth of the time as this state-of-the-art competitor. Also, our computational time is far superior to the re-sampling schemes required to compute the stability paths.
Overall, our Algorithmic Regularization Path for sparse linear regression reveals major computational advantages for finding a sequence of sparse models that approximate the active sets of regularization paths. Additionally, empirical evidence suggests that our methods may also enjoy some important statistical advantages in terms of variable selection that we will explore in future work.

Reduced-Rank Multi-Task Learning
Our ADMM Algorithmic Regularization Path applies generally to many convex penalty types beyond the 1 -norm. Here, we demonstrate our method in conjunction with a reduced-rank multi-task learning problem also called multi-response regression. This problem has been studied by [27] among many others.
Suppose we observe n iid samples measured on p covariates and for q outcomes, yielding a covariate matrix, X ∈ ℜ n×p , and a response matrix Y ∈ ℜ n×q . Then, our goal is to fit the following statistical model: Y = X B +ε, where B is the p × q coefficient matrix which we seek to learn, and ε is independent noise. As often the number of covariates is large relative to the sample size, pq n, many have suggested to regularize the coefficient matrix B by assuming it has a lowrank structure, rank(B) < p ∧ q. Thus, our model space of sparse solutions is given by the space of all possible reduced-rank solutions. Exploring this space is an NP hard computational problem; thus, many have employed the nuclear norm penalty, B * = ∑ p∧q j=1 σ j (B), which is the sum (or 1 -norm) of the singular values of B, σ (B), and the tightest convex relaxation of the rank constraint. Thus, we arrive at the following optimization problem: Here, · F is the Frobenious norm, λ ≥ 0 is the regularization parameter controlling the rank of the solution and · * is the nuclear norm penalty. For model selection then, one seeks to explore the sequence of low-rank solutions obtained as λ varies. To develop our Algorithm Path for approximating this sequence of low-rank solutions, let us consider the ADMM sub-problems for solving (4). The augmented Lagrangian, sub-problems, dual updates are analogous to that of the sparse linear regression example, and hence we omit these here. Examining the Z-subproblem, however, recall that this is the proximal operator for the nuclear norm penalty: Z k = arg min Output {Z k : k = 1, · · · , K} as the algorithmic regularization path.
Then, following the framework of the sparse linear regression example, our ADMM Algorithmic Regularization Path for the reduced-rank mutli-task learning (regression) is outlined in Algorithm 4. Notice that the algorithm has the same basic steps as in the previous example except that solving the proximal operator for the Z sub-problem entails singular value thresholding. This step is the most computationally time consuming aspect of the algorithm as K total SVDs must be computed to approximate the sequence of solutions. Also note that similarly to the sparse regression example, the inversion needed, (X T X /n + I) −1 , can be precomputed by using the matrix inversion identities as previously discussed and cached as a convenient factorization; hence, this is computationally feasible even when p n.
To demonstrate the computational advantages of our approach, we conduct a small simulation study comparing our method to the two most commonly used algorithms for reduced-rank regression: proximal gradient descent and ADMM. First, we generate data according to the model: Y = X B +ε, where X 200×100 is generated as independent standard Gaussians, B 100×100 is an image of the Batman symbol, and ε 200×100 is independent standard Gaussian noise. We set the signal in the coefficient matrix to be a low-rank image of the Batman symbol, rank(B) = 38, which  can be well-approximated by further reduced rank images. We applied our Algorithmic Regularization Paths to this simulated example with 500 logarithmicallyspaced values of γ. Results are given in Figure 3 and show that our Algorithm Path smoothly explores the model space of reduced rank solutions and nicely approximates the true signal as a low-rank batman image. We also conduct a timing comparison to implementations of proximal gradient descent and ADMM algorithms using warm-starts for this same example; results are given in Table 2. Here, we see that our approach requires much fewer SVD computations and is much faster than both algorithms, especially the ADMM algorithm. Additionally, both the ADMM and proximal gradient algorithm employed 100 logarithmically spaced values of the regularization parameter, λ . With this, however, we see that not all possible ranks of the model space are considered, with proximal gradient and ADMM considering 57 and 51 ranks out of 100 respectively. In contrast, our Algorithmic Regularization Path yields a sequence of sparse solutions at a much finer resolution, considering 90 out of the 100 possible ranks. Thus, for proximal gradient and ADMM algorithms to consider the same range of possible sparsity levels (ranks), a greater number of problems would have to be solved over a much finer grid of regularization parameters, further inflating compute times.
Overall, our approach yields substantial computational savings for computing a sequence of sparse solutions for reduced rank regression compared to other stateof-the-art methods for this problem.

Convex Clustering
Our final example applies the ADMM Algorithmic Regularization path to an example with fusion type or non-separable penalties, namely a recently introduced convex formulation of cluster analysis [5,17,19]. Given n points y 1 , . . . , y n in ℜ p , we pose the clustering problem as follows. Assign to each point y i its own cluster center β β β i ∈ ℜ p . We then seek an assignment of β β β i that minimizes the distances between y i and β β β i and seeks sparsity between cluster center pairs β β β i and β β β j . Computing all possible cluster assignments, however, is an NP hard problem. Hence, the following relaxation poses finding the cluster assignments as a convex optimization problem: where λ is a positive regularization parameter, and w i j is a nonnegative weight. When λ = 0, the minimum is attained when β β β i = y i , and each point occupies a unique cluster. As λ increases, the cluster centers begin to coalesce. Two points y i and y j with β β β i = β β β j are said to belong to the same cluster. For sufficiently large λ all points coalesce into a single cluster at y, the mean of the y i . Because the objective in (5) is strictly convex and coercive, it possesses a unique minimizer for each value of λ . This is in stark contrast to other typical criteria used for clustering, which often rely on greedy algorithms that are prone to get trapped in suboptimal local minima. Because of its coalescent behavior, the resulting solution path can be considered a convex relaxation of hierarchical clustering [17]. This problem generalizes the fused LASSO [34], and as with other fused LASSO problems, penalizing affine transformations of the decision variable makes minimization challenging in general. The one exception is when a 1-norm is used instead of the 2-norm in the fusion penalty terms. In this case, the problem reduces to a weighted one-dimensional total variation denoising problem. Under other norms, including the 2-norm, the situation, is salvageable if we adopt a splitting strategy discussed earlier in Section 1.1 for dealing with fusion type or non-separable penalties. Briefly, we consider using the 2-norm in the fusion penalty to be most broadly applicable since the solutions to the convex clustering problem become invariant to rotations in the data. Consequently, clustering assignments will also be guaranteed to be rotationally invariant.
Let the variables z i j ∈ ℜ p record the differences between the ith and jth points. We denote the collections of variables {β β β i } n i=1 and {z z z i j } i< j by β β β and z respectively. Then the original problem can be reformulated as: Consider the ADMM algorithm derived in [5] for solving (6). Let u i j ∈ ℜ p denote the Lagrange multiplier for the i jth equality constraint. Let u denote the collection of variables {u i j } i< j . The augmented Lagrangian is given by: Then, the three ADMM subproblems are given by: Splitting the variables in this manner allows us to solve a series of straightforward subproblems. Updating β β β involves solving a ridge regression problem. Despite the fact that the quadratic penalty term is not separable in the β β β , after some algebraic maneuvering, which is detailed in [5], it is possible to explicitly write down the updates for each β β β separately: Updating z requires minimizing an objective that separates in each of the z i j , This step can be computed explicitly using the block-wise soft-thresholding operator, the proximal operator of the group LASSO [41], namely, where a + = max(a, 0) and τ ≥ 0 controls the amount of shrinkage towards zero. For model selection purposes, one typically studies the sequence of cluster assignments given by coalescent patterns of β β β , or the sparse patterns in the first differences of β β β , as λ varies. We then seek to quickly approximate this sequence of active sets given by the coalescent patterns of β β β with our Algorithmic Regularization Paths, summarized in Algorithm 5.
As in the general case, we can use iterates of the z-subproblem to approximate a sparse sequence of cluster assignments. Given z, we can determine a clustering assignment in time that is linear in the number of data points n. We simply apply breadth-first search to identify the connected components of the following graph induced by the z. The graph identifies a node with every data point and places an edge between the ith and jth node if and only if z i j = 0. Each connected component corresponds to a cluster.
We now illustrate on a simulated "halfmoon" data set of n = 200 points in ℜ 2 , that computing our Algorithm Path can lead to non-trivial computational cost savings for obtaining a sequence of clustering assingments. We first detail some preliminaries. Although we do not take the space to discuss it here, in practice the choice of weights is very important. This topic is explored in [5], and we use the sparse kernel weights which were shown to work well empirically in that paper. We created a geometric sequence of parameters λ (k) and γ (k) , namely given a fixed multiplicative factor t > 1, we set λ (k+1) = tλ (k) . The sequence γ (k) was constructed similarly, although we study our Algorithm Paths for several multiplicative factors, t ∈ {1.1, 1.05, 1.01}.
In contrast to the regularization path, the Algorithm Path does not require any convergence checks since only one step is taken at each grid point. Nonetheless, we only report the number of rounds of ADMM updates taken by each approach. The Algorithm Path took 259, 1294, and 2, 536 rounds of updates for the three step sizes considered; in contrast, the regularization path even for a very modest tolerance level, 10 −4 , required a grand total of 30,008 rounds of updates, substantially more than our approach. Figure 4 shows the ADMM Algorithmic Regularization paths and regularization path respectively for this simulated example. For each data point i we plot the sequence of the segments between consecutive estimates of its center, namely β β β k+1 i and β β β k i . These paths begin to overlap and merge into "trunks" when center estimates for close-by data points begin to coincide as the parameters λ (k) and γ (k) becomes sufficiently large. For sufficiently small step sizes for the regularization levels the Algorithm Path and regularization path are strikingly similar, as expected and demonstrated previously in our other examples. For larger step sizes, however, the paths differ markedly, but still appear to capture the same clustering assignments. Overall, although the simulated data is relatively small, computing the whole regularization path, even for a modest stopping tolerance can, requires an order of magnitude more iterations and computational time than the Algorithm Path.

Discussion
In this chapter, we have presented a novel framework for approximating the sequence of active sets associated with regularization paths of sparse statistical learning problems. Instead of solving optimization problems over a grid of penalty parameters as in traditional regularization paths, our algorithm performs a series of one-step approximations to an ADMM algorithm employing warm-starts with the goal of estimating a good sequence of sparse models. Our approach has a number of advantages including easy implementation, exploration of the sparse model space at a fine resolution, and most importantly fast compute times; we have demonstrated these advantages through several sparse statistical learning examples.
In our demonstrations, we have focused simply on computing the full sequence of active sets corresponding to the regularization path which is the critical computationally intensive step in the process of model selection. Once the sequence of sparse models has been found, common methods for model selection such as AIC, BIC, cross-validation and stability selection, can be employed to choose the optimal model. We note that with regularization paths, model selection procedures typically choose the optimal λ which indexes the optimal sparse model. For our Algorithm Paths which do not directly solve regularized statistical problems, model selection procedures should be used to choose the optimal iteration, k, and the corresponding sparse model given by the active set of z k . While this chapter has focused on finding the sequence of sparse models via our Algorithm Paths, we plan to study using these paths in conjunction with common model selection procedures in future work.
As the ADMM algorithm has been widely used for sparse statistical learning problems, the mechanics are in place for broad application of our Algorithm Paths which utilize the three standard ADMM subproblems. Indeed, our approach could potentially yield substantial computational savings for any ADMM application where the β β β and z can be solved efficiently. Furthermore, there has been much recent interest in distributed versions of ADMM algorithms [39,26]. Thus, there is the potential to use these in conjunction with our problem to distribute computation in the β β β and z subproblems and further speed computations for Big-Data problems. Also, we have focused on developing our Algorithm Path for sparse statistical learning problems that can be written as a composite of a smooth loss function and a non-smooth, convex penalty. Our methods, however, can be easily extended to study constrained statistical learning problems, such as that of the support vector machines. Finally, our framework utilizes the ADMM splitting method, but the strategies we develop could also be useful for computing a sequence of sparse models using other operator splitting algorithms.
Our work raises many questions from statistical and optimization perspectives. Further work needs to be done to characterize and study the mathematical properties of the Algorithm Paths as well as relate them to existing optimization procedures and algorithms. For example, ADMM is just one of many variants of proximal methods [29]. We suspect that other variants, such as proximal gradient descent, used to fit sparse models will also benefit from an Algorithm Path approach in expediting the model selection procedure. We leave this as future work.
In our demonstrations in Section 3, we suggested empirically that our Algorithm Paths with a tiny step size closely approximate the sequence of active sets associated with regularization paths. Further work needs to be done to verify this connection mathematically. Along these lines, a key practical question is how to choose the appropriate step size for increasing the amount of regularization as the algorithm progresses. As we have demonstrated, changing the step size yields paths with very different solutions and behaviors that warrant further investigation. For now, our recommendation is to employ a fairly small step size as these well-approximate the traditional regularization paths in all of the examples we have studied. Additionally, our approach may be related to other new proposals for computing regularization paths based on partial differential equations, for example [32]; these potential connections merit further investigation.
Our work also raises a host of interesting statistical questions as well. The sparse regression example suggested that Algorithm Paths may not simply yield computational savings, but may also perform better in terms of variable selection. This raises an interesting statistical prospect that we plan to carefully study in future work.
To conclude, we have introduced a novel approach to approximating the sequence of active sets associated with regularization paths for large-scale sparse statistical learning procedures. Our methods yield substantial computational savings and raise a number of interesting open questions for future research.