Fast Image Recovery Using Variable Splitting and Constrained Optimization

I Introduction

Image restoration/reconstruction is one of the earliest and most classical linear inverse problems in imaging, dating back to the 1960’s . In this class of problems, a noisy indirect observation ${\bf y}$, of an original image ${\bf x}$, is modeled as

In the particular case of image deblurring/deconvolution, ${\bf B}$ is the matrix representation of a convolution operator; if this convolution is periodic, ${\bf B}$ is then a (block) circulant matrix. This type of observation model describes well several physical mechanisms, such as relative motion between the camera and the subject (motion blur), bad focusing (defocusing blur), or a number of other mechanisms which are well modeled by a convolution.

In more general image reconstruction problems, ${\bf B}$ represents some linear direct operator, such as a set of tomographic projections (Radon transform), a partially observed (e.g., Fourier) transform, or the loss of part of the image pixels.

An alternative formulation applies a regularizer directly to the unknown image, leading to criteria of the form

Finally, it should be mentioned that problems (1) and (2) can be seen as the Lagrangians of associated constrained optimization problems: (1) is the Lagrangian of the constrained problem

Specifically, a solution of (3) (for any $\varepsilon$ such that this problem is feasible) is either the null vector, or else is a minimizer of (1), for some $\tau>0$ (see [39, Theorem 27.4]). A similar relationship exists between problems (2) and (4).

I-B Previous Algorithms

For any problem of non-trivial dimension, matrices ${\bf BW}$, ${\bf B}$, and ${\bf W}$ cannot be stored explicitly, and it is costly, even impractical, to access portions (lines, columns, blocks) of them. On the other hand, matrix-vector products involving ${\bf B}$ or ${\bf W}$ (or their conjugate transposes ${\bf B}^{H}$ and ${\bf W}^{H}$) can be done quite efficiently. For example, if the columns of ${\bf W}$ contain a wavelet basis or a tight wavelet frame, any multiplication of the form ${\bf W}{\bf v}$ or ${\bf W}^{H}{\bf v}$ can be performed by a fast wavelet transform algorithm . Similarly, if ${\bf B}$ represents a convolution, products of the form ${\bf B}{\bf v}$ or ${\bf B}^{H}{\bf v}$ can be performed with the help of the fast Fourier transform (FFT) algorithm. These facts have stimulated the development of special purpose methods, in which the only operations involving ${\bf B}$ or ${\bf W}$ (or their conjugate transposes) are matrix-vector products.

To present a unified view of algorithms for handling (1) and (2), we write them in a common form

where ${\bf A}={\bf BW}$, in the case of (1), while ${\bf A}={\bf B}$, for (2).

Notice that if $\phi$ is proper and convex, the function being minimized is proper and strictly convex, thus the minimizer exists and is unique making the function well defined .

Each IST iteration for solving (5) is given by

where $1/\gamma$ is a step size. Notice that ${\bf A}^{H}\left({\bf A}{\bf x}_{k}-{\bf y}\right)$ is the gradient of the data-fidelity term $(1/2)\|{\bf A}{\bf x}-{\bf y}\|_{2}^{2}$, computed at ${\bf x}_{k}$; thus, each IST iteration takes a step of length $1/\gamma$ in the direction of the negative gradient of the data-fidelity term, followed by the application of the shrinkage/thresholding function associated with the regularizer $\phi$.

It has been shown that if $\gamma>\|{\bf A}\|_{2}^{2}/2$ and $\phi$ is convex, the algorithm converges to a solution of (1) . However, it is known that IST may be quite slow, specially when $\tau$ is very small and/or the matrix ${\bf A}$ is very ill-conditioned , , , . This observation has stimulated work on faster variants of IST, which we will briefly review in the next paragraphs.

In the two-step IST (TwIST) algorithm , each iterate depends on the two previous iterates, rather than only on the previous one (as in IST). This algorithm may be seen as a non-linear version of the so-called two-step methods for linear problems . TwIST was shown to be considerably faster than IST on a variety of wavelet-based and TV-based image restoration problems; the speed gains can reach up two orders of magnitude in typical benchmark problems.

Another two-step variant of IST, named fast IST algorithm (FISTA), was recently proposed and also shown to clearly outperform IST in terms of speed . FISTA is a non-smooth variant of Nesterov’s optimal gradient-based algorithm for smooth convex problems , .

A strategy recently proposed to obtain faster variants of IST consists in relaxing the condition $\gamma>\gamma_{\mbox{\tiny min}}\equiv\|{\bf A}\|_{2}^{2}/2$. In the SpaRSA (standing for sparse reconstruction by separable approximation) framework , , a different $\gamma_{t}$ is used in each iteration (which may be smaller than $\gamma_{\mbox{\tiny min}}$, meaning larger step sizes). It was shown experimentally that SpaRSA clearly outperforms standard IST. A convergence result for SpaRSA was also given in .

Finally, when the slowness is caused by the use of a small value of the regularization parameter, continuation schemes have been found quite effective in speeding up the algorithm. The key observation is that IST algorithm benefits significantly from warm-starting, i.e., from being initialized near a minimum of the objective function. This suggests that we can use the solution of (5), for a given value of $\tau$, to initialize IST in solving the same problem for a nearby value of $\tau$. This warm-starting property underlies continuation schemes , , . The idea is to use IST to solve (1) for a larger value of $\tau$ (which is usually fast), then decrease $\tau$ in steps toward its desired value, running IST with warm-start for each successive value of $\tau$.

I-C Proposed Approach

The approach proposed in this paper is based on the principle of variable splitting, which goes back at least to Courant in the 40’s , . Since the objective function (5) to be minimized is the sum of two functions, the idea is to split the variable ${\bf x}$ into a pair of variables, say ${\bf x}$ and ${\bf v}$, each to serve as the argument of each of the two functions, and then minimize the sum of the two functions under the constraint that the two variables have to be equal, so that the problems are equivalent. Although variable splitting is also the rationale behind the recently proposed split-Bregman method , in this paper, we exploit a different type of splitting to attack problem (5). Below we will explain this difference in detail.

The obtained constrained optimization problem is then dealt with using an augmented Lagrangian (AL) scheme , which is known to be equivalent to the Bregman iterative methods recently proposed to handle imaging inverse problems (see and references therein). We prefer the AL perspective, rather than the Bregman iterative view, as it is a standard and more elementary optimization tool (covered in most textbooks on optimization). In particular, we solve the constrained problem resulting from the variable splitting using an algorithm known as alternating direction method of multipliers (ADMM) .

The application of ADMM to our particular problem involves solving a linear system with the size of the unknown image (in the case of problem (2)) or with the size of its representation (in the case of problem (1)). Although this seems like an unsurmountable obstacle, we show that it is not the case. In many problems of the form (2), such as deconvolution, recovery of missing samples, or reconstruction from partial Fourier observations, this system can be solved very quickly in closed form (with $O(n)$ or $O(n\log n)$ cost). For problems of the form (1), we show how exploiting the fact that ${\bf W}$ is a tight Parseval frame, this system can still be solved efficiently (typically with $O(n\log n)$ cost.

We report results of a comprehensive set of experiments, on a set of benchmark problems, including image deconvolution, recovery of missing pixels, and reconstruction from partial Fourier transform, using both frame-based and TV-based regularization. In all the experiments, the resulting algorithm is consistently and considerably faster than the previous state of the art methods FISTA , TwIST , and SpaRSA .

The speed of the proposed algorithm, which we term SALSA (split augmented Lagrangian shrinkage algorithm), comes from the fact that it uses (a regularized version of) the Hessian of the data fidelity term of (5), that is, ${\bf A}^{H}{\bf A}$, while the above mentioned algorithms essentially only use gradient information.

I-D Organization of the Paper

Section II describes the basic ingredients of SALSA: variable splitting, augmented Lagrangians, and ADMM. In Section III, we show how these ingredients are combined to obtain the proposed SALSA. Section IV reports experimental results, and Section V ends the paper with a few remarks and pointers to future work.

II Basic Ingredients

Consider an unconstrained optimization problem in which the objective function is the sum of two functions, one of which is written as the composition of two functions,

which is clearly equivalent to unconstrained problem (8): in the feasible set $\{({\bf u},{\bf v}):g({\bf u})={\bf v}\}$, the objective function in (9) coincides with that in (8). The rationale behind variable splitting methods is that it may be easier to solve the constrained problem (9) than it is to solve its unconstrained counterpart (8).

The splitting idea has been recently used in several image processing applications. A variable splitting method was used in to obtain a fast algorithm for TV-based image restoration. Variable splitting was also used in to handle problems involving compound regularizers; i.e., where instead of a single regularizer $\tau\phi({\bf x})$ in (5), one has a linear combination of two (or more) regularizers $\tau_{1}\phi_{1}({\bf x})+\tau_{2}\phi_{2}({\bf x})$. In and , the constrained problem (9) is attacked by a quadratic penalty approach, i.e., by solving

by alternating minimization with respect to ${\bf u}$ and ${\bf v}$, while slowly taking $\alpha$ to very large values (a continuation process), to force the solution of (10) to approach that of (9), which in turn is equivalent to (8). The rationale behind these methods is that each step of this alternating minimization may be much easier than the original unconstrained problem (8). The drawback is that as $\alpha$ becomes very large, the intermediate minimization problems become increasingly ill-conditioned, thus causing numerical problems (see , Chapter 17).

A similar variable splitting approach underlies the recently proposed split-Bregman methods ; however, instead of using a quadratic penalty technique, those methods attack the constrained problem directly using a Bregman iterative algorithm . It has been shown that, when $g$ is a linear function, i.e., $g({\bf u})={\bf Gu}$, the Bregman iterative algorithm is equivalent to the augmented Lagrangian method , which is briefly reviewed in the following subsection.

II-B Augmented Lagrangian

Consider the constrained optimization problem

The so-called augmented Lagrangian method (ALM) , also known as the method of multipliers (MM) , , consists in minimizing {\cal L}_{A}({\bf z},\mbox{\boldmath\lambda},\mu) with respect to ${\bf z}$, keeping $\lambda$ fixed, then updating $\lambda$, and repeating these two steps until some convergence criterion is satisfied. Formally, the ALM/MM works as follows:

Set $k=0$, choose $\mu>0$, ${\bf z}_{0}$, and \mbox{\boldmath\lambda}_{0}.

{\bf z}_{k+1}\in\arg\min_{{\bf z}}{\cal L}_{A}({\bf z},\mbox{\boldmath\lambda}_{k},\mu)

\mbox{\boldmath\lambda}_{k+1}=\mbox{\boldmath\lambda}_{k}+\mu({\bf b}-{\bf Az}_{k+1})

It is also possible (and even recommended) to update the value of $\mu$ in each iteration , [3, Chap. 9]. However, unlike in the quadratic penalty approach, the ALM/MM does not require $\mu$ to be taken to infinity to guarantee convergence to the solution of the constrained problem (11).

Notice that (after a straightforward complete-the-squares procedure) the terms added to $E({\bf z})$ in the definition of the augmented Lagrangian {\cal L}_{A}({\bf z},\mbox{\boldmath\lambda}_{k},\mu) in (12) can be written as a single quadratic term (plus a constant independent of ${\bf z}$, thus irrelevant for the ALM/MM), leading to the following alternative form of the algorithm (which makes clear its equivalence with the Bregman iterative method ):

Set $k=0$, choose $\mu>0$ and ${\bf d}_{0}$.

${\bf z}_{k+1}\in\arg\min_{{\bf z}}E({\bf z})+\frac{\mu}{2}\|{\bf Az-d}_{k}\|_{2}^{2}$

${\bf d}_{k+1}={\bf d}_{k}+({\bf b}-{\bf Az}_{k+1})$

It has been shown that, with adequate initializations, the ALM/MM generates the same sequence as a proximal point algorithm applied to the Lagrange dual of problem (11) . Moreover, the sequence $\{{\bf d}_{k}\}$ converges to a solution of this dual problem and all cluster points of the sequence $\{{\bf z}_{k}\}$ are solutions of the (primal) problem (11) .

II-C ALM/MM for Variable Splitting

We now show how the ALM/MM can be used to address problem (9), in the particular case where $g({\bf u})={\bf Gu}$, i.e.,

With these definitions in place, Steps 3 and 4 of the ALM/MM (version II) can be written as follows:

The minimization problem (16) is not trivial since, in general, it involves non-separable quadratic and possibly non-smooth terms. A natural to address (16) is to use a non-linear block-Gauss-Seidel (NLBGS) technique, in which (16) is solved by alternatingly minimizing it with respect to ${\bf u}$ and ${\bf v}$, while keeping the other variable fixed. Of course this raises several questions: for a given ${\bf d}_{k}$, how much computational effort should be spent in approximating the solution of (16)? Does this NLBGS procedure converge? Experimental evidence in suggests that an efficient algorithm is obtained by running just one NLBGS step. It turns out that the resulting algorithm is the so-called alternating direction method of multipliers (ADMM) , which works as follows:

Set $k=0$, choose $\mu>0$, ${\bf v}_{0}$, and ${\bf d}_{0}$.

${\bf u}_{k+1}\in\arg\min_{{\bf u}}f_{1}({\bf u})+\frac{\mu}{2}\|{\bf G}{\bf u}-{\bf v}_{k}-{\bf d}_{k}\|_{2}^{2}$

${\bf v}_{k+1}\in\arg\min_{{\bf v}}f_{2}({\bf v})+\frac{\mu}{2}\|{\bf G}{\bf u}_{k+1}-{\bf v}-{\bf d}_{k}\|_{2}^{2}$

${\bf d}_{k+1}={\bf d}_{k}+{\bf Gu}_{k+1}-{\bf v}_{k+1}$

For later reference, we now recall the theorem by Eckstein and Bertsekas, in which convergence of (a generalized version of) ADMM is shown. This theorem applies to problems of the form (8) with $g({\bf u})={\bf G}{\bf u}$, i.e.,

of which (13) is the constrained optimization reformulation.

Then, if (18) has a solution, the sequence $\{{\bf u}_{k}\}$ converges, ${\bf u}_{k}\rightarrow{\bf u}^{*}$, where ${\bf u}^{*}$ is a solution of (18). If (18) does not have a solution, then at least one of the sequences $\{{\bf v}_{k}\}$ or $\{{\bf d}_{k}\}$ diverges.

Notice that the ADMM algorithm defined above generates sequences $\{{\bf u}_{k}\}$, $\{{\bf v}_{k}\}$, and $\{{\bf d}_{k}\}$ which satisfy the conditions in Theorem 1 in a strict sense (i.e., with $\eta_{k}=\mu_{k}=0$). One of the important consequences of this theorem is that it shows that it is not necessary to exactly solve the minimizations in lines 3 and 4 of ADMM; as long as sequence of errors is absolutely summable, convergence is not compromised.

The proof of Theorem 1 is based on the equivalence between ADMM and the so-called Douglas-Rachford splitting method (DRSM) applied to the dual of problem (18). The DRSM was recently used for image recovery problems in . For recent and comprehensive reviews of ALM/MM, ADMM, DRSM, and their relationship with Bregman and split-Bregman methods, see , .

III Proposed Method

We now return to the unconstrained optimization formulation of regularized image recovery, as defined in (5). This problem can be written in the form (18), with

The constrained optimization formulation is thus

At this point, we are in a position to clearly explain the difference between this formulation and the splitting exploited in split-Bregman methods (SBM) for image recovery . In those methods, the focus of attention is a non-separable regularizer that can be written as $\phi({\bf x})=\varphi({\bf D\,x})$, as is the case of the TV norm. The variable splitting used in SBM addresses this non-separability by defining the following constrained optimization formulation:

In contrast, we assume that the Moreau proximal mapping associated to the regularizer $\phi$, i.e., the function \mbox{\boldmath\Psi}_{\tau\phi}(\cdot) defined in (6), can be computed efficiently. The goal of our splitting is not to address the difficulty raised by a non-separable and non-quadratic regularizer, but to exploit second order (Hessian) information of the function $f_{1}$, as will be shown below.

III-B Algorithm and Its Convergence

Inserting the definitions given in (19)–(21) in the ADMM presented in the previous section yields the proposed SALSA (split augmented Lagrangian shrinkage algorithm).

Set $k=0$, choose $\mu>0$, ${\bf v}_{0}$, and ${\bf d}_{0}$.

${\bf x}^{\prime}_{k}={\bf v}_{k}+{\bf d}_{k}$

${\bf x}_{k+1}=\arg\min_{{\bf x}}\|{\bf A\,x}-{\bf y}\|_{2}^{2}+\mu\|{\bf x-x}^{\prime}_{k}\|_{2}^{2}$

${\bf v}^{\prime}_{k}={\bf x}_{k+1}-{\bf d}_{k}$

${\bf v}_{k+1}=\arg\min_{{\bf v}}\tau\phi({\bf v})+\frac{\mu}{2}\|{\bf v}-{\bf v}^{\prime}_{k}\|_{2}^{2}$

${\bf d}_{k+1}={\bf d}_{k}+{\bf x}_{k+1}-{\bf v}_{k+1}$

Notice that SALSA is an instance of ADMM with ${\bf G=I}$; thus, the full column rank condition on ${\bf G}$ in Theorem 1 is satisfied. If the minimizations in lines 4 and 6 are solved exactly, we can then invoke Theorem 1 to guarantee to convergence of SALSA.

In line 4 of SALSA, a strictly convex quadratic function has to be minimized; which leads to the following linear system

As shown in the next subsection, this linear system can be solved exactly (naturally, up to numerical precision), i.e., non-iteratively, for a comprehensive set of situations of interest. The matrix ${\bf A}^{H}{\bf A}+\mu\,{\bf I}$ can be seen as a regularized (by the addition of $\mu{\bf I}$) version of the Hessian of $f_{1}({\bf x})=(1/2)\|{\bf Ax-y}\|_{2}^{2}$, thus SALSA does use second order information of this function. Notice also that (24) is formally similar to the maximum a posteriori (MAP) estimate of ${\bf x}$, from observations ${\bf y}={\bf A}{\bf x}+{\bf n}$ (where ${\bf n}$ is white Gaussian noise of variance $1/\mu$) under a Gaussian prior of mean ${\bf x}^{\prime}_{k}$ and covariance ${\bf I}$.

The problem in line 6 is, by definition, the Moreau proximal mapping of $\phi$ applied to ${\bf v}^{\prime}_{k}$, thus its solution can be written as

If this mapping can be computed exactly in closed form, for example, if $\phi({\bf x})=\|{\bf x}\|_{1}$ thus $\Psi$ is simply a soft threshold, then, by Theorem 1, SALSA is guaranteed to converge. If $\Psi$ does not have a closed form solution and requires itself an iterative algorithm (e.g., if $\phi$ is the TV norm), then convergence of SALSA still holds if one can guarantee that the error sequence $\nu_{k}$ (see Theorem 1) is summable. This can be achieved (at least approximately) if the iterative algorithm used to approximate $\Psi$ is initialized with the result of the previous outer iteration, and a decreasing stopping threshold is used.

𝑘1{\bf x}_{k+1} As stated above, we are interested in problems where it is not feasible to explicitly form matrix ${\bf A}$; this might suggest that it is not easy, or even feasible, to compute the inverse in (24). However, as shown next, in a number of problems of interest, this inverse can be computed very efficiently.

In this case we have ${\bf A=B}$ (see (1), (2), and (5)), where ${\bf B}$ is the matrix representation of a convolution. This is the simplest case, since the inverse $\left({\bf B}^{H}{\bf B}+\mu\,{\bf I}\right)^{-1}$ can be computed in the Fourier domain. Although this is an elementary and well-known fact, we include the derivation for the sake of completeness. Assuming that the convolution is periodic (other boundary conditions can be addressed with minor changes), ${\bf B}$ is a block-circulant matrix with circulant blocks which can be factorized as

where ${\bf U}$ is the matrix that represents the 2D discrete Fourier transform (DFT), ${\bf U}^{H}={\bf U}^{-1}$ is its inverse (${\bf U}$ is unitary, i.e., ${\bf UU}^{H}={\bf U}^{H}{\bf U}={\bf I}$), and ${\bf D}$ is a diagonal matrix containing the DFT coefficients of the convolution operator represented by ${\bf B}$. Thus,

where $(\cdot)^{*}$ denotes complex conjugate and $|{\bf D}|^{2}$ the squared absolute values of the entries of the diagonal matrix ${\bf D}$. Since $|{\bf D}|^{2}+\mu\,{\bf I}$ is diagonal, its inversion has linear cost $O(n)$. The products by ${\bf U}$ and ${\bf U}^{H}$ can be carried out with $O(n\log n)$ cost using the FFT algorithm. The expression in (28) is a Wiener filter in the frequency domain.

III-C2 Deconvolution with Frame-Based Synthesis Prior

In this case, we have a problem of the form (1), i.e., ${\bf A=BW}$, thus the inversion that needs to be performed is $\left({\bf W}^{H}{\bf B}^{H}{\bf BW}+\mu\,{\bf I}\right)^{-1}$. Assuming that ${\bf B}$ represents a (periodic) convolution, this inversion may be sidestepped under the assumption that matrix ${\bf W}$ corresponds to a normalized tight frame (a Parseval frame), i.e., ${\bf W\,W}^{H}={\bf I}$. Applying the Sherman -Morrison- Woodbury (SMW) matrix inversion formula yields

Let’s focus on the term ${\bf F\equiv B}^{H}\left({\bf B}{\bf B}^{H}\!+\mu\,{\bf I}\right)^{-1}{\bf B}$; using the factorization (26), we have

Since all the matrices in ${\bf D^{*}}\left(|{\bf D}|^{2}+\mu\,{\bf I}\right)^{-1}{\bf D}$ are diagonal, this expression can be computed with $O(n)$ cost, while the products by ${\bf U}$ and ${\bf U}^{H}$ can be computed with $O(n\log n)$ cost using the FFT. Consequently, products by matrix ${\bf F}$ (defined in (29)) have $O(n\log n)$ cost.

Defining ${\bf r}_{k}=\left({\bf A}^{H}{\bf y}+\mu\,{\bf x}^{\prime}_{k}\right)=\left({\bf W}^{H}{\bf B}^{H}{\bf y}+\mu\,{\bf x}^{\prime}_{k}\right)$, allows writing (24) compactly as

Notice that multiplication by ${\bf F}$ corresponds to applying an image filter in the Fourier domain. Finally, notice also that the term ${\bf B}^{H}{\bf W}^{H}{\bf y}$ can be precomputed, as it doesn’t change during the algorithm.

The leading cost of each application of (30) will be either $O(n\log n)$ or the cost of the products by ${\bf W}^{H}$ and ${\bf W}$. For most tight frames used in image processing, these products correspond to direct and inverse transforms for which fast algorithms exist. For example, when ${\bf W}^{H}$ and ${\bf W}$ are the inverse and direct translation-invariant wavelet transforms, these products can be computed using the undecimated wavelet transform with $O(n\log n)$ total cost . Curvelets also constitute a Parseval frame for which fast $O(n\log n)$ implementations of the forward and inverse transform exist . Yet another example of a redundant Parseval frame is the complex wavelet transform, which has $O(n)$ computational cost , . In conclusion, for a large class of choices of ${\bf W}$, each iteration of the SALSA algorithm has $O(n\log n)$ cost.

III-C3 Missing Pixels: Image Inpainting

In the analysis prior case (TV-based), we have ${\bf A}={\bf B}$, where the observation matrix ${\bf B}$ models the loss of some image pixels. Matrix ${\bf A=B}$ is thus an $m\times n$ binary matrix, with $m<n$, which can be obtained by taking a subset of rows of an identity matrix. Due to its particular structure, this matrix satisfies ${\bf B}{\bf B}^{T}={\bf I}$. Using this fact together with the SMW formula leads to

Since ${\bf B}^{T}{\bf B}$ is equal to an identity matrix with some zeros in the diagonal (corresponding to the positions of the missing observations), the matrix in (31) is diagonal with elements either equal to $1/(\mu+1)$ or $1/\mu$. Consequently, (24) corresponds simply to multiplying $({\bf B}^{H}{\bf y}+\mu{\bf x}^{\prime}_{k})$ by this diagonal matrix, which is an $O(n)$ operation.

In the synthesis prior case, we have ${\bf A}={\bf BW}$, where ${\bf B}$ is the binary sub-sampling matrix defined in the previous paragraph. Using the SMW formula yet again, and the fact that ${\bf B}{\bf B}^{T}={\bf I}$, we have

As noted in the previous paragraph, ${\bf B}^{T}{\bf B}$ is equal to an identity matrix with zeros in the diagonal (corresponding to the positions of the missing observations), i.e., it is a binary mask. Thus, the multiplication by ${\bf W}^{H}{\bf B}^{H}{\bf B}{\bf W}$ corresponds to synthesizing the image, multiplying it by this mask, and computing the representation coefficients of the result. In conclusion, the cost of (24) is again that of the products by ${\bf W}$ and ${\bf W}^{H}$, usually $O(n\log n)$.

III-C4 Partial Fourier Observations: MRI Reconstruction.

The final case considered is that of partial Fourier observations, which is used to model magnetic resonance image (MRI) acquisition , and has been the focus of much recent interest due to its connection to compressed sensing . In the TV-regularized case, the observation matrix has the form ${\bf A}={\bf B}{\bf U}$, where ${\bf B}$ is an $m\times n$ binary matrix, with $m<n$, similar to the one in the missing pixels case (it is formed by a subset of rows of an identity matrix), and ${\bf U}$ is the DFT matrix. This case is similar to (32), with ${\bf U}$ and ${\bf U}^{H}$ instead of ${\bf W}$ and ${\bf W}^{H}$, respectively. The cost of (24) is again that of the products by ${\bf U}$ and ${\bf U}^{H}$, i.e., $O(n\log n)$ if we use the FFT.

In the synthesis case, the observation matrix has the form ${\bf A}={\bf B}{\bf U}{\bf W}$. Clearly, the case is again similar to (32), but with ${\bf UW}$ and ${\bf W}^{H}{\bf U}^{H}$ instead of ${\bf W}$ and ${\bf W}^{H}$, respectively. Again, the cost of (24) is $O(n\log n)$, if the FFT is used to compute the products by ${\bf U}$ and ${\bf U}^{H}$ and fast frame transforms are used for the products by ${\bf W}$ and ${\bf W}^{H}$.

IV Experiments

In this section, we report results of experiments aimed at comparing the speed of SALSA with that of the current state of the art methods (all of which are freely available online): TwISTAvailable at http://www.lx.it.pt/~bioucas/code/TwIST_v1.zip , SpaRSAAvailable at http://www.lx.it.pt/~mtf/SpaRSA/ , and FISTAAvailable at http://iew3.technion.ac.il/~becka/papers/wavelet_FISTA.zip . We consider three standard and often studied imaging inverse problems: image deconvolution (using both wavelet and TV-based regularization); image restoration from missing samples (inpainting); image reconstruction from partial Fourier observations, which (as mentioned above) has been the focus of much recent interest due to its connection with compressed sensing and the fact that it models MRI acquisition . All experiments were performed using MATLAB for Windows XP, on a desktop computer equipped with an Intel Pentium-IV $3.0$ GHz processor and $1.5$GB of RAM. To compare the speed of the algorithms, in a way that is as independent as possible from the different stopping criteria, we first run SALSA and then the other algorithms until they reach the same value of the objective function. The value of $\mu$ for fastest convergence was found to differ (though not very much) in each case, but a good rule of thumb, adopted in all the experiments, is $\mu=0.1\tau$.

We consider five benchmark deblurring problems , summarized in Table I, all on the well-known Cameraman image. The regularizer is \phi(\mbox{\boldmath\beta})=\|\mbox{\boldmath\beta}\|_{1}, thus \mbox{\boldmath\Psi}_{\tau\phi} is an element-wise soft threshold. The blur operator ${\bf B}$ is applied via the FFT. The regularization parameter $\tau$ is hand tuned in each case for best improvement in SNR, so that the comparison is carried out in the regime that is relevant in practice. Since the restored images are visually indistinguishable from those obtained in , and the SNR improvements are also very similar, we simply report computation times.

In the first set of experiments, ${\bf W}$ is a redundant Haar wavelet frame with four levels. The CPU times taken by each of the algorithms are presented in Table II. In the second set of experiments, ${\bf W}$ is an orthogonal Haar wavelet basis; the results are reported in Table III. To visually illustrate the relative speed of the algorithms, Figures 1 and 2 plot the evolution of the objective function (see Eq. (1)), versus time, in experiments $1$, $2$B, and $3$A, for redundant and orthogonal wavelets, respectively.

IV-B Image Deblurring with Total Variation

The same five image deconvolution problems listed in Table I were also addressed using total variation (TV) regularization (more specifically, the isotropic discrete total variation, as defined in ). The corresponding Moreau proximal mapping is computed using $5$ iterations of Chambolle’s algorithm .

The CPU times taken by SALSA, TwIST, SpaRSA, and FISTA are listed in Table IV. The evolutions of the objective functions (for experiments $1$, $2$B, and $3$A) are plotted in Figure 3.

We can conclude from Tables II, III, and IV that, in image deconvolution problems, both with wavelet-based and TV-based regularization, SALSA is always clearly faster than the fastest of the other competing algorithms.

IV-C MRI Image Reconstruction

We consider the problem of reconstructing the $128\times 128$ Shepp-Logan phantom (shown in Figure 4) from a limited number of radial lines (22, in our experiments, as shown in Figure 4) of its 2D discrete Fourier transform. The projections are also corrupted with circular complex Gaussian noise, with variance $\sigma^{2}\ =\ 0.5\times 10^{-3}$. We use TV regularization (as described in Subsection IV-B), with the corresponding Moreau proximal mapping implemented by $40$ iterations of Chambolle’s algorithm .

Table V shows the CPU times, while Figure 5 plots the evolution of the objective function over time. Figure 4 shows the estimate obtained using SALSA (the others are, naturally, visually indistinguishable). Again, we may conclude that SALSA is considerably faster than the other three algorithms, while achieving comparable values of mean squared error of the reconstructed image.

IV-D Image Inpainting

Finally, we consider an image inpainting problem, as explained in Section III-C. The original image is again the Cameraman, and the observation consists in loosing $40\%$ of its pixels, as shown in Figure 6. The observations are also corrupted with Gaussian noise (with an SNR of $40$ dB). The regularizer is again TV implemented by $20$ iterations of Chambolle’s algorithm.

The image estimate obtained by SALSA is shown in Figure 6, with the original also shown for comparison. The estimates obtained using TwIST and FISTA were visually very similar. Table VI compares the performance of SALSA with that of TwIST and FISTA and Figure 7 shows the evolution of the objective function for each of the algorithms. Again, SALSA is considerably faster than the alternative algorithms.

V Conclusions

References