Fast Linearized Bregman Iteration for Compressive Sensing and Sparse Denoising

Introduction

Let $A\in R^{m\times n}$, with $n>m$ and $f\in R^{m}$, be given. The aim of a basis pursuit problem is to find $u\in R^{n}$ by solving the constrained minimization problem:

where $J(u)$ is a continuous convex function.

Basis pursuit arises from many applications. In particular, there has been a recent burst of research in compressive sensing, which involves solving (1.1), (1.2). This was led by Candes et.al. , Donoho, , and others, see , and for extensive references. Compressive sensing guarantees, under appropriate circumstances, that the solution to (1.1), (1.2) gives the sparsest solution satisfying $Au=f$. The problem then becomes one of solving (1.1), (1.2) fast. Conventional linear programming solvers are not tailored for the large scale dense matrices $A$ and the sparse solutions $u$ that arise here. To overcome this, a linearized Bregman iterative procedure was proposed in and analyzed in , and . In , true, nonlinear Bregman iteration was also used quite successfully for this problem.

Bregman iteration applied to (1.1), (1.2) involves solving the constrained optimization problem through solving a small number of unconstrained optimization problems:

In we used a method called the fast fixed point continuation solver (FPC) which appears to be efficient. Other solvers of (1.3) could be used in this Bregman iterative regularization procedure.

Here we will improve and analyze a linearized Bregman iterative regularization procedure, which, in its original incarnation, , , involved only a two line code and simple operations and was already extremely fast and accurate.

In addition, we are interested in the denoising properties of Bregman iterative regularization, for signals, not images. The results for images involved the BV norm, which we may discretize for $n\times n$ pixel images as:

We usually regard the success of the ROF TV based model

The paper is organized as follows: In section 2 we describe Bregman iterative algorithms, as well as the linearized version. We motivate these methods and describe previously obtained theoretical results. In section 3 we introduce an improvement to the linearized version, call “kicking” which greatly speeds up the method, especially for solutions $u$ with a large dynamic range. In section 4 we present numerical results, including sparse recovery for $u$ having large dynamic range, and the recovery of signals in large amounts of noise. In another work in progress we apply these ideas to denoising very blurry and noisy signals remarkably well including sparse recovery for $u$. By blurry we mean situations where $A$ is perhaps a subsampled discrete convolution matrix whose elements decay to zero with $n$, e.g. random rows of a discrete Gaussian.

Bregman and Linearized Bregman Iterative Algorithms

The Bregman distance , based on the convex function $J$, between points $u$ and $v$, is defined by

where $p\in\partial J(v)$ is an element in the subgradient of $J$ at the point $v$. In general $D_{J}^{p}(u,v)\not=D_{J}^{p}(v,u)$ and the triangle inequality is not satisfied, so $D_{J}^{p}(u,v)$ is not a distance in the usual sense. However it does measure the closeness between $u$ and $v$ in the sense that $D_{J}^{p}(u,v)\geq 0$ and $D_{J}^{p}(u,v)\geq D_{J}^{p}(w,v)$ for all points $w$ on the line segment connecting $u$ and $v$. Moreover, if $J$ is convex, $D_{J}^{p}(u,v)\geq 0$, if $J$ is strictly convex $D_{J}^{p}(u,v)>0$ for $u\not=v$ and if $J$ is strongly convex, then there exists a constant $a>0$ such that

To solve (1.1) Bregman iteration was proposed in . Given $u^{0}=p^{0}=0$, we define:

It was proven in that if $J(u)\in C^{2}(\Omega)$ and is strictly convex in $\Omega$, then $\|Au^{k}-f\|$ decays exponentially whenever $u^{k}\in\Omega$ for all $k$. Furthermore, when $u^{k}$ converges, its limit is a solution of (1.1). It was also proven in that when $J(u)=|u|_{1}$, i.e. for problem (1.1) and (1.2), or when $J$ is a convex function satisfying some additional conditions, the iteration (2.7) leads to a solution of (1.1) in finitely many steps.

As shown in , see also , , the Bregman iteration (2.7) can be written as:

This was referred to as “adding back the residual” in . Here $f^{0}=0,u^{0}=0$. Thus the Bregman iteration uses solutions of the unconstrained problem

as a solver in which the Bregman iteration applies this process iteratively.

Since there is generally no explicit expression for the solver of (2.7) or (2.8), we turn to iterative methods. The linearized Bregman iteration which we will analyze, improve and use here is generated by

In the special case considered here, where $J(u)=\mu\|u\|_{1}$, then we have the two line algorithm

This linearized Bregman iterative algorithm was invented in and used and analyzed in , and . In fact it comes from the inner-outer iteration for (2.7). In it was shown that the linearized Bregman iteration (2.10) is just one step of the inner iteration for each outer iteration. Here we repeat the arguments also in , which begin by summing the second equation in (2.10) arriving at (using the fact that $u^{0}=p^{0}=0$):

This gives us (2.12), and allows us to rewrite its first equation as:

i.e. we are adding back the “linearized noise”, where $v^{k+1}$ is defined in (2.11).

In and some interesting analysis was done for (2.10), (and some for (2.14)). This was done first for $J(u)$ continuously differentiable in (2.10) and the gradient $\partial J(u)$ satisfying

$\forall u,v\in R^{n},\ \beta>0$. In it was shown that, if (2.15) is true, then both of the sequences $(u^{k})_{k\in N}$ and $(p^{k})_{k\in N}$ defined by (2.10) converge for $0<\delta<\frac{2}{\|AA^{T}\|}$.

In the authors recently give a theoretical analysis, showing that the iteration in (2.11) and (2.12) converges to the unique solution of

They also show the interesting result: let $S$ be the set of all solutions of the Basis Pursuit problem (1.1), (1.2) and let

which is unique. Denote the solution of (2.16) as $u_{\mu}^{*}$. Then

Another theoretical analysis on Linearized Bregman algorithm is given by Yin in , where he shows that Linearized Bregman iteration is equivalent to gradient descent applied to the dual of the problem (2.16) and uses this fact to obtain an elegant convergence proof.

This summarizes the relevant convergence analysis for our Bregman and linearized Bregman models.

Next we recall some results from regarding noise and Bregman iteration.

i.e., until we get too close to the noisy signal in the sense of (2.21). Note, in we took $A$ to be the identity, but these more general results are also proven there. This gives us a stopping criterion for our denoising algorithm.

In practice, we will use (2.21) as our stopping criterion.

Convergence

We begin with the following simple results for the linearized Bregman iteration or the equivalent algorithm (2.10).

If $u^{k}\rightarrow u^{\infty}$, then $Au^{\infty}=f$.

Assume $Au^{\infty}\not=f$. Then $A^{T}(Au^{\infty}-f)\not=0$ since $A^{T}$ has full rank. This means that for some $i$, $(A^{T}(Au^{k}-f))_{i}$ converges to a nonzero value, which means that $v_{i}^{k+1}-v_{i}^{k}$ does as well. On the other hand $\{v^{k}\}=\{u^{k}/\delta+p^{k}\}$ is bounded since $\{u^{k}\}$ converges and $p^{k}\in[-\mu,\mu]$. Therefore $\{v_{i}^{k}\}$ is bounded, while $v_{i}^{k+1}-v_{i}^{k}$ converges to a nonzero limit, which is impossible. ∎

If $u^{k}\rightarrow u^{\infty}$ and $v^{k}\rightarrow v^{\infty}$, then $u^{\infty}$ minimizes $\{J(u)+\frac{1}{2\delta}\|u\|^{2}:Au=f\}$.

Equation (2.16) (from a result in ) implies that $u^{\infty}$ will approach a solution to (1.1), (1.2), as $\mu$ approaches $\infty$.

The linearized Bregman iteration has the following monotonicity property:

If $u^{k+1}\not=u^{k}$ and $0<\delta<2/\|AA^{T}\|$, then

Then the shrinkage operation is such that

Let $Q^{k}=\text{Diag}~{}(q_{i}^{k})$. Then (3.23) can be written as

From (3.23), $Q^{k}$ is diagonal with $0\preceq Q^{k}\preceq I$, so $0\preceq AQ^{k}A^{T}\preceq AA^{T}$. If we choose $\delta>0$ such that $\delta AA^{T}\prec 2I$, then $0\preceq\delta AQ^{k}A^{T}\prec 2I$ or $-I\prec I-\delta AQ^{k}A^{T}\preceq I$ which implies that $\|Au^{k}-f\|$ is not increasing. To get strict decay, we need only show that $AQ^{k}A^{T}(Au^{k}-f)=0$ is impossible if $u^{k+1}\not=u^{k}$. Suppose $AQ^{k}A^{T}(Au^{k}-f)=0$ holds, then from (3.24) we have:

By (3.23), this only happens if $u_{i}^{k+1}=u_{i}^{k}$ for all $i$, which is a contradiction. ∎

We are still faced with estimating how fast the residual decays. It turns out that if consecutive elements of $u$ do not change sign, then $\|Au-f\|$ decays exponentially. By ’exponential’ we mean that the ratio of the residuals of two consecutive iteration converges to a constant, this type of convergence is sometimes called linear convergence. Here we define

(where $\text{sign}(0)=0$ and $\text{sign}(a)=a/|a|$ for $a\neq 0$). Then we have the following:

If $u^{k}\in S\equiv S_{u_{k}}$ for $k\in(T_{1},T_{2})$, then $u^{k}$ converges to $u^{*}$, where $u^{*}\in\arg\min\{\|Au-f\|^{2}:u\in S\}$ and $\|Au^{k}-f\|^{2}$ decays to $\|Au^{*}-f\|^{2}$ exponentially.

. Since $u^{k}\in S$ for $k\in[T_{1},T_{2}]$, we can define $Q\equiv Q^{k}$ for $T_{1}\leq k\leq T_{2}-1$. From (3.23) we see that $Q^{k}$ is a diagonal matrix consisting of zeros or ones, so $Q=Q^{T}Q$. Moreover, it is easy to see that $S=\{x|Qx=x\}$.

Following the argument in Theorem 3.3 we have:

Let $R^{n}=V_{0}\oplus V_{1}$, where $V_{0}$ is the null space of $AQA^{T}$ and $V_{1}$ is spanned by the eigenvectors corresponding to the nonzero eigenvalues of $AQA^{T}$. Let $Au^{k}-f=w^{k,0}+w^{k,1}$, where $w^{k,j}\in V_{j}$ for $j=0,1$. From (3.28) we have

for $T_{1}\leq k\leq T_{2}-1$. Since $w^{k,1}$ is not in the null space of $AQA^{T}$, then (3.27) and (3.28) imply that $\|w^{k,1}\|$ decays exponentially. Let $w^{0}=w^{k,0}$, then $AQA^{T}w^{0}=0$ $AQQA^{T}w^{0}\Rightarrow QA^{T}w^{0}=0$. Therefore, from (3.27) we have

Thus $\|\Delta u^{k}\|$ decays exponentially. This means $\{u^{k}\}$ forms a Cauchy sequence in $S$, so it has a limit $u^{*}\in S$. Moreover

Since $V_{0}$ and $V_{1}$ are orthogonal:

so $\|Au^{k}-f\|^{2}-\|Au^{*}-f\|^{2}$ decays exponentially. The only thing left to show is that

This is equivalent to way that $A^{T}(Au^{*}-f)$ is orthogonal with the hyperspace $\{u:Qu=u\}$. It’s easy to see that since $Q$ is a projection operator, a vector $v$ is orthogonal with $\{u:Qu=u\}$ if and only if $Qv=0$, thus we need to show $QA^{T}(Au^{*}-f)=0$. This is obvious because we have shown that $Au^{*}-f=w^{0}$ and $QA^{T}w^{0}=0$. So we find that $u^{*}$ is the desired minimizer. ∎

Therefore, instead of decaying exponentially with a global rate, the residual of the linearized Bregman iteration decays in a rather sophisticated manner. From the definition of the shrinkage function we can see that the sign of an element of $u$ will change if and only if the corresponding element of $v$ crosses the boundary of the interval $[-\mu,\mu]$. If $\mu$ is relatively large compared with the size of $\Delta v$ (which is usually the case when applying the algorithm to a compressed sensing problem), then at most iterations the signs of the elements of $u$ will stay unchanged, i.e. $u$ will stay in the subspace $S_{u}$ defined in (3.26) for a long while. This theorem tells us that under this scenario $u$ will quickly converge to the point $u^{*}$ that minimizes $\|Au-f\|$ inside $S_{u}$, and the difference between $\|Au-f\|$ and $\|Au^{*}-f\|$ decays exponentially. After $u$ converges to $u^{*}$, $u$ will stay there until the sign of some element of $u$ changes. Usually this means that a new nonzero element of $u$ comes up. After that, $u$ will enter a different subspace $S$ and a new converging procedure begins.

The phenomenon described above can be observed clearly in Fig 1. The final solution of $u$ contains five non-zero spikes. Each time a new spike appears, it converges rapidly to the position that minimizes $\|Au-f\|$ in the subspace $S_{u}$. After that there is a long stagnation, which means $u$ is just waiting there until the accumulating $v$ brings out a new non-zero element of $u$. The larger $\mu$ is, the longer the stagnation takes. Although the convergence of the residual during each phase is fast, the total speed of the convergence suffers much from the stagnation. The solution of this problem will be described in the next section.

Fast Implementation

The iterative formula in Algorithm 1 below gives us the basic linearized Bregman algorithm designed to solve (1.1),(1.2).

This is an extremely concise algorithm, simple to program, involve only matrix multiplication and shrinkage. When $A$ consists of rows of a matrix of a fast transform like FFT which is a common case for compressed sensing, it is even faster because matrix multiplication can be implemented efficiently using the existing fast code of the transform. Also, storage becomes a less serious issue.

We now consider how we can accelerate the algorithm under the problem of stagnation described in the previous section. From that discussion, during a stagnation $u$ converges to a limit $u^{*}$ so we will have $u^{k+1}\approx u^{k+2}\approx\cdots\approx u^{k+m}\approx u^{*}$ for some $m$. Therefore the increment of $v$ in each step, $A^{\top}(f-Au)$, is fixed. This implies that during the stagnation $u$ and $v$ can be calculated explicitly as following

If we denote the set of indices of the zero elements of $u^{*}$ as $I_{0}$ and let $I_{1}=\overline{I_{0}}$ be the support of $u^{*}$, then $v^{k}_{i}$ will keep changing only for $i\in I_{0}$ and the iteration can be formulated entry-wise as:

for $j=1,\cdots,m$. The stagnation will end when $u$ begins to change again. This happens if and only if some element of $v$ in $I_{0}$ (which keeps changing during the stagnation) crosses the boundary of the interval $[-\mu,\mu]$. When $i\in I_{0}$, $v^{k}_{i}\in[-\mu,\mu]$, so we can estimate the number of the steps needed for $v_{i}^{k}$ to cross the boundary $\forall i\in I_{0}$ from (4.30), which is

is the number of steps needed. Therefore, $s$ is nothing but the length of the stagnation. Using (4.29), we can predict the end status of the stagnation by

Therefore, we can kick $u$ to the critical point of the stagnation when we detect that $u$ has been staying unchanged for a while. Specifically, we have the following algorithm: Algorithm 2.

Indeed, this kicking procedure is similar to line search commonly used in optimization problems and modifies the initial algorithm in no way but just accelerates the speed. More precisely, note that the output sequence $\{u^{k},v^{k}\}$ is a subsequence of the original one, so all the previous theoretical conclusions on convergence still hold here.

An example of the algorithm is shown in Fig 2. It is clear that all the stagnation in the original convergence collapses to single steps. The total amount of computation is reduced dramatically.

Numerical Results

In this section, we demonstrate the effectiveness of the algorithm (with kicking) in solving basis pursuit and some related problems.

Consider the constrained minimization problem

where the constraints $Au=f$ are under-determined linear equations with $A$ an $m\times n$ matrix, and $f$ generated from a sparse signal $\bar{u}$ that has a number of nonzeros $\kappa<m$.

Our numerical experiments use two types of $A$ matrices: Gaussian matrices whose elements were generated from i.i.d. normal distributions $\mathcal{N}(0,1)$ (randn(m,n) in MATLAB), and partial discrete cosine transform (DCT) matrices whose $k$ rows were chosen randomly from the $n\times n$ DCT matrix. These matrices are known to be efficient for compressed sensing. The number of rows $m$ is chosen as $m\sim\kappa\log(n/\kappa)$ for Gaussian matrices and $m\sim\kappa\log n$ for DCT matrices (following ).

Note that partial DCT matrices are implicitly stored fast transforms for which matrix-vector multiplications in the forms of $Ax$ and $A^{\top}x$ were computed by the MATLAB commands dct(x) and idct(x), respectively. Therefore, we were able to test on partial DCT matrices of much larger sizes than Gaussian matrices. The sizes $m$-by-$n$ of these matrices are given in the first two columns of Table 1.

Our code was written in MATLAB and was run on a Windows PC with a Intel(R) Core(TM) 2 Duo 2.0GHz CPU and 2GB memory. The MATLAB version is 7.4.

The set of computational results given in Table 1 was obtained by using the stopping criterion

which was sufficient to give a small error $\|u^{k}-\bar{u}\|/\|\bar{u}\|$. Throughout our experiments in Table 1, we used $\mu=1$ to ensure the correctness of the results.

2 Robustness to Noise

In real applications, the measurement $f$ we obtain is usually contaminated by noise. The measurement we have is:

To characterize the noise level, we shall use SNR (signal to noise ratio) instead of $\sigma$ itself. The SNR is defined as follows

3 Recovery of Signal with High Dynamical Range

In this section, we test our algorithm on signals with high dynamical ranges. Precisely speaking, let $\mbox{MAX}=\max\{|\bar{u}_{i}|:i=1,\ldots,n\}$ and $\mbox{MIN}=\min\{|u_{i}|:u_{i}\neq 0,i=1,\ldots,n\}$. The signals we shall consider here satisfy $\frac{\mbox{MAX}}{\mbox{MIN}}\approx 10^{10}$. Our $\bar{u}$ is generated by multiplying a random number in $$with another one randomly picked from$\{1,10,\ldots,10^{10}\}$. Here we adopt the stopping criteria

for the case without noise (Figure 4) and the same stopping criteria as in the previous section for the noisy cases (Figures 5-7). In the experiments, we take the dimension $n=4000$, the number of nonzeros of $\bar{u}$ to be $0.02n$, and $\mu=10^{10}$. Here $\mu$ is chosen to be much larger than before, because the dynamical range of $\bar{u}$ is large. Figure 4 shows results for the noise free case, where the algorithm converges to a $10^{-11}$ residual in less than 300 iterations. Figures 5-7 show the cases with noise (the noise is added the same way as in previous section). As one can see, if the measurements are contaminated with less noise, signals with smaller magnitudes will be recovered well. For example in Figure 5, the SNR$\approx 118$, and the entries of magnitudes $10^{4}$ are well recovered; in Figure 6, the SNR$\approx 97$, and the entries of magnitudes $10^{5}$ are well recovered; and in Figure 7, the SNR$\approx 49$, and the entries of magnitudes $10^{7}$ are well recovered.

4 Recovery of Sinusoidal Waves in Huge Noise

Conclusion

We have proposed the linearized Bregman iterative algorithms as a competitive method for solving the compressed sensing problem. Besides the simplicity of the algorithm, the special structure of the iteration enables the kicking scheme to accelerate the algorithm even when $\mu$ is extremely large. As a result, a sparse solution can always be approached efficiently.

It also turns out that our process has remarkable denoising properties for undersampled sparse signals. We will pursue this in further work.

Our results suggest there is a big category of problem that can be solved by linearized Bregman iterative algorithms. We hope that our method and its extensions could produce even more applications for problems under different scenarios, including very underdetermined inverse problems in partial differential equations.

Acknowledgements

S.O. was supported by ONR Grant N000140710810, a grant from the Department of Defense and NIH Grant UH54RR021813; Y.M. and B.D. were supported by NIH Grant UH54RR021813; W.Y. was supported by NSF Grant DMS-0748839 and an internal faculty research grant from the Dean of Engineering at Rice University.

References