Smoothing proximal gradient method for general structured sparse regression

Introduction

When the structure to be imposed during shrinkage has a relatively simple form, such as nonoverlapping groups over variables (e.g., group lasso [Yuan and Lin 2006]) or a linear-ordering (a.k.a., chain) of variables (e.g., fused lasso [Tibshirani and Saunders 2005]), efficient optimization methods have been developed. For example, under group lasso, due to the separability among groups, a proximal operator The proximal operator associated with the penalty is defined as $\mathop{\arg\min}_{{\bm{\beta}}}\frac{1}{2}\|{\bm{\beta}}-{\mathbf{v}}\|_{2}^{2}+P({\bm{\beta}})$, where ${\mathbf{v}}$ is any given vector and $P({\bm{\beta}})$ is the nonsmooth penalty. associated with the penalty can be computed in closed-form; thus, a number of composite gradient methods [Beck and Teboulle 2009, Liu, Ji and Ye 2009, Nesterov 2007] that leverage the proximal operator as a key step (so-called “proximal gradient method”) can be directly applied. For fused lasso, although the penalty is not separable, a coordinate descent algorithm was shown feasible by explicitly leveraging the linear ordering of the inputs [Friedman et al. 2007].

Unfortunately, these algorithmic advancements have been outpaced by the emergence of more complex structures one would like to impose during shrinkage. For example, in order to handle a more general class of structures such as a tree or a graph over variables, various regression models that further extend the group lasso and fused lasso ideas have been recently proposed. Specifically, rather than assuming the variable groups to be nonoverlapping as in the standard group lasso, the overlapping group lasso [Jenatton, Audibert and Bach 2009] allows each input variable to belong to multiple groups, thereby introducing overlaps among groups and enabling incorporation of more complex prior knowledge on the structure. Going beyond the standard fused lasso, the graph-guided fused lasso extends the original chain structure over variables to a general graph over variables, where the fused-lasso penalty is applied to each edge of the graph [Kim, Sohn and Xing 2009]. Due to the nonseparability of the penalty terms resultant from the overlapping group or graph structures in these new models, the aforementioned fast optimization methods originally tailored for the standard group lasso or fused lasso cannot be readily applied here, due to, for example, unavailability of a closed-form solution of the proximal operator. In principle, generic convex optimization solvers such as the interior-point methods (IPM) could always be used to solve either a second-order cone programming (SOCP) or a quadratic programming (QP) formulation of the aforementioned problems; but such approaches are computationally prohibitive for problems of even a moderate size. Very recently, a great deal of attention has been given to devise practical solutions to the complex structured sparse regression problems discussed above in statistics and the machine learning community, and numerous methods have been proposed [Duchi and Singer 2009, Jenatton et al. 2010, Liu, Yuan and Ye 2010, Mairal et al. 2010, Tibshirani and Taylor 2010, Zhou and Lange 2011]. All of these recent works strived to provide clever solutions to various subclasses of the structured sparsity-inducing penalties; but, as we survey in Section 4, they are still short of reaching a simple, unified and general solution to a broad class of structured sparse regression problems.

In this paper we propose a generic optimization approach, the smoothing proximal gradient (SPG) method, for dealing with a broad family of sparsity-inducing penalties of complex structures. We use the overlapping-group-lasso penalty and graph-guided-fused-lasso penalty mentioned above as our motivating examples. Although these two types of penalties are seemingly very different, we show that it is possible to decouple the nonseparable terms in both penalties via the dual norm; and reformulate them into a common form to which the proposed method can be applied. We call our approach a “smoothing” proximal gradient method because instead of optimizing the original objective function directly as in other proximal gradient methods, we introduce a smooth approximation to the structured sparsity-inducing penalty using the technique from Nesterov 2005. Then, we solve the smoothed surrogate problem by a first-order proximal gradient method known as the fast iterative shrinkage-thresholding algorithm (FISTA) [Beck and Teboulle 2009]. We show that although we solve a smoothed problem, when the smoothness parameter is carefully chosen, SPG achieves a convergence rate of $O(\frac{1}{\varepsilon})$ for the original objective for any desired accuracy $\varepsilon$. Below, we summarize the main advantages of this approach:

It is a first-order method, as it uses only the gradient information. Thus, it is significantly more scalable than IPM for SOCP or QP. Since it is gradient-based, it allows warm restarts, and thereby potentiates solving the problem along the entire regularization path [Friedman et al. 2007].

It is applicable to a wide class of optimization problems with a smooth convex loss and a nonsmooth nonseparable structured sparsity-inducing penalty. Additionally, it is applicable to both uni- and multi-task sparse structured regression, with structures on either (or both) inputs/outputs.

Theoretically, it enjoys a convergence rate of $O(\frac{1}{\varepsilon})$, which dominates that of the standard first-order method such as the subgradient method whose rate is of $O(\frac{1}{\varepsilon^{2}})$.

Finally, SPG is easy to implement with a few lines of MATLAB code.

The rest of this paper is organized as follows. In Section 2 we present the formulation of overlapping group lasso and graph-guided fused lasso. In Section 3 we present the SPG method along with complexity results. In Section 4 we discuss the connections between our method and MM, and comparisons with other related methods. In Section 5 we extend our algorithm to multivariate-task regression. In Section 6 we present numerical results on both simulated and real data sets, followed by conclusions in Section 7. Throughout the paper, we will discuss overlapping-group-lasso and graph-guided-fused-lasso penalties in parallel to illustrate how the SPG can be used to solve the corresponding optimization problems generically.

Background: Linear regression regularized by structured sparsity-inducing penalties

We begin with a basic outline of the high-dimensional linear regression model, regularized by structured sparsity-inducing penalties.

The standard lasso penalty does not assume any structure among the input variables, which limits its applicability to complex high-dimensional scenarios in many applied problems. More structured constraints on the input variables such as groupness or pairwise similarities can be introduced by employing a more sophisticated sparsity-inducing penalty that induces joint sparsity patterns among related inputs. We generically denote the structured sparsity-inducing penalty by $\Omega({\bm{\beta}})$ without assuming a specific form, and define the problem of estimating a structured sparsity pattern of the coefficients as follows:

Overlapping-group-lasso penalty. Given prior knowledge of (possibly overlapping) grouping of variables or features, if it is desirable to encourage coefficients of features within the same group to be shrunk to zero jointly, then a composite structured penalty of the following form can be used:

Apparently, the penalty $\Omega({\bm{\beta}})\equiv\gamma\sum_{g\in\mathcal{G}}w_{g}\|{\bm{\beta}}_{g}\|_{2}$ alone enforces only group-level sparsity but not

where $\tau(r_{ml})$ represent a general weight function that enforces a fusion effect over coefficients $\beta_{m}$ and $\beta_{l}$ of relevant features. In this paper we consider $\tau(r)=|r|$, but any monotonically increasing function of the absolute values of correlations can be used.

The $\operatorname{sign}(r_{ml})$ in (4) ensures that two positively correlated inputs would tend to influence the output in the same direction, whereas two negatively correlated inputs impose opposite effect. Since the fusion effect is calibrated by the edge weight, the graph-guided-fused-lasso penalty in (4) encourages highly inter-correlated inputs corresponding to a densely connected subnetwork in $G$ to be jointly selected as relevant.

It is noteworthy that when $r_{ml}=1$ for all $e=(m,l)\in E$, and $G$ is simply a chain over nodes, we have

which is identical to the standard fused lasso penalty [Tibshirani and Saunders 2005].

Smoothing proximal gradient

Although (2) defines a convex program, of which a globally optimal solution to ${\bm{\beta}}$ is attainable, the main difficulty in solving (2) arises from the nonseparability of elements of ${\bm{\beta}}$ in the nonsmooth penalty function $\Omega({\bm{\beta}})$. As we show in the next subsection, although the overlapping-group-lasso and graph-guided-fused-lasso penalties are seemingly very different, we can reformulate the two types of penalties as a common matrix algebraic form, to which a generic Nesterov smoothing technique can be applied. The key in our approach is to decouple the nonseparable structured sparsity-inducing penalties into a simple linear transformation of ${\bm{\beta}}$ via the dual norm. Based on that, we introduce a smooth approximation to $\Omega({\bm{\beta}})$ using the technique from Nesterov 2005 such that its gradient with respect to ${\bm{\beta}}$ can be easily calculated.

In this section we show that utilizing the dual norm, the nonseparable structured sparsity-inducing penalty in both (3) and (4) can be decoupled; and reformulated into a common form as a maximization problem over the auxiliary variables.

Note that $C$ is a highly sparse matrix with only a single nonzero element in each row and $\sum_{g\in\mathcal{G}}|g|$ nonzero elements in the entire matrix, and, hence, can be stored with only a small amount of memory during the optimization procedure.

Reformulating graph-guided-fused-lasso penalty. First, we rewrite the graph-guided-fused-lasso penalty in (4) as follows:

2 Smooth approximation to structured sparsity-inducing penalty

The common formulation of $\Omega({\bm{\beta}})$ given above [i.e., $\Omega({\bm{\beta}})=\max_{{\bm{\alpha}}\in\mathcal{Q}}{\bm{\alpha}}^{T}C{\bm{\beta}}$] is still a nonsmooth function of ${\bm{\beta}}$, and this makes the optimization challenging. To tackle this problem, using the technique from Nesterov 2005, we construct a smooth approximation to $\Omega({\bm{\beta}})$ as follows:

where $\mu$ is a positive smoothness parameter and $d({\bm{\alpha}})$ is a smoothing function defined as $\frac{1}{2}\|{\bm{\alpha}}\|_{2}^{2}$. The original penalty term can be viewed as $f_{\mu}({\bm{\beta}})$ with $\mu=0$; and one can verify that $f_{\mu}({\bm{\beta}})$ is a lower bound of $f_{0}({\bm{\beta}})$. In order to bound the gap between $f_{\mu}({\bm{\beta}})$ and $f_{0}({\bm{\beta}})$, let $D=\max_{{\bm{\alpha}}\in\mathcal{Q}}d({\bm{\alpha}})$. In our problems, $D=|\mathcal{G}|/2$ for the overlapping-group-lasso penalty and $D=|E|/2$ for the graph-guided-fused-lasso penalty. Then, it is easy to verify that the maximum gap between $f_{\mu}({\bm{\beta}})$ and $f_{0}({\bm{\beta}})$ is $\mu D$:

From Theorem 1 as presented below, we know that $f_{\mu}({\bm{\beta}})$ is a smooth function for any $\mu>0$. Therefore, $f_{\mu}({\bm{\beta}})$ can be viewed as a smooth approximation to $f_{0}({\bm{\beta}})$ with a maximum gap of $\mu D$; and the $\mu$ controls the gap between $f_{\mu}({\bm{\beta}})$ and $f_{0}({\bm{\beta}})$. Given a desired accuracy $\varepsilon$, the convergence result in Section 3.5 suggests $\mu=\frac{\varepsilon}{2D}$ to achieve the best convergence rate.

Now we present the key theorem [Nesterov 2005] to show that $f_{\mu}({\bm{\beta}})$ is smooth in ${\bm{\beta}}$ with a simple form of the gradient.

For any $\mu>0$, $f_{\mu}({\bm{\beta}})$ is a convex and continuously-differentiable function in ${\bm{\beta}}$, and the gradient of $f_{\mu}({\bm{\beta}})$ takes the following form:

where ${\bm{\alpha}}^{\ast}$ is the optimal solution to (10). Moreover, the gradient $\nabla f_{\mu}({\bm{\beta}})$ is Lipschitz continuous with the Lipschitz constant $L_{\mu}=\frac{1}{\mu}\|C\|^{2}$, where

$\|C\|$ is the matrix spectral norm of $C$ defined as $\|C\|\equiv{\max_{\|{\mathbf{v}}\|_{2}\leq 1}}\|C{\mathbf{v}}\|_{2}$.

$f_{\mu}({\bm{\beta}})$ as the Fenchel conjugate of $d(\cdot)$ at $\frac{C{\bm{\beta}}}{\mu}$, the smoothness can be obtained by applying Theorem 26.3 in Rockafellar 1996. The gradient in (11) can be derived from the Danskin’s theorem [Bertsekas 1999] and the Lipschitz constant is shown in Nesterov 2005. The details of the proof are given in the Appendix.

we introduce $d(\alpha)=\frac{1}{2}\alpha^{2}$, let $z_{s}(\alpha,\beta)\equiv\alpha\beta-\frac{1}{2}\alpha^{2}$ and $f_{\mu}(\beta)=\max_{\alpha\in}z_{s}(\alpha,\beta)$. The three-dimensional plot for $z_{s}(\alpha,\beta)$ with $\alpha$ restricted to $$is shown in Figure 1(c). Similarly, we project the surface in Figure 1(c) onto the$\beta-z_{s}$space as shown in Figure 1(d). For fixed$\beta$, the value of$f_{\mu}(\beta)$is the highest point along the$z$-axis. In Figure 1(d), we can see that the sharp point at$\beta=0$is removed and$f_{\mu}(\beta)$ becomes smooth.

To compute the $\nabla f_{\mu}({\bm{\beta}})$ and $L_{\mu}$, we need to know ${\bm{\alpha}}^{\ast}$ and $\|C\|$. We present the closed-form equations for ${\bm{\alpha}}^{\ast}$ and $\|C\|$ for the overlapping-group-lasso penalty and graph-guided-fused-lasso penalty in the following propositions. The proof is presented in the Appendix.

${\bm{\alpha}}^{\ast}$ under overlapping-group-lasso penalty.

Let ${\bm{\alpha}}^{\ast}$, which is composed of $\{{\bm{\alpha}}^{\ast}_{g}\}_{g\in\mathcal{G}}$, be the optimal solution to (10) for the overlapping-group-lasso penalty in (3). For any ${g\in\mathcal{G}}$,

${\bm{\alpha}}^{\ast}$ under graph-guided-fused-lasso penalty.

Let ${\bm{\alpha}}^{\ast}$ be the optimal solution of (10) for the graph-guided-fused-lasso penalty in (4). Then, we have

where $S$ is the projection operator defined as follows:

For any vector ${\bm{\alpha}}$, $S({\bm{\alpha}})$ is defined as applying $S$ on each and every entry of ${\bm{\alpha}}$.

$\|C\|$ is upper-bounded by $\sqrt{2\gamma^{2}\max_{j\in V}d_{j}}$, where

for $j\in V$ in graph $G$, and this bound is tight. Note that when $\tau(r_{e})=1$ for all $e\in E$, $d_{j}$ is simply the degree of the node $j$.

3 Smoothing proximal gradient descent

Given the smooth approximation to the nonsmooth structured sparsity-inducing penalties, now, we apply the fast iterative shrinkage-thresholding algorithm (FISTA) [Beck and Teboulle 2009, Tseng 2008] to solve a generically reformulated optimization problem, using the gradient information from Theorem 1. We substitute the penalty term $\Omega({\bm{\beta}})$ in (2) with its smooth approximation $f_{\mu}({\bm{\beta}})$ to obtain the following optimization problem:

be the smooth part of $\widetilde{f}({\bm{\beta}})$. According to Theorem 1, the gradient of $h({\bm{\beta}})$ is given as

Moreover, $\nabla h({\bm{\beta}})$ is Lipschitz-continuous with the Lipschitz constant,

where $\lambda_{\max}({\mathbf{X}}^{T}{\mathbf{X}})$ is the largest eigenvalue of $({\mathbf{X}}^{T}{\mathbf{X}})$.

adopt FISTA [Beck and Teboulle 2009, Tseng 2008] to minimize $\widetilde{f}({\bm{\beta}})$ as shown in Algorithm 1. Algorithm 1 alternates between the sequences $\{w^{t}\}$ and $\{{\bm{\beta}}^{t}\}$ and $\theta_{t}$ can be viewed as a special “step-size,” which determines the relationship between $\{w^{t}\}$ and $\{{\bm{\beta}}^{t}\}$ as in Step 4 of Algorithm 1. As shown in Beck and Teboulle 2009, such a way of setting $\theta_{t}$ leads to Lemma 1 in the Appendix, which further guarantees the convergence result in Theorem 2.

Rewriting $Q_{L}({\bm{\beta}},{\mathbf{w}}^{t})$ in (1),

Letting ${\mathbf{v}}=({\mathbf{w}}^{t}-\frac{1}{L}\nabla h({\mathbf{w}}^{t}))$, the closed-form solution for ${\bm{\beta}}^{t+1}$ can be obtained by soft-thresholding [Friedman et al. 2007] as presented in the next proposition.

can be obtained by the soft-thresholding operation:

4 Issues on the computation of the Lipschitz constant

When $J$ is large, the computation of $\lambda_{\max}({\mathbf{X}}^{T}{\mathbf{X}})$ and hence the Lipschitz constant $L$ could be very expensive. To further accelerate Algorithm 1, a line search backtracking step could be used to dynamically assign a constant $L_{t}$ for the proximal operator in each iteration [Beck and Teboulle 2009]. More specifically, given any positive constant $R$, let

The key to guarantee the convergence rate of Algorithm 1 is to ensure that the following inequality holds for each iteration:

It is easy to check that when $R$ is equal to the Lipschitz constant $L$, it will satisfy the above inequality for any ${\bm{\beta}}^{t+1}$ and ${\mathbf{w}}^{t}$. However, when it is difficult to compute the Lipschitz constant, instead of using a global constant $L$, we could find a sequence $\{L_{t}\}_{t=0}^{T}$ such that $L_{t+1}$ satisfies the inequality (19) for the $t$th iteration. In particular, we start with any small constant $L_{0}$. For each iteration, we find the smallest integer $a\in\{0,1,2,\ldots\}$ such that by setting $L_{t+1}=\tau^{a}L_{t}$, where $\tau>1$ is a predefined scaling factor, we have

Then we set ${\bm{\beta}}^{t+1}={\bm{\beta}}_{L_{t+1}}({\mathbf{w}}^{t})\equiv\mathop{\arg\min}Q_{L_{t+1}}({\bm{\beta}},{\mathbf{w}}^{t})$.

5 Convergence rate and time complexity

Although we optimize the approximation function $\widetilde{f}({\bm{\beta}})$ rather than the original $f({\bm{\beta}})$ directly, it can be proven that $f(\widehat{{\bm{\beta}}})$ is sufficiently close to the optimal objective value of the original function $f({\bm{\beta}}^{*})$. The convergence rate of Algorithm 1 is presented in the next theorem.

Let ${\bm{\beta}}^{\ast}$ be the optimal solution to (2) and ${\bm{\beta}}^{t}$ be the approximate solution at the $t$th iteration in Algorithm 1. If we require $f({\bm{\beta}}^{t})-f({\bm{\beta}}^{\ast})\leq\varepsilon$ where $f$ is the original objective, and set $\mu=\frac{\varepsilon}{2D}$, then the number of iterations $t$ is upper-bounded by

idea behind the proof of this theorem is to decompose $f({\bm{\beta}}^{t})-f({\bm{\beta}}^{\ast})$ into three parts: (i) $f({\bm{\beta}}^{t})-\widetilde{f}({\bm{\beta}}^{t})$, (ii) $\widetilde{f}({\bm{\beta}}^{t})-\widetilde{f}({\bm{\beta}}^{\ast})$ and (iii) $\widetilde{f}({\bm{\beta}}^{\ast})-f({\bm{\beta}}^{\ast})$. (i) and (iii) can be bounded by the gap of the approximation $\mu D$; and (ii) only involves the function $\widetilde{f}$ and can be upper bounded by $O(\frac{1}{t^{2}})$ as shown in Beck and Teboulle 2009. We obtain (21) by balancing these three terms. The details of the proof are presented in the Appendix. According to Theorem 2, Algorithm 1 converges in $O(\frac{\sqrt{2D}}{\varepsilon})$ iterations, which is much faster than the subgradient method with the convergence rate of $O(\frac{1}{\varepsilon^{2}})$. Note that the convergence rate depends on $D$ through the term $\sqrt{2D}$, and the $D$ depends on the problem size.

Since there is no line search in Algorithm 1, we cannot guarantee that the objective values are monotonically decreasing over iterations theoretically. But empirically, based on our own experience, the objective values always decrease over iterations. One simple

strategy to guarantee the monotone decreasing property is to first compute $\widetilde{{\bm{\beta}}}{}^{t+1}=\mathop{\arg\min}_{{\bm{\beta}}}Q_{L}({\bm{\beta}},{\mathbf{w}}^{t})$ and then set ${\bm{\beta}}^{t+1}=\mathop{\arg\min}_{{\bm{\beta}}\in\{\widetilde{{\bm{\beta}}}{}^{t+1},{\bm{\beta}}^{t}\}}f({\bm{\beta}})$.

As for the time complexity, the main computational cost in each iteration comes from calculating the gradient $\nabla h({\mathbf{w}}_{t})$. Therefore, SPG shares almost the same per-iteration time as the subgradient descent but with a faster convergence rate. In more details, if $J<N$ and ${\mathbf{X}}^{T}{\mathbf{X}}$ and ${\mathbf{X}}^{T}{\mathbf{y}}$ can be pre-computed and stored in memory, the computation of the first part of $\nabla h({\mathbf{w}}_{t})$, $({\mathbf{X}}^{T}{\mathbf{X}}){\mathbf{w}}_{t}-({\mathbf{X}}^{T}{\mathbf{y}})$, takes the time complexity of $O(J^{2})$. Otherwise, if $J>N$, we can compute this part by ${\mathbf{X}}^{T}({\mathbf{X}}{\mathbf{w}}_{t}-{\mathbf{y}})$, which takes the time complexity of $O(JN)$. As for the generic solver, IPM for SOCP for overlapping group lasso or IPM for QP for graph-guided fused lasso, although it converges in fewer iterations [i.e., $\log(\frac{1}{\varepsilon})$], its per-iteration complexity is higher by orders of magnitude than ours as shown in Table 1. In addition to time complexity, IPM requires the pre-storage of ${\mathbf{X}}^{T}{\mathbf{X}}$ and each IPM iteration requires significantly more memory to store the Newton linear system. Therefore, the SPG is much more efficient and scalable for large-scale problems.

6 Summary and discussions

The insight of our work was drawn from two lines of earlier works. The first one is the proximal gradient methods (e.g., Nesterov’s composite gradient method [Nesterov 2007], FISTA [Beck and Teboulle 2009]. They have been widely adopted to solve optimization problems with a convex loss and a relatively simple nonsmooth penalty, achieving $O(\frac{1}{\sqrt{\varepsilon}})$ convergence rate. However, the complex structure of the nonseparable penalties considered in this paper makes it intractable to solve the proximal operator exactly. This is the challenge that we circumvent via smoothing.

Related optimization methods

Recently, many first-order approaches have been developed for various subclasses of overlapping group lasso and graph-guided fused lasso. Below, we provide a survey of these methods:

2 Related work for mixed-norm based group-lasso penalty

Table 4.2 summarizes the applicability, the convergence rate and the per-iteration time complexity for the available first-order methods for different subclasses of group lasso penalties. More specifically, the methods in the first three rows adopt the proximal gradient framework. The first column of these rows gives the solver for the proximal operator. Each entry in Table 4.2 contains the convergence rate and the per-iteration time complexity. For the sake of simplicity, for all methods, we omit the time for computing the gradient of the loss function which is required for all of the methods [i.e., $\nabla g({\bm{\beta}})$ with $O(J^{2})$]. The per-iteration time complexity in the table may come from the computation of the proximal operator or subgradient of the penalty. “N.A.” stands for “not applicable” or no guarantee in the convergence. As we can see from Table 4.2, although our method is not the most ideal one for some of the special cases, our method along with FOBOS [Duchi and Singer 2009] are the only generic first-order methods that can be applied to all subclasses of the penalties.

3 Related work for fused lasso

For the graph-guided-fused-lasso penalty, when the structure is a simple chain, the pathwise coordinate descent method [Friedman et al. 2007] can be applied. For the general graph structure, a first-order method that approximately solves the proximal operator was proposed in Liu, Yuan and Ye 2010. However, the convergence cannot be guaranteed due to the errors introduced in computing the proximal operator over iterations.

Recently, two different path algorithms have been proposed [Tibshirani and Taylor 2010, Zhou and Lange 2011] that can be used to solve the graph-guided fused lasso as a special case. Unlike the traditional optimization methods that solve the problem for a fixed regularization parameter, they solve the entire path of solutions, and, thus, have great practical advantages. In addition, for both methods, updating solutions from one hitting time to another is computationally very cheap. More specifically, a QR decomposition based updating scheme was proposed in Tibshirani and Taylor 2010 and the updating in Zhou and Lange 2011 can be done by an efficient sweep operation.

However, for high-dimensional data with $J\gg N$, the path algorithms can have the following problems:

For a general design matrix ${\mathbf{X}}$ other than the identity matrix, the method in Tibshirani and Taylor 2010 needs to first compute the pseudo-inverse of ${\mathbf{X}}\colon\;{\mathbf{X}}^{+}=({\mathbf{X}}^{T}{\mathbf{X}})^{+}{\mathbf{X}}^{T}$, which could be computationally expensive for large $J$.

The original version of the algorithms in Tibshirani and Taylor 2010 and Zhou and Lange 2011 requires that ${\mathbf{X}}$ has a full column rank. When $J>N$, although one can add an extra $\varepsilon\|{\bm{\beta}}\|_{2}^{2}$ term, this changes the original objective value especially when $\varepsilon$ is large. For smaller $\varepsilon$, the matrix $({\mathbf{X}}^{*})^{T}{\mathbf{X}}^{*}$ with ${\mathbf{X}}^{*}=\bigl[{{\mathbf{X}}\atop\varepsilon I}\bigr]$ is highly ill-conditioned; and hence computing its inverse as the initialization step in Tibshirani and Taylor 2010 is very difficult. There is no known result on how to balance this trade-off.

In both Tibshirani and Taylor 2010 and Zhou and Lange 2011, the authors extend their algorithm to deal with the case when ${\mathbf{X}}$ does not have a full column rank. The extended version requires a Gramm–Schmidt process as the initialization, which could take some extra time.

In Table 2 we present the comparisons for different methods. From our analysis, the method in Zhou and Lange 2011 is more efficient than the one in Tibshirani and Taylor 2010 since it avoids the heavy computation of the pseudo-inverse of ${\mathbf{X}}$. In practice, if ${\mathbf{X}}$ has a full column rank and one is interested in solutions on the entire path, the method in Zhou and Lange 2011 is very efficient and faster than our method. Instead, when $J\gg N$, the path following methods may require a time-consuming preprocessing procedure.

Extensions to multi-task regression with structures on outputs

The structured sparsity-inducing penalties as discussed in the previous section can be similarly used in the multi-task regression setting [Kim and Xing 2010, Kim, Sohn and Xing 2009], where the prior structural information is available for the outputs instead of inputs. For example, in genetic association analysis, where the goal is to discover few genetic variants or single nucleotide polymorphisms (SNPs) out of millions of SNPs (inputs) that influence phenotypes (outputs) such as gene expression measurements, the correlation structure of the phenotypes can be naturally represented as a graph, which can be used to guide the selection of SNPs as shown in Figure 2. Then, the graph-guided-fused-lasso penalty can be used to identify SNPs that are relevant jointly to multiple related phenotypes.

In a sparse multi-task regression with structure on the output side, we encounter the same difficulties of optimizing with nonsmooth and nonseparable penalties as in the previous section, and the SPG can be extended to this problem in a straightforward manner. Due to the importance of this class of problems and its applications, in this section, we briefly discuss how our method can be applied to the multi-task regression with structured-sparsity-inducing penalties.

Overlapping-group-lasso penalty in multi-task regression. We define the overlapping-group-lasso penalty for a structured multi-task regression as follows:

Graph-guided-fused-lasso penalty in multi-task regression. Assuming that a graph structure over the $K$ outputs is given as $G$ with a set of nodes $V=\{1,\ldots,K\}$, each corresponding to an output variable and a set of edges $E$, the graph-guided-fused-lasso penalty for a structured multi-task regression is given as

2 Smoothing proximal gradient descent

Using similar techniques in Section 3.1, $\Omega({\mathbf{B}})$ can be reformulated as

where $\langle\mathbf{U},\mathbf{V}\rangle\equiv\operatorname{Tr}(\mathbf{U}^{T}\mathbf{V})$ denotes a matrix inner product. $C$ is constructed in a similar way as in (7) or (8), just by replacing the index of the input variables with the output variables, and ${\mathbf{A}}$ is the matrix of the auxiliary variables.

Then we introduce the smooth approximation of (25):

where $d({\mathbf{A}})\equiv\frac{1}{2}\|{\mathbf{A}}\|_{F}^{2}$. Following a proof strategy similar to that in Theorem 1, we can show that $f_{\mu}({\mathbf{B}})$ is convex and smooth with gradient $\nabla f_{\mu}({\mathbf{B}})=({\mathbf{A}}^{\ast})^{T}C$, where ${\mathbf{A}}^{\ast}$ is the optimal solution to (26). The closed-form solution of ${\mathbf{A}}^{\ast}$ and the Lipschitz constant for $\nabla f_{\mu}({\mathbf{B}})$ can be derived in the same way.

By substituting $\Omega({\mathbf{B}})$ in (22) with $f_{\mu}({\mathbf{B}})$, we can adopt Algorithm 1 to solve (22) with convergence rate of $O(\frac{1}{\varepsilon})$. The per-iteration time complexity of SPG as compared to IPM for SOCP or QP formulation is presented in Table 3. As we can see, the per-iteration complexity for SPG is linear in $\max(|K|,\sum_{g\in\mathcal{G}}|g|)$ or $\max(|K|,|E|)$, while traditional approaches based on IPM scape at least cubically to the size of outputs $K$.

Experiment

In this section we evaluate the scalability and efficiency of the smoothing proximal gradient method (SPG) on a number of structured sparse regression problems via simulation, and apply SPG to an overlapping group lasso problem on real genetic data.

On an overlapping group lasso problem, we compare the SPG with FOBOS [Duchi and Singer 2009] and IPM for SOCP. We use the state-of-the-art MATLAB package SDPT3 [Tütüncü, Toh and Todd 2003] for SOCP. On a multi-task graph-guided fused lasso problem, we compare the running time of SPG with that of the FOBOS [Duchi and Singer 2009] and IPM for QP. We use the commercial package MOSEK (http://www.mosek.com/) for QP. The graph-guided fused lasso can also be solved by SOCP, but it is less efficient than QP. Note that for FOBOS, since the proximal operator associated with $\Omega({\bm{\beta}})$ cannot be solved exactly, we set the “loss function” to $l({\bm{\beta}})=g({\bm{\beta}})+\Omega({\bm{\beta}})$ and the penalty to $\lambda\|{\bm{\beta}}\|_{1}$. According to Duchi and Singer 2009, for the nonsmooth loss $l({\bm{\beta}})$, FOBOS achieves $O(\frac{1}{\varepsilon^{2}})$ convergence rate, which is slower than our method.

All experiments are performed on a standard PC with 4GB RAM and the software is written in MATLAB. The main difficulty in comparisons is a fair stopping criterion. Unlike IPM, SPG and FOBOS do not generate a dual solution and, therefore, it is not possible to compute a primal-dual gap, which is the traditional stopping criterion for IPM. Here, we adopt a widely used approach for comparing different methods in the optimization literature. Since it is well known that IPM usually gives a more accurate (i.e., lower) objective, we set the objective obtained from IPM as the optimal objective value and stop the first-order methods when the objective is below 1.001 times the optimal objective. For large data sets for which IPM cannot be applied, we stop the first-order methods when the relative change in the objective is below $10^{-6}$. In addition, maximum iterations are set to 20,000.

Since our main focus is on the optimization algorithm, for the purpose of simplicity, we assume that each group in the overlapping group lasso problem receives the same amount of regularization and, hence, set the weights $w_{g}$ for all groups to be 1. In principle, more sophisticated prior knowledge of the importance for each group can be naturally incorporated into $w_{g}$. In addition, we notice that each variable $j$ with the regularization $\lambda|\beta_{j}|$ in $\lambda\|{\bm{\beta}}\|_{1}$ can be viewed as a singleton group. To ease the tuning of parameters, we again assume that each group (including the singleton group) receives the same amount of regularization and, hence, constrain the regularization parameters $\lambda=\gamma$.

The smoothing parameter $\mu$ is set to $\frac{\varepsilon}{2D}$ according to Theorem 2, where $D$ is determined by the problem size. It is natural that for large-scale problems with large $D$, a larger $\varepsilon$ can be adopted without affecting the recovery quality significantly. Therefore, instead of setting $\varepsilon$, we directly set $\mu=10^{-4}$, which provided us with reasonably good approximation accuracies for different scales of problems based on our experience for a range of $\mu$ in simulations. As for FOBOS, we set the stepsize rate to $\frac{c}{\sqrt{t}}$ as suggested in Duchi and Singer 2009, where $c$ is carefully tuned to be $\frac{0.1}{\sqrt{NJ}}$ for univariate regression and $\frac{0.1}{\sqrt{NJK}}$ for multi-task regression.

We simulate data for a univariate linear regression model with the overlapping group structure on the inputs as described below. Assuming that the inputs are ordered, we define a sequence of groups of 100 adjacent inputs with an overlap of 10 variables between two successive groups so that

with $J=90|\mathcal{G}|+10$. We set $\beta_{j}=(-1)^{j}\exp(-(j-1)/100)$ for $1\leq j\leq J$. We sample each element of ${\mathbf{X}}$ from i.i.d. Gaussian distribution, and generate the output data from ${\mathbf{y}}={\mathbf{X}}{\bm{\beta}}+{\bm{\varepsilon}}$, where ${\bm{\varepsilon}}\sim N(0,I_{N\times N})$.

To demonstrate the efficiency and scalability of SPG, we vary $J$, $N$ and $\gamma$ and report the total CPU time in seconds and the objective value in Table 4. The regularization parameter $\gamma$ is set to either $|\mathcal{G}|/5$ or $|\mathcal{G}|/20$. As we can see from Table 4, first, both SPG and FOBOS are more efficient and scalable by orders of magnitude than IPM for SOCP. For larger $J$ and $N$, we are unable to collect the results for SOCP. Second, SPG is more efficient than FOBOS for almost all different scales of the problems. In some entries in Table 4, the Obj. from FOBOS is much larger than other methods. This is because that FOBOS has reached the maximum number of iterations before convergence. Instead, for our simulations, SPG generally converges in hundreds of, or, at most, a few thousand, iterations and never pre-terminates. Third, for SPG, a smaller $\gamma$ leads to faster convergence. This result is consistent with Theorem 2, which shows that the number of iterations is linear in $\gamma$ through the term $\|C\|$. Moreover, we notice that a larger $N$ does not increase the computational time for SPG. This is also consistent with the time complexity analysis, which shows that for linear regression, the per-iteration time complexity is independent of $N$.

However, we find that the solutions from IPM are more accurate and, in fact, it is hard for first-order approaches to achieve the same precision as IPM. Assuming that we require $\varepsilon=10^{-6}$ for the accuracy of the solution, it takes IPM about $O(\log(\frac{1}{\varepsilon}))\approx 14$ iterations to converge, while it takes $O(\frac{1}{\varepsilon})=10^{6}$ iterations for SPG. This is the drawback for any first-order method. However, in many real applications, we do not require the objective to be extremely accurate (e.g., $\varepsilon=10^{-3}$ is sufficiently accurate in general) and first-order methods are more suitable. More importantly, first-order methods can be applied to large-scale high-dimensional problems while IPM can only be applied to small or moderate scale problems due to the expensive computation necessary for solving the Newton linear system.

2 Simulation study II: Multi-task graph-guided fused lasso

We simulate data using the following scenario analogous to the problem of genetic association mapping, where we are interested in identifying a small number of genetic variations (inputs) that influence the phenotypes (outputs). We use $K=10$, $J=30$ and $N=100$. To simulate the input data, we use the genotypes of the 60 individuals from the parents of the HapMap CEU panel [The International HapMap Consortium 2005], and generate genotypes for an additional 40 individuals by randomly mating the original 60 individuals. We generate the regression coefficients ${\bm{\beta}}_{k}$’s such that the outputs $\mathbf{y}_{k}$’s are correlated with a block-like structure in the correlation matrix. We first choose input-output pairs with nonzero regression coefficients as we describe below. We assume three groups of correlated output variables of sizes 3, 3 and 4. We randomly select inputs that are relevant jointly among the outputs within each group, and select additional inputs relevant across multiple groups to model the situation of a higher-level correlation structure across two subgraphs as in Figure 3(a). Given the sparsity pattern of ${\mathbf{B}}$, we set all nonzero

To compare SPG with FOBOS and IPM for QP in solving such a structured sparse regression problem, we vary $K$, $J$, $N$ and present the computation time in seconds in Figures 4(a)–(c), respectively. We select the regularization parameter $\gamma$ using separate validation data, and report the CPU time for the graph-guided fused lasso with the selected $\gamma$. The input/output data and true regression coefficient matrix ${\mathbf{B}}$

are generated in a way similar as above. More precisely, we assume that each group of correlated output variables is of size 10. For each group of the outputs, we randomly select $10\%$ of the input variables as relevant. In addition, we randomly select $5\%$ of the input variables as relevant to every two consecutive groups of outputs and $1\%$ of the input variables as relevant to every three consecutive groups. We set the $\rho$ for each data item so that the number of edges is 5 times the number of the nodes (i.e., $|E|=5K$). Figure 4 shows that SPG is substantially more efficient and can scale up to very high-dimensional and large-scale data sets. Moreover, we notice that the increase of $N$ almost does not affect the computation time of SPG, which is consistent with the complexity analysis in Section 3.5.

3 Real data analysis: Pathway analysis of breast cancer data

In this section we apply the SPG to an overlapping group lasso problem with a logistic loss on real-world data collected from breast cancer tumors [Jacob, Obozinski and Vert 2009, van de Vijver 2002]. The main goal is to demonstrate the importance of employing structured sparsity-inducing penalties for performance enhancement in real life high-dimensional regression problems, thereby further exhibiting and justifying the needs of efficient solvers such as SPG for such problems.

The data are given as gene expression measurements for 8,141 genes in 295 breast-cancer tumors (78 metastatic and 217 nonmetastatic). A lot of research efforts in biology have been devoted to identifying biological pathways that consist of a group of genes participating in a particular biological process to perform a certain functionality in the cell. Thus, a powerful way of discovering genes involved in a tumor growth is to consider groups of interacting genes in each pathway rather than individual genes independently [Ma and Kosorok 2010]. The overlapping-group-lasso penalty provides us with a natural way to incorporate this known pathway information into the biological analysis, where each group consists of the genes in each pathway. This approach can allow us to find pathway-level gene groups of significance that can distinguish the two tumor types. In our analysis of the breast cancer data, we cluster the genes using the canonical pathways from the Molecular Signatures Database [Subramanian et al. 2005], and construct the overlapping-group-lasso penalty using the pathway-based clusters as groups. Many of the groups overlap because genes can participate in multiple pathways. Overall, we obtain 637 pathways over 3,510 genes, with each pathway containing 23.47 genes on average and each gene appearing in four pathways on average. Instead of analyzing all 8,141 genes, we focus on these 3,510 genes which belong to certain pathways. We set up the optimization problem of minimizing the logistic loss with the overlapping-group-lasso penalty to classify the tumor types based on the gene expression levels, and solve it with SPG.

Since the number of positive and negative samples are imbalanced, we adopt the balanced error rate defined as the average error rate of the two classes. See http://www.modelselect.inf.ethz.ch/evaluation.php for more details. We split the data into the training and testing sets with the ratio of $2:1$, and vary the $\lambda=\gamma$ from large to small to obtain the full regularization path.

We perform functional enrichment analysis on the selected pathways, using the functional annotation tool [Huang, Sherman and Lempicki 2009], and verify that the selected pathways are significant in their relevance to the breast-cancer tumor types. For example, in a highly sparse model obtained with the group-lasso penalty at the very left end of Figure 5(b), the selected gene markers belong to only seven pathways, and many of these pathways appear to be reasonable candidates for an involvement in breast cancer. For instance, all proteins in one of the selected pathways are involved in the activity of proteases, whose function is to degrade unnecessary or damaged proteins through a chemical reaction that breaks peptide bonds. One of the most important malignant properties of cancer involves the uncontrolled growth of a group of cells, and protease inhibitors, which degrade misfolded proteins, have been extensively studied in the treatment of cancer. Another interesting pathway selected by the overlapping group lasso is known for its involvement in nicotinate and nicotinamide metabolism. This pathway has been confirmed as a marker for breast cancer in previous studies [Ma and Kosorok 2010]. In particular, the gene ENPP1 (ectonucleotide pyrophosphatase/phosphodiesterase 1) in this pathway has been found to be overly expressed in breast tumors [Abate et al. 2005]. Other selected pathways include the one related to ribosomes and another related to DNA polymerase, which are critical in the process of generating proteins from DNA and relevant to the property of uncontrolled growth in cancer cells.

We also examine the number of selected pathways that give the lowest error rate in Figure 5. At the error rate of 29.23%, 125 pathways (696 genes) are selected. It is interesting to notice that among these 125 pathways, one is closely related to apoptosis, which is the process of programmed cell death that occurs in multicellular organisms and is widely known to be involved in uncontrolled tumor growth in cancer. Another pathway involves the genes BRCA1, BRCA2 and ATR, which have all been associated with cancer susceptibility.

Conclusions and future work

In this paper we investigated an optimization problem for estimating the structured-sparsity pattern in regression coefficients under a general class of structured sparsity-inducing penalties. Many of the structured sparsity-inducing penalties including the overlapping-group-lasso penalties and graph-guided-fused-lasso penalty share a common set of difficulties in optimization such as nonseparability and nonsmoothness. We showed that the optimization problems with these penalties can be transformed into a common form, and proposed a general optimization approach, called the smoothing proximal gradient method, for efficiently solving the optimization problem of this common form. Our results show that the proposed method enjoys both desirable theoretical guarantee and practical scalability under various difficult settings involving complex structure constraints, multi-task and high-dimensionality.

There are several future directions for this work. First, it is known that reducing $\mu$ over iterations leads to better empirical results. However, in such a scenario, the convergence rate is harder to analyze. Moreover, since the method is only based on gradient, its online version with the stochastic gradient descent can be easily derived. However, proving the regret bound will require a more careful investigation.

Appendix

We first introduce the concept of Fenchel conjugate.

The Fenchel conjugate of a function $\varphi({\bm{\alpha}})$ is the function $\varphi^{\ast}({\bm{\beta}})$ defined as

that $d({\bm{\alpha}})=\frac{1}{2}\|{\bm{\alpha}}\|^{2}$ with the $\operatorname{dom}({\bm{\alpha}})=\mathcal{Q}$. According to Definition 1, the conjugate of $d(\cdot)$ at $\frac{C{\bm{\beta}}}{\mu}$ is $d^{\ast}(\frac{C{\bm{\beta}}}{\mu})=\sup_{{\bm{\alpha}}\in\mathcal{Q}}({\bm{\alpha}}^{T}\frac{C{\bm{\beta}}}{\mu}-d({\bm{\alpha}}))$ and, hence,

According to Theorem 26.3 in Rockafellar 1996, “a closed proper convex function is essentially strictly convex if and only if its conjugate is essentially smooth.” Since $d({\bm{\alpha}})$ is a closely proper strictly convex function, its conjugate is smooth. Therefore, $f_{\mu}({\bm{\beta}})$ is a smooth function.

Now we apply Danskin’s theorem [Proposition B.25 in Bertsekas 1999] to derive $\nabla f_{\mu}({\bm{\beta}})$. Let $\phi({\bm{\alpha}},{\bm{\beta}})={\bm{\alpha}}^{T}C{\bm{\beta}}-\mu\,d({\bm{\alpha}})$. Since $d(\cdot)$ is a strongly convex function, $\mathop{\arg\max}_{{\bm{\alpha}}\in\mathcal{Q}}\phi({\bm{\alpha}},{\bm{\beta}})$ has a unique optimal solution and we denote it as ${\bm{\alpha}}^{\ast}$. According to Danskin’s theorem,

As for the proof of the Lipschitz constant of $f_{\mu}({\bm{\beta}})$, readers may refer to Nesterov 2005.

.2 Proof of Proposition 1

.3 Proof of Proposition 2

Similar to the proof technique of Proposition 1, we reformulate the problem of solving ${\bm{\alpha}}^{\ast}$ as a Euclidean projection:

According to the construction of the matrix $C$, we have, for any vector ${\mathbf{v}}$,

By the simple fact that $(a\pm b)^{2}\leq 2a^{2}+2b^{2}$ and the inequality holds as equality if and only if $a=\pm b$, for each edge $e=(m,l)\in E$, the value $(v_{m}-\operatorname{sign}(r_{ml})v_{l})^{2}$ is upper bounded by $2v_{m}^{2}+2v_{l}^{2}$. Hence, when $\|{\mathbf{v}}\|_{2}=1$, the right-hand side of (29) can be further bounded by

Note that this upper bound is tight because the first inequality in (.3) is tight.

.4 Proof of Theorem 2

Based on the result from Beck and Teboulle 2009, we have the following lemma:

For the function $\widetilde{f}({\bm{\beta}})=h({\bm{\beta}})+\lambda\|{\bm{\beta}}\|_{1}$, where $h({\bm{\beta}})$ is an arbitrary convex smooth function and its gradient $\nabla h({\bm{\beta}})$ is Lipschitz continuous with the Lipschitz constant $L$, we apply Algorithm 1 to minimize $\widetilde{f}({\bm{\beta}})$ and let ${\bm{\beta}}^{t}$ be the approximate solution at the $t$th iteration. For any ${\bm{\beta}}$, we have the following bound:

In order to use the bound in (31), we use the similar proof scheme as in Lan, Lu and Monteiro 2011 and decompose $f({\bm{\beta}}^{t})-f({\bm{\beta}}^{\ast})$ into three terms:

According to the definition of $\widetilde{f}$, we know that for any ${\bm{\beta}}$

where $D\equiv\max_{{\bm{\alpha}}\in\mathcal{Q}}d({\bm{\alpha}})$. Therefore, the first term in (.4), $f({\bm{\beta}}^{t})-\widetilde{f}({\bm{\beta}}^{t})$, is upper-bounded by $\mu D$, and the last term in (.4) is less than

or equal to 0 [i.e., $\widetilde{f}({\bm{\beta}}^{\ast})-f({\bm{\beta}}^{\ast})\leq 0$]. Combining (31) with these two simple bounds, we have

By setting $\mu=\frac{\varepsilon}{2D}$ and plugging this into the right-hand side of (.4), we obtain

If we require the right-hand side of (34) to be equal to $\varepsilon$ and solve it for $t$, we obtain the bound of $t$ in (21).

Acknowledgments

We would like to thank Yanjun Qi for the help of preparation and verification of breast cancer data, and Javier Peña for the discussion of the related first-order methods. We would also like to thank the anonymous reviewers and the Associate Editor for their constructive comments on improving the quality of the paper.

References