Previous Up Next

Chapter 3  Reconstruction Methods

While in a majority of biomedical modalities, images are produced directly, in MRI they are only obtained after a reconstruction process. In this chapter, we present the current approaches to MRI reconstruction. We recall in Section 3.1 the continuous model of the data-collection process in MRI and derive its discretized version. In Section 3.2, we review the different approaches that lead to linear reconstructions and we introduce the key concepts of the inverse problem formalism. This approach maps image reconstruction into an optimization problem with the possibility to impose a priori constraints to distinguish the solution from other possible candidates and improve reconstruction quality. We legitimate in Section 3.3 the use of sparsity-promoting priors in MRI and explain how they can be imposed via a proper regularization term. Finally, we review in Section 3.4 the algorithmic procedures that are theoretically capable of achieving the desired reconstructions while being suited to the practical constraints encountered in MRI.

3.1  Model for MRI

3.1.1  Physics

A radio-frequency pulse is emitted to initiate nuclear magnetic resonance (NMR). It excites the spins in a 2-D plane or a 3-D volume, depending of the type of acquisition format. After excitation, the excited spins behave as radio-frequency emitters and have their precessing frequency and phase modified depending on their positions. This is achieved thanks to the time-varying magnetic gradient fields that are applied during the relaxation, defining a trajectory k in the k-space domain. The modulated part of the signal received by a coil of sensitivity Si(r) is given by

mi(k) = 
 


2
Si(r)ρ(r)e−2jπk·rdr.     (1)

The signal ρ is referred to as object. This signal is proportional to the spin density, but might also depend upon other local characteristics. More details on the derivation of relation (1) are provided in Chapter 2.

For an array of R receiving coils with sensitivities denoted by S1SR and a k-space trajectory sampled at N points kn, we represent the measurements concatenated in a global RN× 1 vector

m=

m1,1,…,mN,1
,… 
m1,i,…,mN,i
,… 
m1,R,…,mN,R

.

3.1.2  Model for the Original Data

Spatial Discretization of the Object

From here on, we consider that the Fourier domain and, in particular, the sampling points kn, are scaled to make the Nyquist sampling interval unity. This can be done without any loss of generality if the space domain is scaled accordingly. Therefore, we model the object as a linear combination of pixel-domain basis functions ϕp that are shifted replicates of some generating function ϕ, so that

     
ρ=
 
p∈ ℤ2
 c[pp, with 
    (2)
ϕp(r) =
ϕ
rp
. 
    (3)

Given a sampled version of the coil sensitivity si[p], the sensitivity-weighted object is modeled by

Siρ=
 
p∈ ℤ2
 si[p]c[pp.     (4)

The standard implicit choice for ϕ is Dirac’s delta even if it is hardly justified from an approximation theoretic point of view. Different discretizations have been proposed, for example by Sutton et al. [38] with ϕ as a boxcar function or later by Delattre et al. [39] with B-splines. It is only recently [40] that the details have been worked out to get back the image for general ϕ that are non-interpolating, which is the case for instance for B-splines of degree greater than 1. The image to be reconstructed—i.e., the sampled version of the object ρ(p)—is obtained by filtering the coefficients c[p] with the discrete filter

P
ejω
=
 
h∈ ℤ2
 
^
ϕ
 

ω+2πh
,     (5)

where ϕ^ denotes the Fourier transform of ϕ.

Since a finite field of view (FOV) determines sets of coefficients c and si with a finite number M of elements, we handle them as a vectors c and si, keeping the discrete coordinates p as implicit indexing.

Wavelet Discretization

Due to sparsity properties that are discussed later in this chapter, it might be preferable to represent the object in terms of wavelet coefficients. In the wavelet formalism, some constraints apply to ϕ. It must be a scaling function that satisfies the properties for a multiresolution [41]. In that case, the wavelets can be defined as linear combinations of the ϕp and the object is equivalently characterized by its coefficients in the orthonormal wavelet basis. We refer to Mallat’s reference book [42] for a full review on wavelets. There exists a discrete wavelet transform (DWT) that bijectively maps the coefficients c to the wavelet coefficients w that represent the same object ρ in a continuous wavelet basis. In the rest of the chapter, we represent this DWT by the synthesis matrix W. Note that the matrix-vector multiplications c = Ww and w = W−1c have efficient filterbank implementations.

3.1.3  Matrix Representation of the Model

The data-formation model (1) and the object parameterization (4) are combined to model the measurement corresponding to every point kn sampled in k-space. Accordingly, the measurement vector m is related to the coefficients c through the linear relation

m = Ec,     (6)

where the MRI encoding matrix E is formed as

E = 
IRE0

diag
s1
, … , diag
sR

T,     (7)

with the symbol ⊗ standing for the Kronecker product, IR representing the R× R identity matrix, and E0 being the encoding matrix for the same MRI scan with a single homogeneous receiving coil

E0 = diag





^
ϕ
 
(2πk1),…,
^
ϕ
 
(2πkN)






v1,…,vN
T.     (8)

There, vn are vectors indexed like c such that vn[p] = exp(−2jπkn·p).

Due to the presence of noise and other scanner inaccuracies, the introduction of a new term b, accounting for an additive perturbation, makes the data-formation model

m = E c+b     (9)

more realistic. Equivalently, if the parameters of interest are the wavelet coefficients w, the model writes

m = M w+b     (10)

with M = EW.

In MRI, the major source of noise is a radio-frequency signal originating from the thermal motion in the object under investigation. When observed with a receiving array of coils, this noise presents non-negligible correlations across channels. In other terms, the R× R channel cross-correlation matrix Θ has non-null off-diagonal entries. Accordingly, the additive perturbation is generally modeled as the realization of a centered multivariate Gaussian process bN(0, Ψ) with covariance matrix Ψ=ΘIN.

3.2  Linear Solutions

The problem of imaging is to recover the M coefficients c (or equivalently w) from the N corrupted measurements m. In this section, we review the popular approaches that lead to reconstructions that depend linearly upon the observations. We show that they are functionally equivalent. Two of these approaches rely on a stochastic interpretation of the problem, where the matrices Ψ and Υ are the known covariance matrices of the noise b and the object c, respectively. The corresponding global variances are given by vn = Tr(Ψ)/N=Tr(Θ) and vs = Tr(Υ)/M. We define the normalized covariance matrices as Ψ0=Ψ/vn and Υ0=Υ/vs. Most linear solutions involve a balancing parameter λ which is necessarily positive and can be interpreted in terms of the signal-to-noise ratio λ−1=Tr(Υ)/Tr(Ψ)=M vs/(N vn).

3.2.1  Pseudoinverse

Depending on the scanner settings, the encoding matrix E is generally neither square nor invertible. In such cases, the Moore-Penrose pseudoinverse offers a solution to the reconstruction problem. The reconstruction matrix is then defined as

E = 
 
lim
є→ 0+
 

EHE+є IM
−1EH.     (11)

The Hermitian matrix transpose is used, denoted by the superscript H, because in MRI matrices have complex-valued entries. The problem of inverting a non-square matrix is tackled by considering the backprojected1 problem

EHm = EHEc,     (12)

because the matrix EHE is square.

Considering the singular value decomposition E=UΣVH, where Σ is an RN× M matrix whose diagonal entries are the singular values σn, one gets E=UH, with singular values

σn=

0 if σn= 0,
1/σn otherwise.

The major concern with pseudoinverse reconstruction resides in the propagation of noise. Indeed, very small but non-null singular values lead to drastic amplification of the corresponding noise components. This effect is quantified by the condition number, defined as κ(E)=maxn σn/minn σn. This number, which is greater or equal to 1, is also representative of the numerical challenge faced when inverting E. A linear inverse problem is termed “ill-conditioned” when the corresponding condition number is large. When the null space of E is not limited to {0}, the problem is said to be ill-posed.

The aim of regularized reconstruction schemes is to improve reconstruction with respect to the pseudoinverse approach by limiting the propagation of noise in the images.

It is remarkable that (12) rewrites as c(||mEc ||22)=0, with c standing for the gradient operator with respect to c. The More-Penrose pseudoinverse provides a least-squares solution c=Em to the reconstruction problem because it ensures EHm=EHEc. This least-squares solution makes sense when the noise term b is independent and identically distributed.

Instead, when the noise correlation matrix Ψ0 is available, this knowledge can be exploited using the weighted pseudoinverse

EX = 
 
lim
є→ 0+
 

EHXE+є IM
−1EHX     (13)

with the weighting matrix X=Ψ0. The interest of that type of solution is that it takes into account noise correlations and that it relies less on the noisier samples. Thanks to the relation EHXEEXm=EHXm, the weighted pseudoinverse provides a (weighted) least-squares solution.

3.2.2  Quadratic Regularization

The approach proposed by Phillips [43] and Twomey [44], for finite dimensional problems, and by Tikhonov [45], for infinite dimensional problems, defines the reconstruction as the minimization of the functional

⎪⎪
⎪⎪
m − E c  ⎪⎪
⎪⎪
X2 + λ ⎪⎪
⎪⎪
Rc ⎪⎪
⎪⎪
2,     (14)

where the notation ||· ||X with X positive-definite stands for a weighted norm such that ||v ||X2=vHXv. The functional is a trade-off between a fidelity term, which enforces consistency with the measurements, and a regularization term, which penalizes non-regular solutions with respect to the regularization matrix R. The tuning parameter λ balances the influence of these two terms. The role of the regularization term is to limit the amplification of noise that can be dramatic for ill-conditioned problems (in MRI, see for instance [46]). In practice, it is often designed with a derivation operator to favor smooth solutions. Similar to the weighted pseudoinverse solution, the weighting matrix can be chosen as X=Ψ0 yielding a reconstruction matrix that gives importance to the samples in inverse proportion to their level of noise. Another common choice is to take X diagonal such as to compensate for an inhomogeneous k-space sampling density [37]. This choice facilitates the reconstruction.

The minimization of a quadratic functional yields a linear solution. Indeed, by taking the gradient of the functional and setting it to zero, we find that the reconstruction matrix writes

FQUAD = 

EHXERHR
−1EHX.     (15)

3.2.3  Maximum a posteriori

Here, the reconstruction problem is tackled within a stochastic framework. The unknowns c and b are modeled as realizations of centered multivariate Gaussian distributions: cN(0, Υ) and bN(0, Ψ).

According to the numerical model (9), the measurements also follow a multivariate Gaussian distribution mN(0, EΥEH+Ψ).

The maximum a posteriori solution (MAP) c is the vector that maximizes the posterior distribution given the measurements m. Using Bayes’ theorem, the probability density function of the posterior distribution of c writes

p(c ∣ m) ∝ p(m ∣ cp(c).

In the present stochastic setting, the probability density function can be expanded in

p(c ∣ m) ∝ exp
⎪⎪
⎪⎪
mEc ⎪⎪
⎪⎪
Ψ2
exp
⎪⎪
⎪⎪
c ⎪⎪
⎪⎪
Υ2
.     (16)

Finally, the MAP solution is the vector c that minimizes the functional

⎪⎪
⎪⎪
mEc ⎪⎪
⎪⎪
Ψ02⎪⎪
⎪⎪
c ⎪⎪
⎪⎪
Υ02.     (17)

We introduced the normalized covariance matrices in the later expression in order to have the parameter λ, which is the inverse of the signal to noise ratio, appear explicitly.

Similarly to the previous approaches, the functional to be minimized is composed of quadratic terms. As a consequence, the solution is linear, characterized by the reconstruction matrix

FMAP = 

EHΨ0EΥ0
−1EHΨ0.     (18)

3.2.4  Linear Minimum Mean Squared Error estimator

The Gaussian model used in the MAP approach is hardly justified for true MRI images c. This assumption can be substituted by the constraint that the reconstruction is affine with respect to the measurements. Accordingly, we write the reconstructed image Fm+g.

To determine adequate parameters F and g, one can rely on the two first order statistics of the unknown data; that are, the expectation vectors c and b, and covariance matrices Υ and Ψ. According to the data-formation model (9), the expectation and covariance of the reconstruction error e=Fm+gc are given by

E
e
F(Ec+b) + g − c     (19)

and

E

e−E
e


e−E
e

H

= (FEI)Υ(FEI)H+FΨFH + E
e
E
e
H.     (20)

An unbiased reconstruction2 is obtained when g=cF(Ec+b). For the choice of F, one would reasonably like to minimize the variance of the reconstruction error. Given that the estimator is unbiased, the variance also corresponds to the expectation of the mean-square error. It is is given by the trace of the covariance matrix

Var
e
= Tr
(FEI)Υ(FEI)H
+Tr
FΨFH
.     (21)

Interestingly, this relation reveals two distinct contributions to the error.

The matrix F that minimizes the error variance, also referred to as mean-square error, can be computed using matrix calculus. Using the normalized covariance matrices, it writes

FMMSE = Υ0EH(0EHΨ0)−1.     (22)

3.2.5  Connections

First, Equations (15) and (18) show that quadratic regularization and MAP approaches are equivalent provided that X=Ψ0 and Υ0=RHR.3

Second, the three following equalities reveal the connection between MAP and LMMSE solutions, (18) and (22), in the case where both matrices Υ0 and Ψ0 are invertible:

EHΨ0−10EHΥ0−1Υ0EH=EHΨ0−10EHEHΨ0−1Ψ0
(EHΨ0−1EΥ0−1)Υ0EH=EHΨ0−1(0EHΨ0)
Υ0EH

0EHΨ0
−1
=

EHΨ0−1EΥ0−1
−1EHΨ0−1.

Last, the weighted pseudoinverse solution with X=Ψ0 corresponds to the other solutions in the limiting case where λ tends to 0. This is also the case for the regular Moore-Penrose pseudoinverse when the noise is independent and identically distributed; that is to say Ψ0=IRN/R. As already mentioned, the pseudoinverse solutions are only valid when noise propagation is negligible. This situation occurs with well-conditioned (κ(E)≈ 1) reconstruction problems that are largely overdetermined (MRN) and/or subject to very little noise (Tr(Υ)≫Tr(Ψ)).

About the invertibility of Ψ0 and Υ0:

There is no particular reason for Ψ0 to be singular. Most of the time, the correlation between pixels in the image are not modeled; this translates in a matrix Υ0 which is diagonal. When no signal is expected from some pixels of the image, (for instance outside a predetermined ROI), it could be tempting to set to 0 the corresponding entries in Υ0, resulting in a singular matrix. However, a reasonable problem setting would exclude such entries in the unknown vector c, restoring the invertibility of Υ0.

3.3  Non-Quadratic Regularizations

We just saw that the linear approaches to reconstruction can be derived from the solution of some optimization problems. The corresponding functionals were quadratic, yielding closed-form solutions. In this section, we consider other approaches that are popular in MRI and which involve non-quadratic regularization terms.

3.3.1  Inverse Problem Formalism

The solution c is defined as the minimizer of a cost function that involves two terms: the data fidelity F(b) and the regularization R(c) that penalizes undesirable solutions. This is summarized as

c= arg 
 
min
c
  F(m − Ec) + λ  R(c),     (23)

where the regularization parameter λ≥ 0 balances the two constraints. In MRI, the noise term b=mEc is usually assumed to be the realization of a Gaussian process with normalized covariance matrix Ψ0. From a Bayesian point of view, this justifies the choice F(b)=||b ||Ψ02=bHΨ0b as a proper log-likelihood term. A more practical motivation for this choice is that a quadratic fidelity term yields a simple closed-form gradient that greatly facilitates the design and performance of reconstruction algorithms.

When the k-space sampling is dense enough and the signal-to-noise ratio is high, the quadratic regularization terms (presented in the previous section) yield satisfying reconstructions. But, the constraints to reduce the scan duration favor setups with reduced SNR and k-space trajectories that present regions of low sampling density. In these situations where the reconstruction problem is more challenging, the reconstructed image can often be enhanced by the use of a more suitable regularization term R(c).

3.3.2  Total Variation

Total Variation (TV) was introduced as an edge-preserving denoising method by Rudin et al. [47]. It is now a very popular approach to tackle image enhancement problems.

The TV regularization term corresponds to the sum of the Euclidean norms of the gradient of the object. In practice, it is defined as R(w)=||∇c ||1. In this context, the operator returns pixelwise the ℓ2-norm of finite differences. The use of TV regularization is particularly appropriate for piecewise-constant objects such as the Shepp-Logan (SL) phantom used for simulations in tomography and MRI. Textured and noisy images exhibit a much larger total variation.

Sparsity-Promoting Regularization

Another popular idea is to exploit the fact that the object can be well represented by few non-zero coefficients (sparse representation) in an orthonormal basis of M functions φp. Formally, we write that

It is well-documented that typical MRI images admit sparse representation in bases such as wavelets or block DCT [6]. We illustrate this property in Figure 3.1.


Figure 3.1: Sparse approximation errors of a typical MRI brain image using different orthonormal transforms. The Mean Squared Error is represented as a function of the percentage of coefficients kept. Wavelet transforms achieve the best sparse approximations.

The ℓ1-norm is a good measure of sparsity with interesting mathematical properties (e.g., convexity). Thus, among the candidates that are consistent with the measurements, we favor a solution whose wavelet coefficients have a small ℓ1-norm. Specifically, the solution is formulated as

w= arg 
 
min
w
  C(w),     (24)

with

 C(w) = ⎪⎪
⎪⎪
m − M w ⎪⎪
⎪⎪
22⎪⎪
⎪⎪
w ⎪⎪
⎪⎪
1.     (25)

This is the general solution for wavelet-regularized inverse problems considered by [19] as well as by many other authors.

3.4  Algorithms

MRI gives rise to a large-scale inverse problem in the sense that the number of degrees of freedom—that is to say, the unknown pixel values—is large. Consequently, the matrices are generally too large to be stored in memory not to mention the fact that direct matrix multiplication involves too many operations.4 We summarize in this section the strategies that make the reconstruction in MRI feasible with reasonable computer requirements and acceptable computation times.

3.4.1  Matrix-Vector Multiplications

The matrix-vector multiplications y=E0x and y=E0Hx are two basic operations in MRI reconstruction. They can be implemented efficiently using the FFT algorithm. For non-Cartesian samples kn, the gridding method, based on FFT and interpolation, can provide accurate computations (see [48] for instance). Algorithms 1 and 2 describe the implementation of the operations y=E0x and y=E0Hx, respectively.

Algorithm 1   Matrix-Vector multiplication y=E0x, according to (8).
Inputs:
x, p1,…, pM, k1,…, kN, and ϕ^
Return: y
Algorithm 2   Matrix-Vector multiplication y=E0Hx, according to (8).
Inputs:
x, k1,…, kN, p1,…, pM, and ϕ^
Return: y

An interesting work by Wajer [49] identifies E0HE0 as a convolution matrix associated to the kernel

G[p] = 
N
n=1
 


^
ϕ
 
(2πkn)


2



 
exp
2jπ kn·p
.     (26)

When the kernel is precomputed for the lattice points belonging to the set S={ pqp∈FOV,  q∈FOV}, one can avoid the use of Algorithms 1 and 2. An efficient implementation of the operation y=E0HE0x, which uses zero-padded multidimensional FFTs, is described in Algorithm 3.

Algorithm 3   Matrix-Vector multiplication y=E0HE0x, according to (8) and (26).
Inputs:
x, S and G;
Return: y

Most of the time, in parallel MRI, the covariance matrices are block diagonal. In that case, they are sparse matrices and one can benefit from the related efficient memory storage and matrix operations. As already mentioned, Ψ0 is fully characterized by the channel cross-correlation matrix Θ0=Θ/vn such that Ψ0=Θ0IN. Its pseudoinverse or inverse is then given by Ψ0=Θ0IN. The matrix-vector multiplications with EHΨ0 and EHΨ0E are implemented as described in Algorithms 4 and 5, respectively.

Algorithm 4   Matrix-Vector multiplication y=EHΨ0x, according to (7).
Inputs:
x, s1,…, sR, Θ0 and E0H
Return: y
Algorithm 5   Matrix-Vector multiplication y=EHΨ0Ex, according to (7).
Inputs:
x, s1,…, sR, Θ0 and E0HE0
Return: y

3.4.2  Conjugate Gradient

The conjugate gradient method (CG) [50] is an iterative algorithm that is among the most efficient in solving large-scale linear problems Ac=b, characterized by symmetric and positive-definite matrices A. The only operations involving the matrix A are matrix-vector multiplications Ax. In parallel MRI, it is the method of reference [37] to perform linear reconstructions. The quadratic-regularized solution characterized by the reconstruction matrix in (15) is computed with CG solving the linear problem defined by the matrix A=EHXERHR and vector b=EHXm.

The idea of the method is to decompose the solution in a basis of mutually conjugate vectors; that is to say c=∑i αipi, with piHApj=0 for ij. At iteration i, the estimate is ci=∑jiαjpj and the corresponding residue writes ri=bAci. For the next direction, the choice pi+1=ri−∑ji(pjHAri)pj/||pj ||A ensures the conjugacy constraint. In this direction, the coefficient αi+1=Re(pi+1HAri)/||pi+1 ||A is optimal with respect to the cost C(c)=cHAccHbbHc. An efficient implementation of the method is described in Algorithm 6.

Algorithm 6   CG solving Ac=b with A symmetric and positive-definite.
Inputs:
A, b, and c0 (optional, default: c0=0)
Initialization:
r0=bAc0, p0=r0, and i=0
Repeat until desired tolerance is reached:
Return: ci

The CG algorithm theoretically converges within a finite number of iterations. In practice, this result is compromised by the propagation of round-off errors. In the context of MRI, the property of practical interest is the linear convergence rate achieved by CG. Indeed, the distance to the desired solution decreases as a power of the iteration number, with the convergence rate

0≤ r(A)=

κ(A)
−1

/

κ(A)
+1

<1.

When the condition number κ(A) is large, the rate r(A) gets close to the unity, characterizing a slower convergence. Using the weighted norm ||x ||A=√xHAx, the distance is upperbounded by

⎪⎪
⎪⎪
cic ⎪⎪
⎪⎪
A≤ 2⎪⎪
⎪⎪
c0c ⎪⎪
⎪⎪
Ar(A)i.     (27)

With the regular Euclidean distance, the bound is looser

⎪⎪
⎪⎪
cic ⎪⎪
⎪⎪
2 ≤ 2κ(A)⎪⎪
⎪⎪
c0c ⎪⎪
⎪⎪
2r(A)i.     (28)

3.4.3  Iteratively Reweighted Least-Squares

The Iteratively Reweighted Least-Squares algorithm (IRLS), which is also known as the positive form of half-quadratic minimization [51], can be used to compute the solutions defined as

c= arg 
 
min
c
⎪⎪
⎪⎪
mEc ⎪⎪
⎪⎪
X2 + λ⎪⎪
⎪⎪
Rc ⎪⎪
⎪⎪
pp.     (29)

In this context, the functional is strictly convex for p>1. This condition ensures the unicity of the minimizer.

The principle of IRLS is to design an upperbounding quadratic proxy for the regularization term, tailored to the neighborhood of ci. In practice, one chooses the functional

 Zi(c)=p/2 ⎪⎪
⎪⎪
Rc ⎪⎪
⎪⎪
Di2 + (1−p/2) ⎪⎪
⎪⎪
Rci ⎪⎪
⎪⎪
pp,     (30)

where Di is a diagonal matrix with entries |(Rci)n|p−2. It has the following desirable properties

  1. Zi(ci)=||Rci ||pp,
  2. c Zi(ci) = c ||Rc ||pp(c=ci),
  3. Zi(c)>||Rc ||pp for all cci and p<2,
  4. arg minc ||mEc ||X2 + λ Zi(c)=(EHXE+(λ p/2)RHDiR)−1EHXm.

An implementation of the IRLS is described in Algorithm 7.

Algorithm 7   IRLS solving c=arg minc||mEc ||X2 + λ||Rc ||pp.
Inputs:
A=EHXE, a=EHXm, R, p, λ, and c0 Repeat until desired tolerance is reached: Return: ci

Let us remember that for p≤ 1 the minimization problem might not admit a unique solution. When the minimizer c is unique, it is also the unique fixed-point of the algorithm. As long as p<2, the sequence of functional values C(ci)=||mEci ||X2 + λ||Rci ||pp generated by the IRLS is monotonically decreasing. This guarantees the convergence since the sequence is lower-bounded by the finite quantity C=minc||mEc ||X2 + λ||Rc ||pp.

The IRLS algorithm can be simply adapted in order to solve the minimization with mixed-norm regularization terms. A particular case is the total variation penalty which corresponds to the ℓ1-norm of the pixel-wise ℓ2-norm of the spatial gradient [52]. The IRLS algorithm for TV regularization was first proposed by Wohlberg and Rodríguez [53]. It is described in Algorithm 8.

Algorithm 8   IRLS solving c=arg minc||mEc ||X2 + λ||c ||TV.
Inputs:
A=EHXE, a=EHXm, λ, and c0
Define the finite difference matrices
Rd along every spatial dimension d
Repeat until desired tolerance is reached: Return: ci

Duality-based algorithms proved to be an efficient alternative to achieve TV regularization [54,55].

3.4.4  Iterative Shrinkage/Thresholding Algorithm

The Iterative Shrinkage/Thresholding Algorithm (ISTA)[18,17,19], also known as thresholded Landweber (TL), aims at minimizing the functional

 C(w)=⎪⎪
⎪⎪
mMw ⎪⎪
⎪⎪
X2 + λ⎪⎪
⎪⎪
w ⎪⎪
⎪⎪
1.     (31)

Here, we use the notation w because ISTA is often applied on wavelet coefficients.

An important observation to understand ISTA is to see that the nonlinear shrinkage operation, sometimes called soft-thresholding, solves a minimization problem[56], with

     
Tλ(u=
 
|u|−min(λ /2,|u|)
· sgn(u) 
    (32)
=
 arg 
 
min
w∈ ℂ
 |uw|2 + λ |w|.
 

Figure 3.2: The shrinkage function Tλ.

By separability of norms, this applies component-wise to vectors of ℂN

Tλ(u) = arg 
 
min
w
 ⎪⎪
⎪⎪
uw ⎪⎪
⎪⎪
22 + λ ⎪⎪
⎪⎪
w ⎪⎪
⎪⎪
1.

This means that the ℓ1-regularized denoising problem (i.e., when M and X are identity matrices) is precisely solved by a shrinkage operation.

The ISTA generates a sequence of estimates wi that converges to the minimizer w of (31) when it is unique. The idea is to define at each step a new functional C′(w,wi) whose minimizer wi+1 will be the next estimate

wi+1 = arg 
 
min
w
  C′(w,wi).     (33)

Two constraints must be considered for the definition of C.

  1. It is sufficient for the convergence of the algorithm that C′(w,wi) is an upper bound of C(w) with equality at w=wi; this guarantees that the sequence { C(wi)} is monotonically decreasing.
  2. The inner minimization (33) should be performed by a simple shrinkage operation to ensure the rapidity and accuracy of the algorithm.

In accordance with Constraint 1, C can take the generic quadratically augmented form

 C′(w,wi) =  C(w) + ⎪⎪
⎪⎪
wwi ⎪⎪
⎪⎪
ΛMHXM2,     (34)

with the constraint that (ΛMHXM) is positive definite, where the weighting matrix Λ plays the role of a tuning parameter.

Then, ISTA corresponds to the trivial choice Λ=L/2I, with the value of L chosen to be greater or equal to the Lipschitz constant of the gradient of ||Mw ||X2, so that L≥ 2λmax( MHXM).

Let us define a = MHXm, A = MHXM, and

zi=wi+2(aAwi)/L.     (35)

Then, using standard linear algebra, we can write

wi+1=
arg 
 
min
w
⎪⎪
⎪⎪
wzi ⎪⎪
⎪⎪
22+(2λ/L)⎪⎪
⎪⎪
w ⎪⎪
⎪⎪
ℓ1
 =
T2λ/L
zi
.

This shows that Constraint 2 is automatically satisfied.

Note that both the intermediate variable zn in (35) and the threshold values will vary depending on L.

Algorithm 9   ISTA solving w=arg minw||mMw ||X2 + λ||w ||1.
Inputs:
A=MHXM, a=MHXm, w0, and L Repeat until desired tolerance is reached: Return: wi

Beck and Teboulle[20, Thm. 3.1] showed that this algorithm decreases the cost function in direct proportion to the number of iterations i.

Proposition 1   Let { wi} be the sequence generated by Algorithm 9 with L≥ 2λmax( A). Then, for any i>i0∈ ℕ,
 C(wi)− C(w)≤ 
L
2(ii0)
⎪⎪
⎪⎪
wi0w ⎪⎪
⎪⎪
22.     (36)

Selecting L as small as possible will clearly favor the speed of convergence. It also raises the importance of a “warm” starting point.

Among the variants of ISTA, FISTA, proposed by Beck and Teboulle[20], ensures state-of-the-art convergence properties while preserving a comparable computational cost. Thanks to a controlled over-relaxation at each step, FISTA quadratically decreases the cost function, with

 C(wi)− C(w) ≤ 
2L
(i+1)2
⎪⎪
⎪⎪
w0w ⎪⎪
⎪⎪
22.     (37)

More details on FISTA, as a particular case of FWISTA with the trivial choice Λ=L/2I, can be found in Section 5.2.3.

An implementation of FISTA is given Algorithm (10).

Algorithm 10   FISTA solving w=arg minw||mMw ||X2 + λ||w ||1
Inputs
A=MHXM, a=MHXm, w0, and L Initialization: i=0, v0 = w0, t0=1
Repeat until desired tolerance is reached:
Return: wi

1
Backprojection is a term used in tomography that refers to the multiplication by the transpose of the encoding matrix.
2
That is to say E[ e] = 0.
3
Let us mention that Υ being a covariance matrix, it is necessarily Hermitian symmetric. Its pseudoinverse is also Hermitian symmetric and admits the same eigenvectors.
4
Take a single channel MRI problem with M=256× 256 unknowns and N=256× 256 measurements. The corresponding encoding matrix E is 65 536× 65 536. With double precision floats that are required for accurate calculations, the storage of this complex-valued matrix would require 64 GiB of RAM memory which is a too large value for current personal computers. Performing direct matrix-vector multiplications involves an overwhelming amount of scalar multiplications and additions.

Previous Up Next