An Augmented Lagrangian Approach to the Constrained Optimization Formulation of Imaging Inverse Problems

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 deconvolution, ${\bf B}$ is the matrix representation of a convolution operator. This type of 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 .

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

The problem of estimating ${\bf x}$ from ${\bf y}$ is called a linear inverse problem (LIP); for most scenarios of practical interest, this is an ill-posed LIP (IPLIP), i.e., matrix ${\bf B}$ is singular and/or very ill-conditioned. Consequently, this IPLIP requires some sort of regularization (or prior information, in Bayesian inference terms). One way to regularize the problem of estimating ${\bf x}$, given ${\bf y}$, consists in a constrained optimization problem of the form

In recent years, an explosion of interest in problems of the form (1) was sparked by the emergence of compressive sensing (CS) , . The theory of CS provides conditions (on matrix ${\bf B}$ and the degree of sparseness of the original ${\bf x}$) under which a solution of (1), for $\phi({\bf x})=\|{\bf x}\|_{1}$, is an optimal (in some sense) approximation to the “true” ${\bf x}$.

I-B Analysis and Synthesis Formulations

This formulation is referred to as the synthesis approach , , since it is based on a synthesis equation: ${\bf x}$ is synthesized from its representation coefficients ({\bf x}={\bf W}\mbox{\boldmath\beta}) which are the object of the estimation criterion. Naturally, the estimate of ${\bf x}$ is \widehat{\bf x}={\bf W}\widehat{\mbox{\boldmath\beta}}. Of course, (2) has the form (1) with ${\bf BW}$ replacing ${\bf B}$.

An alternative formulation applies a regularizer directly to the unknown image, leading to criteria of the form (1), usually called analysis approaches, since they are based on a regularizer that analyzes the image itself, rather than the coefficients of a representation thereof. Arguably, the best known and most often used regularizer in analysis approaches to image restoration is the total variation (TV) norm , .

Wavelet-based analysis approaches are also possible and have the form

where ${\bf P}$ is some linear operator (a matrix) . In this paper, we always assume that ${\bf P}$ is the analysis operator associated with 1-tight (Parseval) frame, thus ${\bf P}^{H}{\bf P}={\bf I}$ .

I-C Previous Algorithms

If the regularizers are convex, problems (1)–(3) are convex, but the very high dimension (at least $\geq 10^{4}$, often $\gg 10^{5}$) of ${\bf x}$ and ${\bf y}$ precludes the direct application of off-the-shelf optimization algorithms. This difficulty is further amplified by the following fact: for any problem of non-trivial dimension, matrices ${\bf B}$, ${\bf W}$, or ${\bf P}$ 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 these matrices (or their conjugate transposes, denote by $(\cdot)^{H}$) can be computed 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 transform algorithm . Similarly, if ${\bf B}$ represents a convolution, products by ${\bf B}$ or ${\bf B}^{H}$ can be performed with the help of the fast Fourier transform (FFT). These facts have stimulated the development of special purpose methods, in which the only operations involving matrices are matrix-vector products.

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 . For some choices of $\phi$, the corresponding \mbox{\boldmath\Psi}_{\tau\phi} have well known closed forms. For example, if $\phi({\bf x})\equiv\|{\bf x}\|_{1}$, then \mbox{\boldmath\Psi}_{\tau\phi}({\bf y})=\mbox{soft}({\bf y},\tau), where $\mbox{soft}(\cdot,\tau)$ denotes the component-wise application of the soft-threshold function $y\mapsto\mbox{sign}(y)\max\{|y|-\tau,0\}$.

In ,, we proposed a new algorithm called split augmented Lagrangian shrinkage algorithm (SALSA), to solve unconstrained optimization problems of the form (4) based on variable splitting , . The idea is to transform the unconstrained problem (4) into a constrained one via a variable splitting “trick”, and then attack this constrained problem using an augmented Lagrangian (AL) method . AL is known to be equivalent to the Bregman iterations 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 more standard and elementary tool (covered in most optimization textbooks). On several benchmark experiments (namely image deconvolution, recovery of missing pixels, and reconstruction from partial Fourier observations) using either frame-based or TV-based regularization, SALSA was found to be faster than the previous state-of-the-art methods FISTA , TwIST , and SpaRSA .

Although it is usually easier to solve an unconstrained problem than a constrained one, formulation (1) has an important advantage: parameter $\varepsilon$ has a clear meaning (it is proportional to the noise standard deviation) and is much easier to set than parameter $\tau$ in (4). Of course, one may solve (1) by solving (4) and searching for the “correct” value of $\tau$ that makes (4) equivalent to (1). Clearly, this is not efficient, as it involves solving many instances of (4). Obtaining fast algorithms for solving (1) is thus an important research front.

The Bregman iterative algorithm (BIA) was recently proposed to solve (1) with $\varepsilon=0$, but is not directly applicable when $\varepsilon>0$ . To deal with the case of $\varepsilon>0$, it was suggested that the BIA for $\varepsilon=0$ is used and stopped when $\|{\bf B}{\bf x}-{\bf y}\|_{2}\leq\varepsilon$ , . Clearly, that approach is not guaranteed to find a good solution, since it depends strongly on the initialization; e.g., if the algorithm starts at a feasible point, it will immediately stop, although the point may be far from a minimizer of $\phi$.

I-D Proposed Approach

In this paper, we introduce an algorithm for solving optimization problems of the form (1). The original constrained problem (1) is transformed into an unconstrained one by adding the indicator function of the feasible set, the ellipsoid $\{{\bf x}:\|{\bf Bx}-{\bf y}\|\leq\varepsilon\}$, to the objective in (1). The resulting unconstrained problem is then transformed into a different constrained problem, by the application of a variable splitting operation; finally, the obtained constrained problem is dealt with using the alternating direction method of multipliers (ADMM) , , , which belongs to the family of augmented Lagrangian (AL) techniques . Since (as SALSA), the proposed method uses variable splitting and AL optimization, we call it C-SALSA (for constrained-SALSA).

C-SALSA is experimentally shown to efficiently solve image recovery problems, such as MRI reconstruction from CS-type partial Fourier observations using TV regularization, and image deblurring using wavelet-based or TV regularization, faster than SPGL1 and NESTA.

I-E Organization of the Paper

The paper is organized as follows. Section II describes the basic ingredients of C-SALSA: variable splitting, augmented Lagrangians, and the ADMM. Section III contains the derivation leading to C-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

The rationale behind VS is that it may be easier to solve the constrained problem (7) than it is to solve its equivalent unconstrained counterpart (6).

VS has been recently used in several image processing applications. A VS method was used in to obtain a fast algorithm for TV-based restoration. Variable splitting was also used in to handle problems involving compound regularizers. In and , the constrained problem (7) is attacked using 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 (8) to approach that of (7), which in turn is equivalent to (6). The rationale of these methods is that each step of this alternating minimization may be much easier than the original unconstrained problem (6). The drawback is that as $\alpha$ grows, the intermediate minimization problems become increasingly ill-conditioned, thus causing numerical problems .

A similar VS 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 , which is known to be equivalent to the augmented Lagrangian method , , .

II-B Augmented Lagrangian

Consider the constrained optimization problem with linear equality constraints

Set $k=0$, choose $\mu>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 Az}_{k+1}-{\bf b})

until some stopping criterion is satisfied.

It is also possible (and even recommended) to update the value of $\mu$ in each iteration , . 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 (9).

After a straightforward complete-the-squares procedure, the terms added to $E({\bf z})$ in the augmented Lagrangian {\cal L}_{A}({\bf z},\mbox{\boldmath\lambda}_{k},\mu) can be written as a single quadratic term (plus a constant independent of ${\bf z}$, thus irrelevant to 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 Az}_{k+1}-{\bf b})$

until some stopping criterion is satisfied.

II-C ALM/MM for Variable Splitting and ADMM

The constrained problem (7) can be written as (9) by defining $E({\bf z})\equiv f_{1}({\bf u})+f_{2}({\bf v})$ and setting

With these definitions in place, Steps 3 and 4 of the ALM/MM (version II) become

The minimization problem yielding $\left({\bf u}_{k+1},{\bf v}_{k+1}\right)$ is not trivial since, in general, it involves a non-separable quadratic term and possibly non-smooth terms. A natural approach is to use a non-linear block-Gauss-Seidel (NLBGS) technique which alternates between minimizing with respect to ${\bf u}$ and ${\bf v}$ while keeping the other fixed. Of course this raises several questions: for a given ${\bf d}_{k}$, how much computational effort should be spent in this problem? Does the NLBGS procedure converge? The simplest answer to these questions is given in the form of the so-called alternating direction method of multipliers (ADMM) , , , which is simply an ALM/MM in which only one NLBGS step is performed in each outer iteration.

Set $k=0$, choose $\mu>0$, ${\bf v}_{0}$, ${\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})$

until some stopping criterion is satisfied.

For later reference, we now recall the theorem by Eckstein and Bertsekas [23, Theorem 8] in which convergence of (a generalized version of) ADMM is shown.

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

Notice that the ADMM as defined above (if each step is implemented exactly) generates sequences $\{{\bf u}_{k}\}$, $\{{\bf v}_{k}\}$, and $\{{\bf d}_{k}\}$ that satisfy the conditions in Theorem 1 in a strict sense (i.e., with $\eta_{k}=\nu_{k}=0$). The remaining key condition for convergence is then that ${\bf G}$ has full column rank. One of the important corollaries of this theorem is that it is not necessary to exactly solve the minimizations in lines 3 and 4 of ADMM; as long as the sequence of errors are absolutely summable, convergence is not compromised.

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

II-D A Variant of ADMM

Consider a generalization of problem (6), where instead of two functions, there are $J$ functions, that is,

Moreover, the fact that $f_{1}=0$ turns Step 3 of the algorithm into a simple quadratic minimization problem, which has a unique solution if ${\bf G}$ has full column rank:

where \mbox{\boldmath\zeta}_{k}={\bf v}_{k}+{\bf d}_{k} (and, naturally, \mbox{\boldmath\zeta}_{k}^{(j)}={\bf u}_{k}^{(j)}+{\bf d}_{k}^{(j)}) and the second equality results from the particular structure of ${\bf G}$ in (14).

Furthermore, our particular way of mapping problem (13) into problem (6) allows decoupling the minimization in Step 4 of ADMM into a set of $J$ independent ones. In fact,

Clearly, the minimizations with respect to ${\bf u}^{(1)},\dots,{\bf u}^{(J)}$ are decoupled, thus can be solved separately, leading to

Since this algorithm is exactly an ADMM, and since all the functions $g_{j}$, for $j=1,...,J$, are closed, proper, and convex, convergence is guaranteed if ${\bf G}$ has full column rank. Actually, this full column rank condition is also required for the inverse in (II-D) to exist. Finally, notice that the update equations in (32) can be written as

where the \mbox{\boldmath\Psi}_{g_{j}/\mu} are, by definition, the Moreau proximal mappings of $g_{1}/\mu,...,g_{J}/\mu$.

In summary, the variant of ADMM (herein referred to as ADMM-2) that results from the formulation just presented is described in the following algorithmic framework.

Set $k=0$, choose $\mu>0$, ${\bf v}_{0}^{(1)}$, …, ${\bf v}_{0}^{(J)}$, ${\bf d}_{0}^{(1)}$, …, ${\bf d}_{0}^{(J)}$.

do \mbox{\boldmath\zeta}_{k}^{(i)}={\bf v}_{k}^{(i)}+{\bf d}_{k}^{(i)}

{\bf u}_{k+1}=\biggl{[}\sum_{j=1}^{J}({\bf H}^{(j)})^{H}{\bf H}^{(j)}\biggr{]}^{-1}\sum_{j=1}^{J}\bigl{(}{\bf H}^{(j)}\bigr{)}^{H}\mbox{\boldmath\zeta}^{(j)}_{k}

do {\bf v}_{k+1}^{(i)}=\mbox{\boldmath\Psi}_{g_{i}/\mu}\left({\bf H}^{(i)}{\bf u}_{k+1}-{\bf d}^{(i)}_{k}\right)

${\bf d}_{k+1}^{(i)}={\bf d}_{k}^{(i)}-{\bf H}^{(i)}{\bf u}_{k+1}+{\bf v}_{k+1}^{(i)}$

until some stopping criterion is satisfied.

III Proposed Method

We now apply the algorithmic framework described in the previous section to the basic problem (1) (which includes (2) as a special case), as well as the analysis formulation (3).

For the constrained optimization problem (1), the feasible set is the ellipsoid

which is possibly infinite in some directions (since ${\bf B}$ may be singular). Problem (1) can be written as an unconstrained problem, with a discontinuous objective,

Notice that $E(\varepsilon,{\bf I},{\bf y})$ is simply a closed $\varepsilon$-radius Euclidean ball centered at ${\bf y}$.

Problem (35) has the form (13) with $J=2$ and

Instantiating ADMM-2 to this particular case requires the definition of the Moreau proximal maps associated with $g_{1}\equiv\phi$ and $g_{2}\equiv\iota_{E(\varepsilon,{\bf I},{\bf y})}$. Concerning $\phi$, the regularizer, we assume that \mbox{\boldmath\Psi}_{\tau\phi}(\cdot) (see (5)) can be computed efficiently. This is of course the case of $\phi({\bf x})\equiv\|{\bf x}\|_{1}$, for which \mbox{\boldmath\Psi}_{\tau\phi} is simply a soft threshold. If $\phi$ is the TV norm, we may use one the fast algorithms available to compute the corresponding denoising function , . The Moreau proximal map of $g_{2}\equiv\iota_{E(\varepsilon,{\bf I},{\bf y})}$ is defined as

which is obviously independent of $\mu$ and is simply the orthogonal projection of ${\bf s}$ on the closed $\varepsilon$-radius ball centered at ${\bf y}$:

We are now in a position to instantiate ADMM-2 for solving (35) (equivalently (1)). The resulting algorithm, which we call C-SALSA-1, is as follows.

Set $k=0$, choose $\mu>0$, ${\bf v}_{0}^{(1)}$, ${\bf v}_{0}^{(2)}$, ${\bf d}_{0}^{(2)}$, ${\bf d}_{0}^{(2)}$.

${\bf r}_{k}={\bf v}_{0}^{(1)}+{\bf d}_{0}^{(1)}+{\bf B}^{H}\left({\bf v}_{0}^{(2)}+{\bf d}_{0}^{(2)}\right)$

{\bf u}_{k+1}=\biggl{(}{\bf I}+{\bf B}^{H}{\bf B}\biggr{)}^{-1}{\bf r}_{k}

{\bf v}_{k+1}^{(1)}=\mbox{\boldmath\Psi}_{\phi/\mu}\left({\bf u}_{k+1}-{\bf d}^{(1)}_{k}\right)

{\bf v}_{k+1}^{(2)}=\mbox{\boldmath\Psi}_{\iota_{E(\varepsilon,{\bf I},{\bf y})}}\left({\bf B}{\bf u}_{k+1}-{\bf d}^{(2)}_{k}\right)

${\bf d}_{k+1}^{(1)}={\bf d}_{k}^{(1)}-{\bf u}_{k+1}+{\bf v}_{k+1}^{(1)}$

${\bf d}_{k+1}^{(2)}={\bf d}_{k}^{(2)}-{\bf B}{\bf u}_{k+1}+{\bf v}_{k+1}^{(2)}$

until some stopping criterion is satisfied.

The issue of how to efficiently solve the linear system of equations in line 4 of C-SALSA-1 will be addressed in Subsection III-C.

Convergence of C-SALSA-1 is guaranteed by Theorem 1 since it is an instance of ADMM with

which is a full column rank matrix, and both $\phi$ and $\iota_{E(\varepsilon,{\bf I},{\bf y})}$ are closed, proper, convex functions.

Finally, notice that to apply C-SALSA-1 to problem (2) we simply have to replace ${\bf B}$ with ${\bf BW}$.

III-B Problem (3)

Problem (3) can also be written as an unconstrained problem

The resulting ADMM algorithm, called C-SALSA-2, is similar to C-SALSA-1, with only a few minor differences.

Set $k=0$, choose $\mu>0$, ${\bf v}_{0}^{(1)}$, ${\bf v}_{0}^{(2)}$, ${\bf d}_{0}^{(2)}$, ${\bf d}_{0}^{(2)}$.

${\bf r}_{k}={\bf P}\left({\bf v}_{k}^{(1)}+{\bf d}_{k}^{(1)}\right)+{\bf B}^{H}\left({\bf v}_{k}^{(2)}+{\bf d}_{k}^{(2)}\right)$

{\bf u}_{k+1}\!=\!\biggl{(}{\bf P}^{H}{\bf P}+{\bf B}^{H}{\bf B}\biggr{)}^{-1}{\bf r}_{k}

{\bf v}_{k+1}^{(1)}=\mbox{\boldmath\Psi}_{\phi/\mu}\left({\bf P}{\bf u}_{k+1}-{\bf d}^{(1)}_{k}\right)

{\bf v}_{k+1}^{(2)}=\mbox{\boldmath\Psi}_{\iota_{E(\varepsilon,{\bf I},{\bf y})}}\left({\bf B}{\bf u}_{k+1}-{\bf d}^{(2)}_{k}\right)

${\bf d}_{k+1}^{(1)}={\bf d}_{k}^{(1)}-{\bf P}{\bf u}_{k+1}+{\bf v}_{k+1}^{(1)}$

${\bf d}_{k+1}^{(2)}={\bf d}_{k}^{(2)}-{\bf B}{\bf u}_{k+1}+{\bf v}_{k+1}^{(2)}$

until some stopping criterion is satisfied.

In this paper, we assume that ${\bf P}$ is the analysis operator of a 1-tight (Parseval) frame, thus ${\bf P}^{H}{\bf P}={\bf I}$ and line 4 of C-SALSA-2 is similar to line 4 of C-SALSA-1:

The issue of how to efficiently solve this linear system will be addressed in the next subsection.

Since both $\phi$ and $\iota_{E(\varepsilon,{\bf I},{\bf y})}$ are closed, proper, convex functions, convergence of C-SALSA-2 holds (by Theorem 1) if

is a full column rank matrix. This is of course true if ${\bf P}$ is itself a full column rank matrix, which is the case if ${\bf P}$ is the analysis operator of a tight frame .

III-C Solving (49)

As mentioned in Subsection I-C, in most imaging problems of interest, it may not be feasible to explicitly form matrix ${\bf B}$. This might suggest that it is not easy, or even feasible, to compute the inverse of $\left({\bf I}+{\bf B}^{H}{\bf B}\right)$. However, as shown next, in a number of problems of interest, this inverse can be computed very efficiently with $O(n\log n)$ cost.

In the case of analysis formulations of the form (1) or (3) to image deconvolution problems, matrix ${\bf B}$ represents a 2D convolution. Consequently, matrix ${\bf B}$ can be factorized as ${\bf B=U}^{H}{\bf DU}$, where ${\bf U}$ is the unitary matrix (${\bf U}^{H}={\bf U}^{-1}$) representing the discrete Fourier transform (DFT) and ${\bf D}$ is diagonal. Thus,

where $|{\bf D}|^{2}$ is the matrix with squared absolute values of the entries of ${\bf D}$. Since $|{\bf D}|^{2}+{\bf I}$ is diagonal, its inversion cost is $O(n)$. Products by ${\bf U}$ and ${\bf U}^{H}$ have $O(n\log n)$ cost, using the FFT algorithm.

III-C2 Deconvolution with Synthesis Formulation

In this case, as seen in Section I-B, we have ${\bf BW}$ instead of ${\bf B}$, and even if ${\bf B}$ is a convolution, ${\bf BW}$ is not diagonalizable by the DFT. To sidestep this difficulty, we assume that ${\bf W}$ contains a 1-tight (Parseval) frame (i.e., ${\bf W\,W}^{H}={\bf I}$). Using the Sherman -Morrison- Woodbury (SMW) matrix inversion lemma,

thus line 4 of C-SALSA-1 and C-SALSA-2 can be written as

Since ${\bf B}$ is a convolution, ${\bf B=U}^{H}{\bf DU}$, thus multiplying by ${\bf F}$ corresponds to applying an image filter in the Fourier domain

which has $O(n\log n)$ cost, since all the matrices in ${\bf D^{*}}\left(|{\bf D}|^{2}+{\bf I}\right)^{-1}{\bf D}$ are diagonal and the products by ${\bf U}$ and ${\bf U}^{H}$ are carried out via the FFT. The cost of (52) will thus 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, there are fast $O(n\log n)$ algorithms to compute the products by ${\bf W}^{H}$ and ${\bf W}$ . For example, in the case of 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 formulation of this problem, the observation matrix ${\bf B}$ models the loss of some image pixels; it 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}^{H}={\bf I}$. Using this fact together with the SMW formula leads to

Since ${\bf B}^{H}{\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 (55) is diagonal with elements either equal to $1$ or $1/2$. Consequently, line 4 of C-SALSA-1 and C-SALSA-2 corresponds to multiplying this diagonal matrix by ${\bf r}_{k}$, obviously with $O(n)$ cost.

In the frame-based synthesis formulation, we have ${\bf BW}$ instead of ${\bf B}$. Using the SMW formula yet again, and the facts that ${\bf B}{\bf B}^{H}={\bf I}$ and ${\bf WW}^{H}={\bf I}$, we have

As noted in the previous paragraph, ${\bf A}^{H}{\bf A}$ is equal to an identity matrix with zeros in the diagonal, i.e., a binary mask. Thus, the multiplication by ${\bf W}^{H}{\bf A}^{H}{\bf A}{\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 line 4 of C-SALSA-1 and C-SALSA-2 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)

Finally, we consider the case of partial Fourier observations, which is used to model MRI acquisition and has been the focus of recent interest due to its connection to compressed sensing , . In the analysis formulation, ${\bf B}={\bf M}{\bf U}$, where ${\bf M}$ is an $m\times n$ binary matrix ($m<n$) again, formed by a subset of rows of the identity, and ${\bf U}$ is the DFT matrix. Due to its particular structure, matrix ${\bf M}$ satisfies ${\bf M}{\bf M}^{H}={\bf I}$; this fact together with the matrix inversion lemma leads to

where ${\bf M}^{H}{\bf M}$ is equal to an identity with some zeros in the diagonal. Consequently, the cost of line 4 of C-SALSA-1 and C-SALSA-2 is again that of the products by ${\bf U}$ and ${\bf U}^{H}$, i.e. $O(n\log n)$ using the FFT.

In the synthesis case, the observation matrix has the form ${\bf M}{\bf U}{\bf W}$. Clearly, the case is again similar to (56), but with ${\bf UW}$ and ${\bf W}^{H}{\bf U}^{H}$ instead of ${\bf W}$ and ${\bf W}^{H}$, respectively. Again, the cost of line 4 of C-SALSA-1 and C-SALSA-2 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}$.

III-D Computational Complexity

As shown in the previous section, the cost of line 4 of C-SALSA-1 and C-SALSA-2 is $O(n\log n)$. The other lines of the algorithms simply involve: (a) matrix-vector products involving ${\bf B}$, ${\bf W}$, ${\bf P}$, or their conjugate transposes, which have $O(n\log n)$ cost; (b) vector additions, which have $O(n)$ cost; and (c) the computation of the Moreau proximal maps (lines 5 and 6 of C-SALSA-1 and C-SALSA-2). In the case of the projections on a ball (line 6), it is clear from (42) that the cost is $O(n)$.

Finally, we consider the computational cost of the Moreau proximal map of the regularizer $\phi$ (line 5 of C-SALSA-1 and C-SALSA-2). In some cases, this map can be computed exactly in closed form; for example, if $\phi({\bf x})\equiv\|{\bf x}\|_{1}$, then \mbox{\boldmath\Psi}_{\tau\phi} is simply a soft threshold and the cost is $O(n)$. In other cases, the Moreau proximal map does not have a closed form solution; for example, if $\phi({\bf x})\equiv\mbox{TV}({\bf x})$, the corresponding \mbox{\boldmath\Psi}_{\tau\phi} has to be computed using one of several available iterative algorithms , . Most of these iterative algorithms can be implemented with $O(n)$ cost, although with a factor that depends on the number of iterations. In our implementation of C-SALSA we use Chambolle’s algorithm .

In summary, for a wide choice of regularizers and frame representations, the C-SALSA algorithms have $O(n\log n)$ computational complexity.

IV Experiments

In this section, we report results of experiments aimed at comparing the speed of C-SALSA with that of the current state of the art methods (that are freely available online): SPGL1Available at http://www.cs.ubc.ca/labs/scl/spgl1 , and NESTAAvailable at http://www.acm.caltech.edu/~nesta .

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, on a Windows XP desktop computer with an Intel Pentium-IV $3.0$ GHz processor and $1.5$ GB of RAM. The number of calls to the operators ${\bf B}$ and ${\bf B}^{H}$, the number of iterations, CPU times, and MSE values presented are the averages values over 10 runs of each experiment. The number of calls reported for each experiment is the average over the 10 instances, with the minimum and maximum indicated in the parentheses. Since the stopping criteria of the implementations of the available algorithms differ, to compare the speed of the algorithms in a way that is as independent as possible from these criteria, the experimental protocol that we followed was the following: we first run one of the algorithms with its stopping criterion, and then run C-SALSA until the constraint in (1) is satisfied and the MSE of the estimate is below that obtained by the other algorithms.

The value of $\varepsilon$ in (1) used in all cases was $\sqrt{m+8\sqrt{m}}\sigma$, where $m$ is the number of observations, and $\sigma$ is the noise standard deviation. The parameter $\mu$ was hand-tuned for fastest convergence.

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. We compare C-SALSA against SPGL1 and NESTA in the synthesis case, and against only NESTA in the analysis case, since SPGL1 is hardwired with $\|{\bf x}\|_{1}$ as the regularizer, and not $\|{\bf Px}\|_{1}$. Since the restored images are visually indistinguishable from those obtained in , and the SNR improvements are also very similar, we simply compare the speed of the algorithms, that is, the number of calls to the operators ${\bf B}$ and ${\bf B}^{H}$, the number of iterations, and the computation time.

In the first set of experiments, ${\bf W}$ is a redundant Haar wavelet frame with four levels. For the synthesis case, the CPU times taken by each of the algorithms are presented in Table II. Table III presents the corresponding results for the case with the analysis prior. In the second set of experiments, ${\bf W}$ is an orthogonal Haar wavelet basis; the results are reported in Table IV for the synthesis case, and in Table V for the analysis case. To visually illustrate the relative speed of the algorithms, Figure 1 plots the evolution of the constraint $\|{\bf Bu}_{k}-{\bf y}\|$, versus time, in experiments $1$, for the synthesis prior case, with redundant wavelets.

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 .

Table VI compares the performance of C-SALSA and NESTA, in terms of speed. The evolutions of the objective function and the constraint for experiment $1$ are plotted in Figure 2.

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

IV-C MRI Image Reconstruction

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

Table VII shows the number of calls, number of iterations, and CPU times, while Figure 4 plots the evolution of the objective function and constraint over time. Figure 3 shows the estimate obtained using C-SALSA (the estimate NESTA is, naturally, visually indistinguishable). Again, we may conclude that C-SALSA is faster than NESTA, while achieving comparable values of mean squared error of the reconstructed image.

A related example that we will consider here is the reconstruction of images composed of random squares, from their partial Fourier measurements, with TV regularization (see , section $6.4$). The dynamic range of the signals (the amplitude of the squares) varies from 20 dB to 80 dB. The size of each image is $128\times 128$, the number of radial lines in the DFT measurement mask is $27$ (corresponding to $m/n\approx 0.2$), and the Gaussian noise standard deviation is $\sigma=0.1$.

Figure 5 shows the original image with a dynamic range of $40$ dB and the estimate obtained using C-SALSA. Figure 6 shows the evolution over time of the objective and the error constraint for C-SALSA and NESTA, while Table VIII compares the two algorithms with respect to the number of calls to ${\bf A}$ and ${\bf A}^{H}$, number of iterations, CPU time, and MSE obtained, over 10 random trials. It is clear from Table VIII that C-SALSA uses considerably fewer calls to the operators ${\bf A}$ and ${\bf A}^{H}$ than NESTA.

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 losing $40\%$ of its pixels, as shown in Figure 7. The observations are also corrupted with Gaussian noise (with an SNR of $40$ dB). The regularizer is again TV, implemented by $10$ iterations of Chambolle’s algorithm.

The image estimate obtained by C-SALSA is shown in Figure 7, with the original also shown for comparison. The estimate obtained using NESTA was visually very similar. Table IX compares the performance of the two algorithms, and Figure 8 shows the evolution of the objective function for each of them.

V Conclusions

References