Convergence of Common Proximal Methods for L1-Regularized Least Squares

1 Introduction

Many problems in machine learning and data fitting can be cast as a least squares problem with a regularization term to limit overfitting. Using an ℓ1-norm regularization has been found to be particularly effective in many applications like feature selection [Tibshirani, 1996], compressed sensing [Chen et al., 1998], sparse coding [Gregor and LeCun, 2010], and discovery of graph connectivity [Hsieh et al., 2011]. In this paper, we use a model problem consisting of a linear least squares problem with ℓ1-regularization, also known as a LASSO problem:

where A ∈ Rm×n is a short-flat matrix (i.e. n > m) with full row rank, b is a given vector, and λ is a positive scalar. The ℓ1 regularizer tends to produce a sparse solution, avoiding overfitting while reducing the computational cost [Tibshirani, 1996].

Most general convex solvers such as interior point methods [Ben-Tal and Nemirovski, 2001] do not scale to the largescale data problems encountered in practice, but many algorithms have been proposed recently to take advantage of the special structure in (1). In this paper we compare four of the most popular of these methods: the Alternating Direction Method of Multipliers (ADMM) [Boyd et al., 2011], Iterative Shrinkage Thresholding Algorithm (ISTA) [Parikh and Boyd, 2014] and its accelerated version Fast ISTA (FISTA) [Beck and Teboulle, 2009], and the cyclic Coordinate Descent method (CD) [Saha and Tewari, 2010].

There has been a recent flurry of activity on bounds for the convergence rates of these methods. A global convergence rate bound of O(1/k) has been shown for ADMM [Eckstein and Bertsekas, 1990; Deng and Yin, 2012; He and Yuan, 2012] and ISTA, while FISTA enjoys a bound of O(1/k2 ) [He and Yuan, 2012]. Moreover, for ISTA, there have been recent local convergence results under some strict sparsity conditions, e.g. [Bredies and Lorenz, 2008; Liang et al., 2014]. The convergence of CD in general was shown in [Tseng and Yun, 2009], and an O(1/k) convergence rate for the model LASSO problem was shown in [Saha and Tewari, 2010]. In this paper, we focus on local convergence behavior, by modelling each iteration as a matrix recurrence. We find that each iteration passes through several phases or regimes before reaching a final regime with linear convergence, as first found in [Boley, 2013] for ADMM on general LPs and QPs. The analysis shows that ISTA can sometimes be faster than FISTA when close enough to the solution.

The paper is organized as follows. Section 2 gives some basic preliminaries of the LASSO problem. Sections 3 & 4 give the linear convergence analysis of different iterations. Numerical Examples are shown in Section 5 and conclusions are drawn in Section 6.

2 Preliminaries

2.1 Optimality condition of the LASSO problem

The first order KKT optimality conditions for the LASSO problem (1) are

where each component of ν satisfies

Here the "sign" function is defined by sign(x) = 1, 0,-1 according to whether x > 0, x = 0, x < 0, respectively.

2.2 Uniqueness

There are various sufficient and necessary conditions for the uniqueness of the LASSO problem or its variants. For example, [Osborne et al., 1999; Cand`es and Plan, 2009; Fuchs, 2005] showed different sufficient conditions and [Tibshirani, 2013] studied the necessary conditions for the LASSO problem. In fact, the problem (1) needs to have a unique solution in many situations. For example, in compressed sensing, having non-uniqueness solutions will result in unreliable recovery given the data. We refer readers to [Tibshirani, 2013; Zhang et al., 2012] and references therein for a discussion of the uniqueness of the LASSO problem. In some cases, we will assume uniqueness of the LASSO solution.

3 Convergence of ADMM, ISTA and FISTA

3.1 Auxiliary Variables with Local Monotonic Behavior

In this section, we show that ADMM, ISTA and FISTA can be all transformed into a matrix recurrence form in a similar way. We distinguish the iterates of the different algorithms using the notation x, ̂x and ˜x to denote the iterates of ADMM, ISTA and FISTA respectively.

ADMM as a Matrix Recurrence

The ADMM is constructed by splitting the primal x variables into two separate variables such that the minimum with respect to each individual variable can be easily computed, and then imposing an equality constraint between the two variables. A typical splitting for LASSO problem is to use variable x for the least squares loss function and z for the l1-norm regularizer. Then the modified LASSO problem becomes

with augmented Lagrangian function for the resulting problem being

where ρ is a penalty parameter, µ is the dual variable. Let u = µ/ρ, then the formal ADMM iterates can be represented

where Shrσ(s) = (1 - σ/|s|)+s and Thrσ(s) = s - Shrσ(s). We replace the iterates z[k] , u[k] with two equivalent "auxiliary" iterates carrying the same information. One variable, namely w[k] , exhibits smooth behavior, with linear convergence locally around a fixed point, and the other variable d[k] is a discrete ternary "flag" vector indicating which of the three cases of the shrinkage operator applies to each component. Specifically, for all k, we let the common iterate be w[k] = z[k] + u[k] and d[k] be a vector defined elementwise as d[k] i = sign(Shrλ/ρ (w[k] i )) and the flag matrix

D[k] = diag(d[k] ). Then one can derive

where

(7)

The ADMM update with x[k+1] , z[k+1] , u[k+1] now can be modified in matrix form in terms of D[k+1] and w[k+1] as below.

Note that step 1 of Alg. 1 is written as a homogeneous matrix recurrence, which will be used to characterize its convergence property.

where we denote M[k] aug as (M [k] h[k] 0 1 ), the augmented matrix of M [k] , in this paper.

FISTA and ISTA as a Matrix Recurrence

The ISTA and FISTA iteration [Daubechies et al., 2004; Beck and Teboulle, 2009] are given as follows, where t[0] = t[1] = 1 and L is the Lipschitz constant equal to ‖AT A‖2.

If τ [k] = 0 (or equivalently t[k] = 1) for all k, leading to v[k+1] = ˜x[k] , then this reduces to the ISTA updates, specifi-

Similar to ADMM, we use auxiliary variables ˜w[k] , ˜D[k] to replace variable ˜x[k] to carry the FISTA iterations. Set ˜w[k] = (I - 1/LATA)v[k] + 1/LATb and the vector ˜d[k] is defined elementwise as ˜ d[k] i = sign(Shrλ/L ( ˜ w[k] i )) and the flag matrix

˜D[k] = diag(˜d[k] ). The iteration can be described as

where we denote

Therefore, the FISTA update (9) can be written using the new auxiliary variables ˜w and ˜D as follows.

Algorithm 2: One pass of modified FISTA

Step 2 of above procedure can also be formulated as a homogeneous matrix recurrence analogous to (8) for ADMM with a larger (approximately double) dimension:

where we denote N [k] = (P [k] Q[k-1] I 0 ) and ˜h[k] =

(¯h[k] 0 ) and N[k] aug = (N [k] ˜h[k] 0 1 ) in the remainder of this paper. For ISTA, Alg. 2 reduces to

where R[k] = (I - 1/LATA)( ̂D[k] )2 and ̂h[k] = -(I - 1/LATA)λ/L ̂d[k] + 1/LATb. We denote R[k] aug as (R[k] ̂h[k] 0 1 ), the augmented matrix of R[k] , in this paper.

3.2 Spectral Properties and Four Regimes

It is seen that M[k] aug, R[k] aug and N[k] aug play key roles in the convergence. We summarize their respective properties.

Lemma 3.1 Suppose D[k-1] = D[k] = D[k+1] , (and same for ̂D[k] and ˜D[k] ). Then the iteration matricesM [k] (7), R[k] (13), and N [k] (11) & (12) (for ADMM, ISTA and FISTA, respectively) have the following properties:

• (a). ‖M [k] ‖2 ≤ 1, ‖R[k]‖2 ≤ 1, ‖N [k] ‖2 ≤ 1.

(b). All eigenvalues or M [k] ,R[k] ,N [k] lie in the closed disk in the complex plane with center 1/2 and radius 1/2, denoted as D(1/2, 1/2), and the eigenvalues of R[k] are real.

• (c). The eigenvalue 1, if it exists, must have a complete set of eigenvectors (no Jordan blocks larger than 1 × 1).

Proof. We temporarily omit the pass number [k] . For M (ADMM) we have:

(a). Observe 2D2 - I is an orthogonal matrix by (2D2 - I)(2D2 - I) = I. Hence, ‖M‖2 = ‖M(2D2 - I)‖2 = ‖D2 - I + (1/ρATA + I)-1‖2 ≤ 1.

eigenvalues of M lie in the disk D(1/2, 1/2), which is entirely in the open right half plane plus the origin.

The proof of (c) and all the cases for ISTA and FISTA follow the same lines as in [Boley, 2013; Tao et al., 2015] and hence are omitted. □

Lemma 3.1 gives rise to the four possible "regimes" associated with the ADMM, ISTA and FISTA iterations, depending on the flag matrix and the eigenvalues of operators Maug, Raug, Naug. We treat separately the case where the flag matrix remains the same at each iteration, in which there are three possible regimes, and treat all the transitional cases together in their own fourth regime.

When flag matrix remains unchanged from iteration k to k +1, (D[k] =D[k+1] (or ̂D[k] = ̂D[k+1] , ˜D[k] = ˜D[k+1] )):

Regime [A]. The spectral radius of M [k] (or R[k] , N [k] ) is strictly less than 1. If close enough to the optimal solution (if it exists), the result is linear convergence to that solution. The convergence rate depends on the second largest eigenvalue of Maug (or Raug, Naug, resp.) according to the theory of power method for the matrix eigenvalue problem [Golub and Loan, 2013].

Regime [B]. Also known as "constant step regime", M [k] (or R[k] , N [k] ) has an eigenvalue equal to 1 but which yields

mented matrix. Then the iteration process tends to a constant step (eigenvector of M [k] (or R[k] , N [k] ) for eigenvalue 1). The iteration will continue in this way until the discrete flag matrix D changes. Such a change is guaranteed to occur due to the global convergence of the algorithms.

Regime [C]. M [k] (or R[k] , N [k] ) has an eigenvalue equal to 1, but the augmented matrix still has a complete set of eigenvectors for eigenvalue 1 (this eigenvalue has no 2 × 2 Jordan block). The iterates will converge to an eigenvector for eigenvalue 1 with a linear rate as in Regime [A]. This cannot occur in the final regime if the original optimization problem has a unique solution.

When flag matrix changes at next iteration, i.e. D[k] ̸= D[k+1] , then we have:

Regime [D]. The operator M [k+1] (or R[k+1] , N [k+1] ) will be different from M [k] (or R[k] , N [k] ) due to different flag matrix.

3.3 Local Linear Convergence

Now we present part of our main results below. Essentially, when close enough to the optimal solution, we can give a guarantee that eventually the flag matrix will not change. Then the iterations of ADMM, ISTA and FISTA eventually behave like the power method for eigenvalues such that they converge linearly to the optimal solution. The proof of the following theorem follows [Boley, 2013] (for ADMM) and [Tao et al., 2015] (for F/ISTA).

Theorem 3.2 Suppose the LASSO problem (1) has a unique solution and this solution has strict complementarity: that is for every index i, w∗ i ̸= ±λ/ρ for ADMM, ̂w∗ i ̸= ±λ/L for ISTA, ˜w∗ i ̸= ±λ/L for FISTA. Then eventually the ADMM iteration used specifically for the LASSO problem (6), the ISTA

iteration (10), and the FISTA iteration (9) all reach a stage where the iteration converges linearly to that unique solution.

The eigenanalysis in terms of regimes established in [Tao et al., 2015] allows one not only to study the local convergence behavior but also to explain the whole iteration behavior. The following two propositions show that how FISTA can be viewed as an accelerated ISTA process and FISTA may slow down compared to ISTA when close enough to the optimal solution.

Proposition 3.3 In regime [B], the constant step vector for ISTA is v, where v = Rv is a scaled eigenvector of R for eigenvalue 1, while the constant step vector for FISTA is 1 1-τ [k] v. Since τ [k] → 1, the constant stepsize of FISTA is much larger, which yields a speedup.

Proposition 3.4 In regime [A] or [C], let βmax, γmax be the largest eigenvalue of R and N respectively, then FISTA is faster than ISTA if 1 > βmax > τ [k] > 0 but slower if 1 > τ [k] > βmax > 0. Since βmax is a fixed value for one specific regime, with the τ growing to 1, ISTA will be faster than FISTA toward the end. Besides, when τ > βmax, γmax must be one of a complex conjugate pair.

The proof follows the same lines as [Tao et al., 2015]. In practice, if βmax is well separated from 1, then it is advantageous to make FISTA iterations switch to ISTA once it reaches the final regime. This idea is implemented in Section 5 called Hybrid F/ISTA.

4 Convergence of CD

4.1 Local linear convergence

The cyclic CD is widely used due to its easy update rule. Essentially it goes through and updates all of the components in a cyclic fashion instead of updating them simultaneously as the gradient descent method. In this section, we establish a natural relationship between the iterations of the cyclic CD and of the Gauss-Seidel method when close enough to the optimal solution of problem (1) so that linear convergence is guaranteed eventually under some mild conditions.

We denote y as the iterates of CD, y[k] i as the i-th coordinate of y[k] and Ai as the i-th column of A. All coordinates other than i are denoted as -i. The CD updates for problem (1) is as follows.

Algorithm 3: One pass of CD

Start with y[k] .

for coordinate i = 1, 2, · · · , n,

end for

Result is y[k+1] to the next pass.

Our result is motivated by the following key observation.

Lemma 4.1 Suppose the solution to problem (1) has strict complementarity as in Theorem 3.2, then there exists a K such that for all k > K, the CD iterate y[k] is close enough to the optimal solution y∗ that sign(y[k] ) is fixed.

Proof. We first define the following index set based on y∗

and E as the complement set of E. Consequently,

Based on the notation, it also can be seen that AE is composed of the columns corresponding to the nonzero element of y∗E . For the simplicity of our analysis and without loss of generality, we split matrix A based on E and permute the columns so that A = [AE , AE ]. Then the optimality condition of problem (1) then can be rewritten as

where dE is composed of elements {±1} based on sign(y∗E ), and dE is composed of elements strictly between -1 and +1 by the assumption of strict complementarity. Let δ be the largest entry in absolute value found in dE so that

Now consider the case of E. Let y[k] be a vector very close to y∗ such that ‖y[k] - y∗‖∞ < ϵ1. Then |yl - y∗ l | < ϵ1 ∀l ∈ E and |yj| < ϵ1 ∀j ∈ E. If ϵ1 is sufficiently small, then (b) above induces

so that Shrλ(y[k] E ) = 0. Here c1 is some constant magnification factor depending only on A. Since y[k] converges to y∗ , it would imply that the future iterates, starting from y[k] , would have zeros in the E part.

Next consider the case of E. Let ϵ2 = minE{|y∗E |} - c2 for a positive constant c2 sufficiently small to make ϵ2 > 0. And for each l ∈ E, we define a ball around each y∗ l that B(y∗ l ) = {yl : ‖yl -y∗ l ‖∞ ≤ ϵ2}.

According to Theorem 16 in [Saha and Tewari, 2010], cyclic CD converges to LASSO problem at the rate of O(1/k). Hence there must exist an iteration number K that

which implies y[k] E = 0 and y[k] l ∈ B(y∗ l ), ∀l ∈ E. In other words, for each component i, y[k] i must fall in one of three cases: yi < 0, yi = 0 and yi > 0 and never jump out to another case. □

Theorem 4.2 Suppose the optimal solution y∗ of problem (1) is sparse such that ATE AE is positive definite, with E defined in (14). Then eventually the cyclic coordinate descent iteration reaches a stage where it converges linearly to the solution.

Proof. From Lemma 4.1, when close enough to the optimal solution, sign(y[k] ) is fixed. So y[k] E remained zero for future iterations. As for y[k] E , we can simplify the resulting CD

updates in Alg. 3 by eliminating y[k] j (∀j ∈ E) so that ∀l ∈ E,

Assuming λ is large enough such that optimal solution y∗ is sparse and ATE AE is positive definite, then (15) is equivalent to the Gauss-Seidel method applied to

where elements of dE is fixed and equal to either 1 or -1. The iteration must converge linearly by the theory of GaussSeidel method [Golub and Loan, 2013]. □

Remark: we note here that ATE AE being positive definite is a mild assumption in practice. Since the optimal solution is sparse, |E| < m can be satisfied easily for λ not too small, and in many applications this is sufficient to guarantee that all columns of AE are linearly independent. If not, one can increase λ to increase the sparsity.

4.2 Comparison with ISTA

In this part, we show that CD should converge faster than ISTA when both CD and ISTA iterations are in their final regimes from the viewpoint of preconditioning. The next lemma shows the equivalence of the ISTA iteration and the classical Richardson iteration [Kelley, 1995].

Lemma 4.3 When ISTA iteration reaches the final regime, the regime of linear convergence, then the ISTA iteration is equivalent to the Richardson iteration for solving the linear system (16).

Proof. When ISTA reaches the final regime, according to Theorem 3.2 ̂w = (I - 1/LATA)̂x + 1/LATb would fall in into three cases: ̂w < -L, -L < ̂w < L , ̂w > L , and never jump out to another case. Hence the shrinkage operator in updating step (10) is fixed so that ISTA reduces to

where ̂dE = sign(x∗E ). And ̂x[k+1] E = 0. The resulting ISTA updates (17) is exactly Richardson iteration [Kelley, 1995] for solving (16). □

Combining the result of Theorem 4.2 and Lemma 4.3, if both CD and ISTA reach their final regimes, then CD is equivalent to Gauss-Seidel iteration and ISTA is equivalent to Richardson iteration for solving the same linear system (16). Let T and U be the diagonal and strict upper triangular part of ATE AE and L is the lipschitz constant. Gauss-Seidel iteration can be written as the preconditioned Richardson iteration as below (and hence generally faster):

with preconditioner L(T +U T )-1 [Kelley, 1995].

5 Numerical Examples

We consider examples of compressed sensing to show different convergence behaviors to support our analysis. Suppose there exists a true sparse signal represented by a ndimensional vector xs with s non-zero elements. We observe the image of xs under the linear transformation Axs, where A is the so-called measurement matrix. Our observation thus should be

where ϵ is the observation noise. The goal is to recover the sparse vector xs from the measurement matrix A and observation b. This is a form of feature selection, where we discover the most relevant features. We let A ∈ Rm×n be Gaussian matrix whose elements are i.i.d distributed as N(0, 1), ϵ be a vector whose elements are i.i.d distributed as N(0, σ2 ) with σ = 10-3 .

The numerical results are summarized in Table 1. We compare the convergence behavior of ADMM, ISTA, FISTA and CD in terms of the total number of iterations (Total #), number of iterations to reach final linear regime (Final #) and the linear rate in the final regime (Rate), with different λ and sparsity s. Based on Propositions 3.3 & 3.4, we also report the behavior of a hybrid F/ISTA that follows FISTA during initial iterations to reach the final regime and then switch to ISTA until convergence. In general, after comparison with the linear rates, we see that, CD is often much faster than the other iterations during the final linear regime.

Fig. 1 illustrates the methods' behavior for the instance marked ∗∗ in Table 1, to show how the theorems established before explains the behaviors in practice. Fig. 1(Left) shows the convergence behavior of all algorithms. Observe that all of them pass through a few transitions in the early part of the iterations and then settle on their final linear regimes. Fig. 1(Right) shows the eigenvalues during the final regimes for all algorithms. One can see that the eigenvalues of operators Maug,Raug,Naug lie strictly inside the circle D(1/2, 1/2) (except for 0 & 1), consistent with Lemma 3.1. Based on this instance, we also observe the following:

(a). It costs FISTA many fewer steps (245 iterations) than ISTA (2823 iterations) to get to the final regime. The main reason is that FISTA has much larger constant steps in regime

[B] so that it can more quickly reach the end of the stagnating regime [B], as suggested by Proposition 3.3. In fact, one can show that the difference between consecutive iterates of ISTA remains a constant for many iterations while FISTA does not. Moreover, since the rate of ISTA is 0.959, well separated from 1, Proposition 3.4 predicts switching to ISTA in the final regime should converge faster than standard FISTA. Indeed, hybrid F/ISTA converges in 411 iterations compared to 502 iterations for FISTA.

(b). FISTA oscillates in the final regime. This is because the second largest eigenvalue of FISTA operator Naug is a pair of complex conjugates (cf. Lemma 3.1 & Proposition 3.4). Hence, according to the theory of the power method, the convergence will oscillate between the two conjugate complex numbers. Since all the eigenvalues of ISTA operator Raug are real, ISTA iterates do not oscillate.

(c). Implied by the observation in Section 4.2, in the final regime, CD is equivalent to a preconditioned ISTA-like Richardson method for solving the same linear system. In this instance, we can see the preconditioning plays a big role and the eigenvalues of CD are much smaller than ISTA shown in Fig. 1(Right).

(d). We remark that the cost per iteration is about the same (O(n2 )) for all the methods. ISTA, FISTA, and CD require the equivalent of a matrix-vector product. ADMM requires an LU-factorization of (ATA+ ρI), but only once at the beginning, and requires only a forward-back substitution in each iteration. Hence iteration counts closely reflect total times.

6 Conclusion

In this paper, we show the locally linear convergence of ADMM, ISTA, FISTA and CD applied to the LASSO problem. We model ADMM, ISTA and FISTA as the matrix recurrence form and connect them with the power method. We also establish a connection between CD method and GaussSeidel method. By spectral analysis, we show that all of the iterations normally pass through several regimes of different types and eventually settle on a "linear regime" in which the iterates converge linearly.

Such analysis provides a way to study the behavior through the whole iteration process. We explain why ISTA often appears to stagnate during the initial iterations, and why FISTA oscillates towards the end. Besides, we illustrate how the unaccelerated ISTA can sometimes be faster when close enough to the solution compared to FISTA and propose the idea of switching to ISTA in certain circumstances.

Acknowledgements

This research was supported in part by NSF grant IIS1319749.