Accelerated Projected Gradient Method for Linear Inverse Problems with Sparsity Constraints

Introduction

In many practical problems, one cannot observe directly the quantities of most interest; instead their values have to be inferred from their effect on observable quantities. When this relationship between observable $y$ and interesting quantity $f$ is (approximately) linear, as it is in surprisingly many cases, the situation can be modeled mathematically by the equation

where $A$ is a linear operator mapping a vector space $\mathcal{K}$ (which we assume to contain all possible “objects” $f$) to a vector space $\mathcal{H}$ (which contains all possible data $y$). The vector spaces $\mathcal{K}$ and $\mathcal{H}$ can be finite– or infinite–dimensional; in the latter case, we assume that $\mathcal{K}$ and $\mathcal{H}$ are (separable) Hilbert spaces, and that $A:\mathcal{K}\to\mathcal{H}$ is a bounded linear operator. Our main goal consists in reconstructing the (unknown) element $f\in\mathcal{K}$, when we are given $y$. If $A$ is a “nice”, easily invertible operator, and if the data $y$ are free of noise, then this is a trivial task. Often, however, the mapping $A$ is ill-conditioned or not invertible. Moreover, typically (1) is only an idealized version in which noise has been neglected; a more accurate model is

Several types of signals appearing in nature admit sparse frame expansions and thus, sparsity is a realistic assumption for a very large class of problems. For instance, natural images are well approximated by sparse expansions with respect to wavelets or curvelets .

Framework and Notations

Orthonormal bases are particular examples of frames, but there also exist many interesting frames in which the $\psi_{\lambda}$ are not linearly independent. Frames allow for a (stable) series expansion of any $f\in\mathcal{K}$ of the form

Several authors have proposed independently an iterative soft-thresholding algorithm to approximate the solution $\bar{x}(\tau)$ . More precisely, $\bar{x}(\tau)$ is the limit of sequences $x^{(n)}$ defined recursively by

We will call the iteration (7) the iterative soft-thresholding algorithm or the thresholded Landweber iteration.

Discussion of the Thresholded Landweber Iteration

We will call this the projected Landweber iteration.

The typical dynamics of this projected Landweber algorithm are illustrated in Fig. 1(a). The norm $\|x^{(n)}\|_{1}$ no longer overshoots $R$, but quickly takes on the limit value (i.e., $\|\bar{x}(\tau)\|_{1}$); the speed of convergence remains very slow, however. In this projected Landweber iteration case, modifying the iterations by introducing an adaptive “descent parameter” $\beta^{(n)}>0$ in each iteration, defining $x^{(n+1)}$ by

does lead, in numerical simulations, to promising, converging results (in which it differs from the soft-thresholded Landweber iteration, where introducing such a descent parameter did not lead to numerical convergence, as noted above).

The typical dynamics of this modified algorithm are illustrated in Fig. 1(b), which clearly shows the larger steps and faster convergence (when compared with the projected Landweber iteration in Fig. 1(a)). We shall refer to this modified algorithm as the projected gradient iteration or the projected steepest descent; it will be the main topic of this paper.

There exist results in the literature on convergence of projected gradient iterations, where the projections are (as they are here) onto convex sets, see, e.g., and references therein. These results treat iterative projected gradient methods in much greater generality than we need: they allow more general functionals than $\mathcal{D}$, and the convex set on which the iterative procedure projects need not be bounded. On the other hand, these general results typically have the following restrictions:

The convergence in infinite-dimensional Hilbert spaces (i.e., $\Lambda$ is countable but infinite) is proved only in the weak sense and often only for subsequences;

In the descent parameters are typically restricted to cases for which $\lim_{n\to\infty}\beta^{(n)}=0$. In , it is shown that the algorithm converges weakly for any choice of $\beta^{(n)}\in\left[\varepsilon,\frac{2-\varepsilon}{\|K\|}\right]$, for $\varepsilon>0$ arbitrarily small. Of most interest to us is the case where the $\beta^{(n)}$ are picked adaptively, can grow with $n$, and are not limited to values below $\frac{2}{\|K\|}$; this case is not covered by the methods of either or .

Before we get to this theorem, we need to build some more machinery first.

We first observe a useful property of the soft-thresholding operator.

A schematic illustration is given in Figure 2.

for all $x\neq b$. Since $\|b\|_{1}=R$, it follows that

The proof is standard for projection operators onto convex sets; we include it because its technique will be used often in this paper.

Switching the role of $x$ and $x^{\prime}$ one finds:

By combining these last two inequalities, one finds:

from which inequality (14) follows. $\Box$

The Projected Gradient Method

We begin with the following characterization of the minimizers of $\mathcal{D}$ on $B_{R}$.

for any $\beta>0$, which in turn is equivalent to the requirement that

It follows from this that, for all $w\in B_{R}$ and for all $\beta>0$,

The minimizer of $\mathcal{D}$ on $B_{R}$ need not be unique. We have, however

Moreover, if the vector $e$ is defined by

Remark 5.6 By changing, if necessary, signs of the canonical basis vectors, we can assume, without loss of generality, that $e_{\lambda}=+1$ for all $\lambda\in\Gamma$. We shall do so from now on. $\Box$

2 Weak convergence to minimizing accumulation points

We shall now impose some conditions on the $\beta^{(n)}$. We shall see examples in Section 6 where these conditions are verified.

Note that the choice $\beta^{(n)}=1$ for all $n$, which corresponds to the projected Landweber iteration, automatically satisfies Condition (B); since we shall show below that we obtain convergence when the $\beta^{(n)}$ satisfy Condition (B), this will then establish, as a corollary, convergence of the projected Landweber iteration algorithm (10) as well. We shall be interested in choosing, adaptively, larger values of $\beta^{(n)}$; in particular, we like to choose $\beta^{(n)}$ as large as possible.

Condition (B) is inspired by the standard length-step in the steepest descent algorithm for the (unconstrained, unpenalized) functional $\|Kx-y\|^{2}$. In this case, one can speed up the standard Landweber iteration $x^{(n+1)}=x^{(n)}+K^{*}(y-Kx^{(n)})$ by defining instead $x^{(n+1)}=x^{(n)}+\alpha K^{*}(y-Kx^{(n)})$, where $\alpha$ is picked so that it gives the largest decrease of $\|Kx-y\|^{2}$ in this direction. This gives

In this linear case, one easily checks that $\alpha$ also equals

in fact, it is this latter expression for $\alpha$ (which inspired the formulation of Condition (B)) that is most useful in proving convergence of the steepest descent algorithm.

Because the definition of $x^{(n+1)}$ involves $\beta^{(n)}$, the inequality (B2), which uses $x^{(n+1)}$ to impose a limitation on $\beta^{(n)}$, has an “implicit” quality. In practice, it may not be straightforward to pick $\beta^{(n)}$ appropriately; one could conceive of trying first a “greedy” choice, such as e.g. $\frac{\|r^{(n)}\|^{2}}{\|Kr^{(n)}\|^{2}}$; if this value works, it is retained; if it doesn’t, it can be gradually decreased (by multiplying it with a factor slightly smaller than 1) until (B2) is satisfied. (A similar way of testing appropriate step lengths is adopted in .)

Proof: First of all, observe that the functional $F_{\beta}(\cdot,x)$ is strictly convex, so that it has a unique minimizer on $B_{R}$; let $\hat{x}$ be this minimizer. Then for all $w\in B_{R}$ and for all $t\in$

Proof: Comparing the definition of $x^{(n+1)}$ in (11) with the statement of Lemma 5.9, we see that $x^{(n+1)}=T_{R}(\beta^{(n)};x^{(n)})$, so that $x^{(n+1)}$ is the minimizer, for $x\in B_{R}$, of $F_{\beta^{(n)}}(x;x^{(n)})$. Setting $\gamma=\frac{1}{r}-1>0$, we have

Therefore, the series $\sum_{n=0}^{\infty}\|x^{(n+1)}-x^{(n)}\|^{2}$ converges and $\lim_{n\to\infty}\|x^{(n+1)}-x^{(n)}\|=0$. $\Box$

In particular, specializing to our subsequence and taking the $\limsup$, we have

Because $\|x^{(n_{j})}-x^{(n_{j}+1)}\|\rightarrow 0$, for $j\to\infty$, and $w-x^{(n_{j}+1)}$ is uniformly bounded, we have

By adding $\beta^{(n_{j})}\langle K^{*}(y-Kx^{(n_{j}+1)}),x^{(n_{j}+1)}-x^{(n_{j})}\rangle$, which also tends to zero as $j\to\infty$, we transform this into

Since the $\beta^{(n_{j})}$ are all in $[1,\bar{\beta}]$, it follows that

where we have used the weak convergence of $x^{(n_{j})}$. This can be rewritten as

so that $x^{\#}$ is a minimizer of $\mathcal{D}$ on $B_{R}$, by Lemma 5.1. $\Box$

3 Strong convergence to minimizing accumulation points

In this subsection we show how the weak convergence established in the preceding subsection can be strengthened into norm convergence, again by a series of lemmas. Since the distinction between weak and strong convergence makes sense only when the index set $\Lambda$ is infinite, we shall implicitly assume this is the case throughout this section.

Proof: Specializing the inequality (39) to $w=x^{\#}$, we obtain

together with $\|Kx^{\#}\|^{2}\leq\liminf_{j\to\infty}\|Kx^{(n_{j})}\|^{2}$ (a consequence of the weak convergence of $Kx^{(n_{j})}$ to $Kx^{\#}$), this implies $\lim_{j\to\infty}\|K(x^{(n_{j})})\|^{2}=\|Kx^{\#}\|^{2}$, and thus $\lim_{j\to\infty}K(x^{(n_{j})})=Kx^{\#}$. $\Box$

Combining this with $\|u^{(j)}-v^{(j)}\|\,_{\overrightarrow{j\to\infty}}\,0$, we obtain

The following proposition summarizes in one statement all the findings of the last two subsections.

4 Uniqueness of the accumulation point

Because $M_{R}$ is convex, this projection operator has the following property:

Because $a^{(n)}$ is a minimizer, we can also apply Lemma 5.1 to $a^{(n)}$ and conclude

With these inequalities, we can prove the following crucial result.

where we have used $Ka^{(n)}=Ka^{(n+1)}$. It follows that

Adding $\frac{1}{2}\beta^{(n)}\|K(b^{(n)}-b^{(n+1)})\|^{2}\leq\frac{r}{2}\|x^{(n)}-x^{(n+1)}\|^{2}$ to (60), we have

We are now ready to state the main result of our work.

Numerical Experiments and Additional Algorithms

We conduct a number of numerical experiments to gauge the effectiveness of the different algorithms we discussed. All computations were done in Mathematica 5.2 on a 2Ghz workstation with 2Gb memory.

We are primarily interested in the behavior, as a function of time (not number of iterations), of the relative error $\|x^{(n)}-\bar{x}\|/\|\bar{x}\|$. To this end, and for a given operator $K$ and data $y$, we need to know in advance the actual minimizer $\bar{x}(\tau)$ of the functional (5).

One can calculate the minimizer exactly (in practice up to computer round-off) with a finite number of steps using the LARS algorithm described in (the variant called ‘Lasso’, implemented independently by us). This algorithm scales badly, and is useful in practice only when the number of non-zero entries in the minimizer $\bar{x}(\tau)$ is sufficiently small. We made our own implementation of this algorithm to make it more directly applicable to our problem (i.e., we do not renormalize the columns of the matrix to have zero mean and unit variance, as it is done in the statistics context ). We also double-check the minimizer obtained in this manner by verifying that it is indeed a fixed point of the iterative thresholding algorithm (7) (up to machine epsilon). We then have an ‘exact’ minimizer $\bar{x}$ together with its radius $R=\|\bar{x}\|_{1}$ (used in the projected algorithms) and, according to Lemma 5.3, the corresponding threshold $\tau=\max_{i}|\bar{r}_{i}|$ with $\bar{r}=K^{\ast}(y-K\bar{x})$ (used in the iterative thresholding algorithm).

The numerical examples below are listed in order of increasing complexity; they illustrate that the algorithms can behave differently for different examples. In these experiments we choose $\beta^{(n)}=\beta^{(n)}_{\mbox{\tiny{st.}}}:=\|r^{(n)}\|^{2}/\|Kr^{(n)}\|^{2}$, (where, as before, $r^{(n)}=K^{\ast}(y-Kx^{(n)})$); $\beta^{(n)}_{\mbox{\tiny{st.}}}$ is the standard descent parameter from the classical linear steepest descent algorithm.

When $K$ is a partial Fourier matrix (i.e., a Fourier matrix with a prescribed number of deleted rows), there is no advantage in using a dynamical step size $\beta_{\mbox{\tiny{st.}}}^{(n)}=\|r^{(n)}\|^{2}/\|Kr^{(n)}\|^{2}$ as this ratio is always equal to 1. This trivially fulfills Condition (B) in Section 5.1. The performance of the projected steepest descent iteration simply equals that of the projected Landweber iterations.

By combining a scaled partial Fourier transform with a rank 1 projection operator, we constructed our second example, in which $K$ is a $1536\times 2049$ matrix, of rank $1536$, with largest singular value equal to 0.99 and all the other singular values between 0.01 and 0.11. Because of the construction of the matrix, the FFT algorithm provides a fast way of computing the action of this matrix on a vector. For the $y$ and $\tau$ that were chosen, the limit vector $\bar{x}_{\tau}$ has $429$ nonzero entries. For this example, the LARS procedure is slower than thresholded Landweber, which in turn is significantly slower than projected steepest descent. To get within a distance of the true minimizer corresponding to a $5\%$ relative error, the projected steepest descent algorithm takes $2\sec$, the thresholded Landweber algorithm $39\sec$, and LARS $151\sec$. (The relatively poor performance of LARS in this case is due to the large number of nonzero entries in the limit vector $\bar{x}_{\tau}$; the complexity of LARS is cubic in this number of nonzero entries.) In this case, the $\beta^{(n)}_{\mbox{\tiny{st.}}}=\|r^{(n)}\|^{2}/\|Kr^{(n)}\|^{2}$ are much larger than 1; moreover, they satisfy Condition (B) of Section 5.1 at every step. We illustrate the results in Figure 3.

The last example is inspired by a real-life application in geoscience , in particular an application in seismic tomography based on earthquake data. The object space consists of the wavelet coefficients of a 2D seismic velocity perturbation. There are $8192$ degrees of freedom. In this particular case the number of data is $1848$. Hence the matrix $K$ has $1848$ rows and $8192$ columns. We apply the different methods to the same noisy data that are used in and measure the time to convergence up to a specified relative error (see Table 1 and Figure 4). This example illustrates the slow convergence of the thresholded Landweber algorithm (7), and the improvements made by a projected steepest descent iteration (11) with the special choice $\beta^{(n)}=\beta^{(n)}_{\mbox{\tiny{st.}}}$ above. In this case, this choice turns out not to satisfy Condition (B) in general. One could conceivably use successive corrections, e.g. by a line-search, to determine, starting from $\beta^{(n)}_{\mbox{\tiny{st.}}}$, values of $\beta^{(n)}$ that would satisfy condition (B), and thus guarantee convergence as established by Theorem 5.18. This would slow down the method considerably. The $\beta^{(n)}_{\mbox{\tiny{st.}}}$ seem to be in the right ballpark, and provide us with a numerically converging sequence. We also implemented the projected Landweber algorithm (10); it is listed in Table 1 and illustrated in Figure 4.

The matrix $K$ in this example is extremely ill-conditioned: its largest singular value was normalized to 1, but the remaining singular values quickly tend to zero. The threshold was chosen, according to the (known or estimated) noise level in the data, so that $\mathcal{D}(\bar{x})/\sigma^{2}=1848$ ( = the number of data points), where $\sigma$ is the data noise level; this is a standard choice that avoids overfitting.

It is worthwhile noticing that for the three algorithms the value of the functional (5) converges much faster to its limit value than the minimizer itself: When the reconstruction error is 10%, the corresponding value of the functional is already accurate up to three digits with respect to the value of the functional at $\bar{x}$. We can imagine that in this case the functional has a long narrow “valley” with a very gentle slope in the direction of the eigenvectors with small (or zero) singular values.

In this particular example, the number of nonzero components of $\bar{x}$ equals $128$. The LARS (exact) algorithm only takes $55$ seconds, which is much faster than any of the iterative methods demonstrated here. However, as illustrated above, by the second example, LARS looses its advantage when dealing with larger problems where the minimizer is not sparse in absolute number of entries, as is the case in, e.g., realistic problems of global seismic tomography. Indeed, the example presented here is a “toy model” for proof-of-concept for geoscience applications. The 3D model will involve millions of unknowns and solutions that may be sparse compared with the total number of unknowns, but not sparse in absolute numbers. Because the complexity of LARS is cubic in the number of nonzero components of the solution, such 3D problems are expected to lie beyond its useful range.

2 Relationship to other methods

The projected iterations (16) and (17) are related to the POCS (Projection on Convex Sets) technique . The projection of a vector $a$ on the solution space $\{x:Kx=y\}$ (a convex set, assumed here to be non-empty; no such assumption was made before because the functional (5) always has a minimum) is given by:

This may be practical in case of a small number of data or when there is structure in $K$, i.e., when $KK^{*}$ is efficiently inverted. Approximating $KK^{\ast}$ by the unit matrix, yields the projected Landweber algorithm (16); approximating $(KK^{*})^{-1}$ by a constant multiple of the unit matrix yields the projected gradient iteration (17) if one chooses the constant equal to $\beta^{(n)}$.

Conclusions

There is no universal method that performs best for any choice of the operator, data, and penalization parameter. As a general rule of thumb we expect that, among the algorithms discussed in this paper for which we have convergence proofs,

the thresholded Landweber algorithm (7) works best for an operator $K$ close to the identity (independently of the sparsity of the limit),

the projected steepest descent algorithm (11) works best for an operator with a relatively nice spectrum, i.e., with not too many zeroes (also independently of the sparsity of the minimizer), and

the exact (LARS) method works best when the minimizer is sparse in absolute terms.

Obviously, the three cases overlap partially, and they do not cover the whole range of possible operators and data. In future work we intend to investigate algorithms that would further improve the performance for the case of a large ill-conditioned matrix and a minimizer that is relatively sparse with respect to the dimension of the underlying space. We intend, in particular, to focus on proving convergence and other mathematical properties of (67).

Acknowledgments

M. F. acknowledges the financial support provided by the European Union’s Human Potential Programme under the contract MOIF-CT-2006-039438. I. L. is a post-doctoral fellow with the F.W.O.-Vlaanderen (Belgium). M.F. and I.L. thank the Program in Applied and Computational Mathematics, Princeton University, for the hospitality during the preparation of this work. I. D. gratefully acknowledges partial support from NSF grants DMS-0245566 and 0530865.

References