Accelerated Linearized Bregman Method

Introduction

In this paper, we are interested in the following optimization problem

Since the development of the new paradigm of compressed sensing , the basis pursuit problem (1.2) has become a topic of great interest. In compressed sensing, $A$ is usually the product of a sensing matrix $\Phi$ and a transform basis matrix $\Psi$ and $b$ is a vector of the measurements of the signal $s=\Psi x$. The theory of compressed sensing guarantees that the sparsest solution (i.e., representation of the signal $s=\Psi x$ in terms of the basis $\Psi$) of $Ax=b$ can be obtained by solving (1.2) under certain conditions on the matrix $\Phi$ and the sparsity of $x$. This means that (1.2) gives the optimal solution of the following NP-hard problem :

where $\|x\|_{0}$ counts the number of nonzero elements of $x$.

Matrix generalizations of (1.3) and (1.2), respectively, are the so-called matrix rank minimization problem

and its convex relaxation, the nuclear norm minimization problem:

The matrix completion problem has a lot of interesting applications in online recommendation systems, collaborative filtering , etc., including the famous Netflix problem . It has been proved that under certain conditions, the solutions of the NP-hard problems (1.4) and (1.6) are given respectively by solving their convex relaxations (1.5) and (1.7), with high probability (see, e.g., ).

The linearized Bregman (LB) method was proposed in to solve the basis pursuit problem (1.2). The method was derived by linearizing the quadratic penalty term in the augmented Lagrangian function that is minimized on each iteration of the so-called Bregman method introduced in while adding a prox term to it. The linearized Bregman method was further analyzed in and applied to solve the matrix completion problem (1.7) in .

Throughout of this paper, we will sometimes focus our analysis on the basis pursuit problem (1.2). However, all of the analysis and results can be easily extended to (1.5) and (1.7). The linearized Bregman method depends on a single parameter $\mu>0$ and, as the analysis in shows, actually solves the problem

rather than the problem (1.2). Recently it was shown in that the solution to (1.8) is also a solution to problem (1.2) as long as $\mu$ is chosen large enough. Furthermore, it was shown in that the linearized Bregman method can be viewed as a gradient descent method applied to the Lagrangian dual of problem (1.8). This dual problem is an unconstrained optimization problem of the form

where the objective function $G_{\mu}(y)$ is differentiable since $g_{\mu}(x)$ is strictly convex (see, e.g., ). Motivated by this result, some techniques for speeding up the classical gradient descent method applied to this dual problem such as taking Barzilai-Borwein (BB) steps , and incorporating it into a limited memory BFGS (L-BFGS) method , were proposed in . Numerical results on the basis pursuit problem (1.2) reported in show that the performance of the linearized Bregman method can be greatly improved by using these techniques.

Our starting point is also motivated by the equivalence between applying the linearized Bregman method to (1.2) and solving the Lagrangian dual problem (1.9) by the gradient descent method. Since the gradient of $G_{\mu}(y)$ can be shown to be Lipschitz continuous, it is well-known that the classical gradient descent method with a properly chosen step size will obtain an $\epsilon$-optimal solution to (1.9) (i.e., an approximate solution $y^{k}$ such that $G_{\mu}(y^{k})-G_{\mu}(y^{*})\leq\epsilon$) in $O(1/\epsilon)$ iterations. In , Nesterov proposed a technique for accelerating the gradient descent method for solving problem of the form (1.9) (see, also, ), and proved that using this accelerated method, the number of iterations needed to obtain an $\epsilon$-optimal solution is reduced to $O(1/\sqrt{\epsilon})$ with a negligible change in the work required at each iteration. Nesterov also proved that the $O(1/\sqrt{\epsilon})$ complexity bound is the best bound that one can get if one uses only the first-order information. Based on the above discussion, we propose an accelerated linearized Bregman (ALB) method for solving (1.8) which is equivalent to an accelerated gradient descent method for solving the Lagrangian dual (1.9) of (1.8). As a by-product, we show that the basic and the accelerated linearized Bregman methods require $O(1/\epsilon)$ and $O(1/\sqrt{\epsilon})$ iterations, respectively, to obtain an $\epsilon$-optimal solution with respect to the Lagrangian for (1.8).

The rest of this paper is organized as follows. In Section 2 we describe the original Bregman iterative method, as well as the linearized Bregman method. We motivate the methods and state some previously obtained theoretical results that establish the equivalence between the LB method and a gradient descent method for the dual of problem (1.8). We present our accelerated linearized Bregman method in Section 3. We also provide a theoretical foundation for the accelerated algorithm and prove complexity results for it and the unaccelerated method. In Section 4, we describe how the LB and ALB methods can be extended to basis pursuit problems that include additional convex constraints. In Section 5, we report preliminary numerical results, on several compressed sensing basis pursuit and matrix completion problems. These numerical results show that our accelerated linearized Bregman method significantly outperforms the basic linearized Bregman method. We make some conclusions in Section 6.

Bregman and Linearized Bregman Methods

The Bregman method was introduced to the image processing community by Osher et al.in for solving the total-variation (TV) based image restoration problems. The Bregman distance with respect to convex function $J(\cdot)$ between points $u$ and $v$ is defined as

where $p\in\partial J(v)$, the subdifferential of $J$ at $v$. The Bregman method for solving (1.1) is given below as Algorithm 1. Note that the updating formula for $p^{k}$ (Step 4 in Algorithm 1) is based on the optimality conditions of Step 3 in Algorithm 1:

It was shown in that the Bregman method (Algorithm 1) converges to a solution of (1.1) in a finite number of steps.

It is worth noting that for solving (1.1), the Bregman method is equivalent to the augmented Lagrangian method in the following sense.

The sequences $\{x^{k}\}$ generated by Algorithm 1 and by the augmented Lagrangian method, which computes for $k=0,1,\cdots$

starting from $\lambda^{0}=0$ are exactly the same.

From Step 4 of Algorithm 1 and the fact that $p^{0}=0$, it follows that $p^{k}=-\sum_{j=1}^{k}A^{\top}(Ax^{j}-b)$. From the second equation in (1) and using $\lambda^{0}=0$, we get $\lambda^{k}=-\sum_{j=1}^{k}(Ax^{j}-b)$. Thus, $p^{k}=A^{\top}\lambda^{k}$ for all $k$. Hence it is easy to see that Step 3 of Algorithm 1 is exactly the same as the first equation in (1) and that the $x^{k+1}$ computed in Algorithm 1 and (1) are exactly the same. Therefore, the sequences $\{x^{k}\}$ generated by both algorithms are exactly the same. ∎

Although there are many algorithms for solving the subproblem (2.5) such as FPC , SPGL1 , FISTA etc., it often takes them many iterations to do so. The linearized Bregman method was proposed in , and used in to overcome this difficulty. The linearized Bregman method replaces the quadratic term $\frac{1}{2}\|Ax-b\|^{2}$ in the objective function that is minimized in Step 3 of Algorithm 1 by its linearization $\langle A^{\top}(Ax^{k}-b),x\rangle$ plus a proximal term $\frac{1}{2\mu}\|x-x^{k}\|^{2}$. Consequently the updating formula for $p^{k}$ is changed since the optimality conditions for this minimization step become:

In Algorithm 2 below we present a slightly generalized version of the original linearized Bregman method that includes an additional parameter $\tau$ that corresponds to the length of a gradient step in a dual problem.

In , it is shown that when $\mu\|A\|^{2}<2$, where $\|A\|$ denotes the largest singular value of $A$, the iterates of the linearized Bregman method (Algorithm 2 with $\tau=1$) converge to the solution of the following regularized version of problem (1.1):

We prove in Theorem 3 below an analogous result for Algorithm 2 for a range of values of $\tau$. However, we first prove, as in , that the linearized Bregman method (Algorithm 2) is equivalent to a gradient descent method

of (2.6), which we express as the following equivalent minimization problem:

To show that $G_{\mu}(y)$ is continuously differentiable, we rewrite $G_{\mu}(y)$ as

is strictly convex and continuously differentiable with gradient $\nabla\Phi_{\mu}(v)=\frac{v-\hat{w}}{\mu}$, and $\hat{w}=\arg\min_{w}\{J(w)+\frac{1}{2\mu}\|w-v\|^{2}\}$ (e.g., see Proposition 4.1 in ). From this it follows that $\nabla G_{\mu}(y)=Aw^{*}-b$. Hence the gradient method (2.7) corresponds to Algorithm 3 below.

Lemma 2 and Theorem 3 below generalize Theorem 2.1 in by allowing a step length choice in the gradient step (2.7) and show that Algorithms 2 and 3 are equivalent. Our proof closely follows the proof of Theorem 2.1 in .

$x^{k+1}$ computed by Algorithm 2 equals $w^{k+1}$ computed by Algorithm 3 if and only if

By comparing Step 3 in Algorithms 2 and 3, it is obvious that $w^{k+1}$ is equal to $x^{k+1}$ if and only if (2.9) holds. ∎

The sequences $\{x^{k}\}$ and $\{w^{k}\}$ generated by Algorithms 2 and 3 are the same.

We prove by induction that equation (2.9) holds for all $k\geq 0$. Note that (2.9) holds for $k=0$ since $p^{0}=x^{0}=0$ and $y^{0}=\tau b$. Now let us assume that (2.9) holds for all $0\leq k\leq n-1$; thus by Lemma 2 $w^{k+1}=x^{k+1}$ for all $0\leq k\leq n-1$. By iterating Step 4 in Algorithm 3 we get

By iterating Step 4 in Algorithm 2 we get

where the last equality follows from (2.10); thus by induction (2.9) holds for all $k\geq 0$, which implies by Lemma 2 that $x^{k}=w^{k}$ for all $k\geq 0$. ∎

Before analyzing Algorithms 2 and 3, we note that by defining $v^{k}=A^{\top}y^{k}$ and algebraically manipulating the last two terms in the objective function in Step 3 in Algorithm 3, Steps 3 and 4 in that algorithm can be replaced by

if we set $v^{0}=\tau A^{\top}b$. Because Algorithms 2 and 3 are equivalent, convergence results for the gradient descent method can be applied to both of them. Thus we have the following convergence result.

Let $J(w)\equiv\|w\|_{1}$. Then $G_{\mu}(y)$ in the dual problem (2.8) is continuously differentiable and its gradient is Lipschitz continuous with the Lipschitz constant $L\leq\mu\|A\|^{2}$. Consequently, if the step length $\tau<\frac{2}{\mu\|A\|^{2}}$, the sequences $\{x^{k}\}$ and $\{w^{k}\}$ generated by Algorithms 2 and 3 converge to the optimal solution of (2.6).

When $J(x)=\|x\|_{1}$, $w^{k+1}$ in (2) reduces to

$G_{\mu}(y)$ is continuously differentiable since $g_{\mu}(x)$ is strictly convex. Since for any point $y$, $\nabla G_{\mu}(y)=Aw-b$, where $w=\mu\cdot\mbox{shrink}(A^{\top}y,1)$, it follows from the fact that the shrinkage operator is non-expansive, i.e.,

for any two points $y^{1}$ and $y^{2}$. Thus the Lipschitz constant $L$ of $\nabla G_{\mu}(\cdot)$ is bounded above by $\mu\|A\|^{2}$.

When $\tau<\frac{2}{\mu\|A\|^{2}}$, we have $\tau L<2$ and thus $|1-\tau L|<1$. It then follows that the gradient descent method $y^{k+1}=y^{k}-\tau\nabla G_{\mu}(y^{k})$ converges and therefore Algorithms 2 and 3 converge to $x_{\mu}^{*}$, the optimal solution of (2.6). ∎

Before developing an accelerated version of the LB algorithm in the next section. We would like to comment on the similarities and differences between the LB method and Nesterov’s composite gradient method and the ISTA method applied to problem (1.1) and related problems. The latter algorithms iterate Step 3 in the LB method (Algorithm 2) with $p^{k}=0$, and never compute or update the subgradient vector $p^{k}$. More importantly, their methods solve the unconstrained problem

Hence, while these methods and the LB method both linearize the quadratic term $\|Ax-b\|^{2}$ while handling the nonsmooth term $\|x\|_{1}$ directly, they are very different.

Similar remarks apply to the accelerated LB method presented in the next section and fast versions of ISTA and Nesterov’s composite gradient method.

The Accelerated Linearized Bregman Algorithm

Based on Theorem 3, i.e., the equivalence between the linearized Bregman method and the gradient descent method, we can accelerate the linearized Bregman method by techniques used to accelerate the classical gradient descent method. In , Yin considered several techniques such as line search, BB step and L-BFGS, to accelerate the linearized Bregman method. Here we consider the acceleration technique proposed by Nesterov in . This technique accelerates the classical gradient descent method in the sense that it reduces the iteration complexity significantly without increasing the per-iteration computational effort. For the unconstrained minimization problem (1.9), Nesterov’s accelerated gradient method replaces the gradient descent method (2.7) by the following iterative scheme:

where the scalars $\alpha_{k}$ are specially chosen weighting parameters. A typical choice for $\alpha_{k}$ is $\alpha_{k}=\frac{3}{k+2}$. If $\tau$ is chosen so that $\tau\leq 1/L$, where $L$ is the Lipschitz constant for $\nabla G_{\mu}(\cdot)$, Nesterov’s accelerated gradient method (3) obtains an $\epsilon$-optimal solution of (1.9) in $O(1/\sqrt{\epsilon})$ iterations, while the classical gradient method (2.7) takes $O(1/\epsilon)$ iterations. Moreover, the per-iteration complexities of (2.7) and (3) are almost the same since computing the gradient $\nabla G_{\mu}(\cdot)$ usually dominates the computational cost in each iteration. Nesterov’s acceleration technique has been studied and extended by many others for nonsmooth minimization problems and variational inequalities, e.g., see .

In the following, we first establish the equivalence between the accelerated linearized Bregman method and the corresponding accelerated gradient descent method (3), which we give explicitly as (5) below applied to the dual problem (2.8). Based on this, we then present complexity results for both basic and accelerated linearized Bregman methods. Not surprisingly, the accelerated linearized Bregman method improves the iteration complexity from $O(1/\epsilon)$ to $O(1/\sqrt{\epsilon})$.

More specifically, the sequence $\{x^{k}\}$ generated by Algorithm 4 is exactly the same as the sequence $\{w^{k}\}$ generated by (5).

Note that the Step 3 of Algorithm 4 is equivalent to

Comparing (3.8) with the first equation in (5), it is easy to see that $x^{k+1}=w^{k+1}$ if and only if

Thus we proved that (3.9) holds for $k=1$. Now let us assume that (3.9) holds for $0\leq k\leq n-1$, which implies $x^{k}=w^{k},\forall 0\leq k\leq n$ since $x^{0}=w^{0}$. We will prove that (3.9) holds for $k=n$.

where the first equality is from Step 4 of Algorithm 4 and the second equality is from (3.9) for $k=n-1$. From Step 6 of Algorithm 4 and (3.12), we have

where the last equality uses Step 5 of Algorithm 4. On the other hand, from (5) we have

where the third equality is from $w^{n}=x^{n}$ and $w^{n-1}=x^{n-1}$, the last equality is from Step 5 of Algorithm 4. Combining (3) and (3) we get that (3.9) holds for $k=n$. ∎

Like the linearized Bregman, we can also use a simpler implementation for accelerated linearized Bregman method in which the main computation at each step is a proximal minimization. Specifically, (5) is equivalent to the following three steps.

Next we prove iteration complexity bounds for both basic and accelerated linearized Bregman algorithms. Since these algorithms are standard gradient descent methods applied to the Lagrangian dual function and these results have been well established, our proofs will be quite brief.

Let the sequence $\{x^{k}\}$ be generated by the linearized Bregman method (Algorithm 2) and $(x^{*},y^{*})$ be the pair of optimal primal and dual solutions for Problem (2.6). Let $\{y^{k}\}$ be the sequence generated by Algorithm 3 and suppose the step length $\tau\leq\frac{1}{L}$, where $L$ is the Lipschitz constant for $\nabla G_{\mu}(y)$. Then for the Lagrangian function

Thus, if we further have $\tau\geq\beta/L$, where $0<\beta\leq 1$, then $(x^{k+1},y^{k})$ is an $\epsilon$-optimal solution to Problem (2.6) with respect to the Lagrangian function if $k\geq\lceil C/\epsilon\rceil$, where $C:=\frac{L\|y^{*}-y^{0}\|^{2}}{2\beta}$.

By using the convexity of function $G_{\mu}(\cdot)$ and the Lipschitz continuity of the gradient $\nabla G_{\mu}(\cdot)$, we get for any $y$,

Setting $y=y^{k-1}$ in (3), we obtain $G_{\mu}(y^{k})\leq G_{\mu}(y^{k-1})$ and thus the sequence $\{G_{\mu}(y^{k})\}$ is non-increasing. Moreover, summing (3) over $k=1,2,\ldots,n$ with $y=y^{*}$ yields

Before we analyze the iteration complexity of the accelerated linearized Bregman method, we introduce a lemma from that we will use in our analysis.

The following theorem gives an iteration-complexity result for the accelerated linearized Bregman method. Our proof of this theorem closely follows the proof of Proposition 2 in .

Let the sequence $\{x^{k}\}$ be generated by accelerated linearized Bregman method (Algorithm 4) and $(x^{*},y^{*})$ be the optimal primal and dual variable for Problem (2.6). Let $\{\alpha_{k}\}$ be chosen as

Let the sequence $\{y^{k}\}$ be defined as in (5) and the step length $\tau\leq\frac{1}{L}$, where $L$ is the Lipschitz constant of $\nabla G_{\mu}(y)$ and $G_{\mu}(\cdot)$ is defined by (3.28). We have

Thus, if we further have $\tau\geq\beta/L$, where $0<\beta\leq 1$, then $(x^{k+1},y^{k})$ is an $\epsilon$-optimal solution to Problem (2.6) with respect to the Lagrangian function (3.26) if $k\geq\lceil\sqrt{C/\epsilon}\rceil$, where $C:=\frac{2L\|y^{*}-y^{0}\|^{2}}{\beta}$.

and denote the linearization of $G_{\mu}(y)$ as

Therefore the second equality in (5) is equivalent to

Define $\hat{y}^{k}:=(1-\theta_{k})y^{k}+\theta_{k}y^{*}$, we have

From (3.37), it is easy to show that $\frac{1-\theta_{k}}{\theta_{k}^{2}}\leq\frac{1}{\theta_{k-1}^{2}}$ for all $k\geq 0$. Thus (3.41) implies that

Summing (3.42) over $k=0,1,\ldots,n-1$, we get

The proof technique and the choice of $\theta_{k}$ used here are suggested in for accelerating the basic algorithm. Other choices of $\theta_{k}$ can be found in . They all work here and give the same order of iteration complexity.

Extension to Problems with Additional Convex Constraints

We now consider extensions of both the LB and ALB methods to problems of the form

to compute a subgradient $p^{k+1}\in\partial J(x^{k+1})$. Fortunately, the Lagrangian dual gradient versions of these algorithms do not suffer from this difficulty. All that is required to extend them to problem (4.1) is to include the constraint $w\in X$ in the minimization step in these algorithms. Note that the gradient of

remains the same. Also it is clear that the iteration complexity results given in Theorems 6 and 8 apply to these algorithms as well.

Being able to apply the LB and ALB methods to problems of the form of (4.1) greatly expands their usefulness. One immediate extension is to compressed sensing problems in which the signal is required to have nonnegative components. Also (4.1) directly includes all linear programs. Applying the LB and ALB to such problems, with the goal of only obtaining approximated optimal solutions, will be the subject of a future paper.

Numerical Experiments

In this section, we report some numerical results that demonstrate the effectiveness of the accelerated linearized Bregman algorithm. All numerical experiments were run in MATLAB 7.3.0 on a Dell Precision 670 workstation with an Intel Xeon(TM) 3.4GHZ CPU and 6GB of RAM.

In this subsection, we compare the performance of the accelerated linearized Bregman method against the performance of the basic linearized Bregman method on a variety of compressed sensing problems of the form (1.2).

For compressed sensing problems, where $J(x)=\|x\|_{1}$, the linearized Bregman method reduces to the two-line algorithm:

Both algorithms are very simple to program and involve only one $Ax$ and one $A^{\top}y$ matrix-vector multiplication in each iteration.

We ran both LB and ALB with the $seed$ used for generating random number in MATLAB setting as $0$. Here we set $n=2000,m=0.4\times n,s=0.2\times m,\mu=5$ for all data sets. We set $\tau=\frac{2}{\mu\|A\|^{2}}$. We terminated the algorithms when the stopping criterion

was satisfied or the number of iterations exceeded 5000. Note that (5.3) was also used in . We report the results in Table 1.

In Table 1, we see that for three out of six problems, LB did not achieve the desired convergence criterion within 5000 iterations, while ALB satisfied this stopping criterion in less than 330 iterations on all six problems. To further demonstrate the significant improvement the ALB achieved over LB, we plot in Figures 3, 3 and 3 the Euclidean norms of the residuals and the relative errors as a function of the iteration number that were obtained by LB and ALB applied to the same data sets. These figures also depict the non-monotonic behavior of the ALB method.

2 Numerical Results on Matrix Completion Problems

There are fast implementations of linearized Bregman and other solvers for solving matrix completion problems. We do not compare the linearized Bregman and our accelerated linearized Bregman algorithms with these fast solvers here. Rather our tests are focused only on comparing ALB with LB and verifying that the acceleration actually occurs in practice for matrix completion problems.

The nuclear norm matrix completion problem (1.7) can be rewritten as

where $[\mathcal{P}_{\Omega}(X)]_{ij}=X_{ij}$ if $(i,j)\in\Omega$ and $[\mathcal{P}_{\Omega}(X)]_{ij}=0$ otherwise. When the convex function $J(\cdot)$ is the nuclear norm of matrix $X$, the Step 3 of Algorithm 2 with inputs $X^{k},P^{k}$ can be reduced to

It is known (see, e.g., ) that (5.5) has the closed-form solution,

where the matrix shrinkage operator is defined as

and $Y=U\mbox{Diag}(\sigma)V^{\top}$ is the singular value decomposition (SVD) of matrix $Y$. Thus, a typical iteration of the linearized Bregman method (Algorithm 2), with initial inputs $X^{0}=P^{0}=0$, for solving the matrix completion problem (5.4) can be summarized as

where the sequence $\alpha_{k}$ is chosen according to Theorem 8.

Note that this stopping criterion was used in . We also set the maximum number of iteration to 2000.

We report the number of iterations needed by LB and ALB to reach (5.14) in Table 2. Note that performing the shrinkage operation, i.e., computing an SVD, dominates the computational cost in each iteration of LB and ALB. Thus, the per-iteration complexities of LB and ALB are almost the same and it is reasonable to compare the number of iterations needed to reach the stopping criterion. We report the relative error $err:=\|X^{k}-M\|_{F}/\|M\|_{F}$ between the recovered matrix $X^{k}$ and the true matrix $M$ in Table 2. We see from Table 2 that ALB needed significantly fewer iterations to meet the stopping criterion (5.14).

In Figures 4 and 5, we plot the Frobenius norms of the residuals and the relative errors obtained by LB and ALB for iteration 1-500 for the tests involving matrices with dimension $n=200,300,400$ and $500$. Note that the non-monotonicity of ALB is far less pronounced on these problems.

Conclusions

In this paper, we analyzed for the first time the iteration complexity of the linearized Bregman method. Specifically, we show that for a suitably chosen step length, the method achieves a value of the Lagrangian of a quadratically regularized version of the basis pursuit problem that is within $\epsilon$ of the optimal value in $O(1/\epsilon)$ iterations. We also derive an accelerated version of the linearized Bregman method whose iteration complexity is reduced to $O(1/\sqrt{\epsilon})$ and present numerical results on basis pursuit and matrix completion problems that illustrate this speed-up.

References