# Rank-One and Transformed Sparse Decomposition for Dynamic Cardiac MRI

Authors: Xianchao Xiu; Lingchen Kong
Journal: BioMed Research International
 (2015)

---

## Abstract

It is challenging and inspiring for us to achieve high spatiotemporal resolutions in dynamic cardiac magnetic resonance imaging (MRI). In this paper, we introduce two novel models and algorithms to reconstruct dynamic cardiac MRI data from under-sampledk
-
t space data. In contrast to classical low-rank and sparse model, we use rank-one and transformed sparse model to exploit the correlations in the dataset. In addition, we propose projected alternative direction method (PADM) and alternative hard thresholding method (AHTM) to solve our proposed models. Numerical experiments of cardiac perfusion and cardiac cine MRI data demonstrate improvement in performance.


---

## Body

## 1. Introduction

Dynamic magnetic resonance imaging (MRI) reconstructs a temporal series of images to resolve the motion or the variation of the imaged object. It is often used in cardiac, perfusion, functional, and gastrointestinal imaging. The dynamic contrast variations in the myocardium are typically imaged 40–60 s in cardiac perfusion imaging; however, most patients cannot maintain a breath-hold for such long durations. The data in MRI, samples ink
-
t space of the spatial Fourier transform of the object, are acquired sequentially in time. Therefore, achieving high spatiotemporal resolutions is challenging and inspiring in dynamic MRI due to the hardware limitations and the risk of peripheral nerve stimulation.Compressed sensing (CS) [1, 2] has been successfully applied to accelerate MRI, for example, [3–7] and references therein. Lustig et al. [3] systematically exploited the sparsity which is implicit in MRI and developed a framework for using CS in MRI. Liang [4] utilized low-rank matrix completion to dynamic MRI by considering each temporal frame as a column of a space-time matrix. Lingala et al. [5] and Zhao et al. [6] combined CS and low-rank matrix completion which finds a solution that is both low-rank and sparse. Gao et al. [8] proposed a different method which decomposes dynamic MRI matrix as a superposition of a low-rank component and a sparse component. It is worth noting that the low-rank component can model the temporally correlated background and the sparse component can model the dynamic information that lies on the top of the background. Otazo et al. [7] presented low-rank and sparse reconstructions for dynamic MRI using joint multicoil reconstruction for Cartesian and non-Cartesian k-space sampling.Most of the above work is concerned about convex relaxation; then it comes to an interesting question whether we can deal with dynamic MRI usingl
0-minimization. Many authors have made great effort; see [9–13], for example. Blumensath and Davies [9] were the first to propose the iterative hard thresholding (IHT) to solve a type of l
0-regularized problems and showed that the IHT method converges to a local minimizer. Beck and Eldar [10] introduced and analyzed three kinds of necessary optimality conditions for l
0-constrained problems. Bahmani et al. [14] studied l
0-constrained optimization in cases where nonlinear models are involved or the cost function is not quadratic. Recently, Lu and Zhang [11] presented a penalty decomposition method for solving a more general class of l
0-constrained and l
0-regularized minimization problems. Lu [12] considered l
0-regularized convex cone problems and then gave an IHT method and its convergence. More recently, Pan et al. [13] provided the first- and second-order necessary and sufficient optimality conditions for l
0-constrained optimization.Motivated by the recent advances in the use of CS in MRI andl
0-minimization, we propose rank-one and transformed sparse models and algorithms to significantly accelerate dynamic MRI. In contrast to classical low-rank and transformed sparse decomposition, this paper has the following merits:(i)
We introduce the use of rank-one matrix in dynamic cardiac MRI data which avoids singular value decomposition (SVD).
(ii)
We assume the dynamic component has a sparse representation in some known sparsifying transformation, not itself.
(iii)
We develop newl
0-constrained model and l
0-regularized model, rather than their convex relaxation.
(iv)
We design two ADM-based algorithms, projected alternative direction method (PADM) and alternative hard thresholding method (AHTM). Numerical experiments of dynamic cardiac MRI illustrate the efficiency of our proposed algorithms.The remainder of this paper is organised as follows. Section2 reviews some existing related work and proposes our new models. Section 3 describes PADM and AHTM to solve our proposed models. Section 4 reports experimental results. Section 5 concludes this paper with some final remarks.


## 2. Model Analysis

We begin this section by introducing the low-rank and transformed sparse decomposition in [7]; that is,(1)
m
i
n
L
∗
+
β
T
S
1
s
.
t
.
E
L
+
S
=
d
,where ·
∗ is defined as the sum of all singular values and ·
1 is defined as the sum of absolute values of all entries. β
∈
R is a real positive weighting parameter to balance the weights of rank and sparsity, and d is the undersampled k
-
t data. L
∈
C
m
×
n is a low-rank matrix, which can model the temporally correlated background. S
∈
C
m
×
n is a sparse matrix, which can model the dynamic information that lies on top of the background. Generally, the dynamic component S has a sparse representation in some known sparsifying transformation T (e.g., temporal frequency domain), hence the idea of minimizing T
S and not S itself. For acquisition with multiple receiver coils, E is given by the frame-by-frame multicoil encoding operator, which performs a multiplication by coil sensitivities followed by a Fourier transform according to the sampling pattern.Compared with a low-rank or sparse constraint, (1) can significantly increase compressibility of dynamic MRI data and thus enable high acceleration factors. However, (1) is the convex relaxation problem. The original optimization problem associated with (1) is formulated as follows:(2)
m
i
n
rank
⁡
L
+
β
T
S
0
s
.
t
.
E
L
+
S
=
d
,and its low-rank and sparse constrained optimization problem is as follows:(3)
m
i
n
E
L
+
S
-
d
2
s
.
t
.
rank
⁡
L
≤
r
T
S
0
≤
s
,where r denotes the rank of L, s denotes the sparsity of T
S, and ·
0 is l
0-norm, which counts the number of nonzero entries.Note that ifE and T reduce to identity matrix, (3) becomes a low-rank and sparse constrained optimization problem:(4)
m
i
n
L
+
S
-
D
2
s
.
t
.
rank
⁡
L
≤
r
S
0
≤
s
,where D
=
E
∗
d is given data. Zhou and Tao [15] developed GoDec method to estimate the low-rank part L and the sparse part S in (4). GoDec alternatively updated L by SVD and S by hard thresholding. To overcome the difficulty of SVD, they proposed bilateral random projections. Even so, it still costs much time.We find sometimes that the background is exactly a rank-one matrix, instead of a general low-rank matrix, so Li et al. [16] proposed to replace low-rank matrix with rank-one matrix which avoids any SVD completely. Xiu and Kong [17] extended it to the case of tensor decomposition and showed that it performs well in surveillance video. Specifically speaking, we set L
=
u
1
T, where u
∈
C
m and 1 denotes the vector in R
n whose entries are all 1, which leads to the following rank-one and transformed sparse matrix decomposition problem:(5)
m
i
n
E
u
1
T
+
S
-
d
2
s
.
t
.
T
S
0
≤
s
,and its unconstrained version(6)
min
⁡
E
u
1
T
+
S
-
d
2
+
λ
T
S
0
,where λ
>
0 is a regularized parameter.


## 3. Algorithm and Convergence

In this section, we present the projected alternative direction method (PADM) for (5) and alternative hard thresholding method (AHTM) for (6). Then we discuss their convergent properties.


### 3.1. Projected Alternative Direction Method for (5)

Let us look at the optimization problem (5). Inspired by ADM-based methods, it is easy to divide (5) into three subproblems:(7)
u
k
=
argmin
u
⁡
u
1
T
+
S
k
-
1
-
M
k
-
1
2
,
S
k
=
argmin
S
⁡
u
k
1
T
+
S
-
M
k
-
1
2
:
T
S
0
≤
s
,
M
k
=
u
k
1
T
+
S
k
-
E
∗
E
u
k
1
T
+
S
k
-
d
,where M
k is obtained by enforcing data consistency, where the aliasing artifacts corresponding to the residual in k
-
t space data E
∗
(
E
(
u
k
1
T
+
S
k
)
-
d
) are subtracted from u
k
1
T
+
S
k. Here, we denote these subproblems in (7) as u-subproblem, S-subproblem, and M-subproblem, respectively.The firstu-subproblem has the solution(8)
u
k
=
m
e
a
n
⁡
M
k
-
1
-
S
k
-
1
,where, for any Z
∈
C
m
×
n, mean(
·
)
∈
C
m is defined as (9)
m
e
a
n
Z
i
=
1
n
∑
j
=
1
n
Z
i
j
,
i
=
1,2
,
…
,
m
.Before establishing the solution of the second S-subproblem, we first give a proposition about sparse projection.Proposition 1.
LetΩ
=
S
:
T
S
0
≤
s. Then P
Ω
(
S
)
=
T
∗
(
T
S
)
max
⁡
(
s
), where P
Ω
(
S
) denotes all but s largest absolute value components of S.Proof.
LettingΩ
=
S
:
T
S
0
≤
s, we have(10)
P
Ω
S
=
argmin
TS
0
≤
s
⁡
X
-
S
2
=
argmin
TS
0
≤
s
⁡
T
X
-
T
S
2
=
argmin
Y
0
≤
s
⁡
T
X
-
Y
2
=
P
B
T
S
,where B
=
Y
:
Y
0
≤
s
=
Y
max
⁡
(
s
) and the third equation is satisfied by introducing Y
=
T
S. Then, it holds that (11)
P
Ω
S
=
T
∗
T
P
Ω
S
=
T
∗
P
B
T
S
=
T
∗
T
S
max
⁡
s
.The proof is completed.From Proposition1, we derive the solution of S-subproblem:(12)
S
k
=
P
Ω
M
k
-
1
-
u
k
1
T
.Hence the scheme is summarized as the following projected alternative direction method (PADM).In Algorithm1, M-subproblem is obtained by enforcing data consistency, so the convergence properties can be analyzed by considering the iterations of u-subproblem and S-subproblem. Now we will present the convergence result of PADM as follows.Algorithm 1:Projected alternative direction method for (5).
Input:
d: undersampled k
-
t space data
E: space-time encoding operator
T: sparsifying transform
Initialize: M
0
=
E
∗
d, S
0
=
0, s
While
u
k
=
m
e
a
n
(
M
k
-
1
-
S
k
-
1
)
S
k
=
P
Ω
(
M
k
-
1
-
u
k
1
T
)
M
k
=
u
k
1
T
+
S
k
-
E
∗
(
E
(
u
k
1
T
+
S
k
)
-
d
)
End
Output: L, STheorem 2.
Let{
(
u
k
,
S
k
)
} be the sequence generated by the above PADM, and (
u
∗
,
S
∗
) is an accumulation point of {
(
u
k
,
S
k
)
}. Then (
u
∗
,
S
∗
) is a local minimizer of (5).Proof.
LetF
(
u
,
S
)
=
E
u
1
T
+
S
-
d
2, for all S
∈
Ω. Indeed, one can observe that(13)
F
u
k
+
1
,
S
k
≤
F
u
,
S
k
,
F
u
k
,
S
k
+
1
≤
F
u
k
,
S
,
∀
S
∈
Ω
.It follows that (14)
F
u
k
+
1
,
S
k
+
1
≤
F
u
k
+
1
,
S
k
≤
F
u
k
,
S
k
.Hence, the sequence {
F
(
u
k
,
S
k
)
} is nonincreasing. Since (
u
∗
,
S
∗
) is an accumulation point of {
(
u
k
,
S
k
)
}, there exists a subsequence K such that l
i
m
k
∈
K
→
∞
(
u
k
,
S
k
)
=
(
u
∗
,
S
∗
). We then observe that {
F
(
u
k
,
S
k
)
}
k
∈
K is bounded, which together with the monotonicity of {
F
(
u
k
,
S
k
)
} implies that {
F
(
u
k
,
S
k
)
} is bounded below and hence l
i
m
k
→
∞
(
u
k
,
S
k
) exists. This observation and the continuity of {
F
(
u
k
,
S
k
)
} yield(15)
lim
k
→
∞
⁡
F
u
k
+
1
,
S
k
+
1
=
lim
k
→
∞
⁡
F
u
k
,
S
k
=
lim
k
∈
K
→
∞
⁡
F
u
k
,
S
k
=
F
u
∗
,
S
∗
.Using these relations and taking limits on both sides of (13) as k
∈
K
→
∞, we have(16)
F
u
∗
,
S
∗
≤
F
u
,
S
∗
,
F
u
∗
,
S
∗
≤
F
u
∗
,
S
,
∀
S
∈
Ω
.In addition, from the definition of Ω, we know that T
S
k
0
≤
s, which immediately implies T
S
∗
0
≤
s. Thus, (
u
∗
,
S
∗
) is a saddle point of (5). Let Δ
S be a variant such that S
∗
+
Δ
S
∈
Ω and let Δ
u be a vector in R
m. It then follows from (16) and the first-order optimality condition that [
∇
u
F
(
u
∗
,
S
∗
)
]
T
Δ
u
≥
0 and [
∇
S
F
(
u
∗
,
S
∗
)
]
T
Δ
S
≥
0. Henceforth, we obtain that(17)
F
u
∗
+
Δ
u
,
S
∗
+
Δ
S
≥
F
u
∗
,
S
∗
+
∇
u
F
u
∗
,
S
∗
T
Δ
u
+
∇
S
F
u
∗
,
S
∗
T
Δ
S
≥
F
u
∗
,
S
∗
,which implies that (
u
∗
,
S
∗
) is a local minimizer of (5).


### 3.2. Alternative Hard Thresholding Method for (6)

We can also divide the optimization problem (6) into the following u-subproblem, S-subproblem, and M-subproblem:(18)
u
k
=
argmin
u
⁡
u
1
T
+
S
k
-
1
-
M
k
-
1
2
,
S
k
=
argmin
S
⁡
u
k
1
T
+
S
-
M
k
-
1
2
+
λ
T
S
0
,
M
k
=
u
k
1
T
+
S
k
-
E
∗
E
u
k
1
T
+
S
k
-
d
.Compared with the above subsection, we find that the u-subproblem and M-subproblem are the same; hereafter we mainly explain how to solve the S-subproblem. By applying the separability of the objective function and the technique of operator splitting, the S-subproblem can be converted into m
n corresponding single variable minimization problems. Therefore, the following proposition will give the solution of the corresponding single variable minimization problem.Proposition 3.
Suppose thatx
∗ is a solution of problem (19)
min
⁡
φ
t
λ
x
=
x
2
-
2
t
x
+
λ
x
0
,where x
∈
R is variable and t
∈
R is a parameter. Then, (20)
x
∗
=
H
λ
0.5
t
=
t
,
i
f
t
>
λ
0.5
,
0
,
i
f
t
<
λ
0.5
,
t
o
r
0
,
o
t
h
e
r
w
i
s
e
.Proof.
Note thatφ
t
λ
(
0
)
=
0 and when x
≠
0, φ
t
λ
(
x
)
=
(
x
-
t
)
2
+
λ
-
t
2.(i)
Ift
>
λ
0.5, φ
t
λ
(
t
)
<
0, then x
∗
=
t.
(ii)
Ift
<
λ
0.5, φ
t
λ
(
t
)
>
0, then x
∗
=
0.
(iii)
Ift
=
λ
0.5, φ
t
λ
(
t
)
=
0, then x
∗
=
t or 0.Thus, we immediately obtain the proposition.It is straightforward to establish the following result.Proposition 4.
For the generall
0-regularized problem, (21)
min
⁡
A
X
-
B
2
+
λ
X
0
,where A
∈
R
p
×
m is some transform matrix and B is a given matrix in R
p
×
n. Then, the solution X
k is described as follows:(22)
X
k
=
H
λ
0.5
X
k
-
1
-
A
T
A
X
k
-
1
-
B
,where H
λ
0.5 is the hard thresholding operator defined as (23)
H
λ
0.5
X
i
j
=
X
i
j
,
i
f
X
i
j
>
λ
0.5
,
0
,
i
f
X
i
j
<
λ
0.5
,
X
i
j
o
r
0
,
o
t
h
e
r
w
i
s
e
.Proof.
The conclusion holds directly from Proposition3 by separating the objective function.From Proposition4, the solution of S-subproblem is characterized as(24)
S
k
=
T
∗
H
λ
0.5
T
M
k
-
1
-
u
k
1
T
.Based on the above argument, we have derived the following alternative hard thresholding method (AHTM) for (6).In the end of this subsection, we will establish the convergence result of Algorithm2.Algorithm 2:Alternative hard thresholding method for (6).
Input:
d: undersampled k
-
t space data
E: space-time encoding operator
T: sparsifying transform
Initialize: M
0
=
E
∗
d, S
0
=
0, λ
While
u
k
=
m
e
a
n
(
M
k
-
1
-
S
k
-
1
)
S
k
=
T
∗
(
H
λ
0.5
(
T
(
M
k
-
1
-
u
k
1
T
)
)
)
M
k
=
u
k
1
T
+
S
k
-
E
∗
(
E
(
u
k
1
T
+
S
k
)
-
d
)
End
Output: L, STheorem 5.
Let{
(
u
k
,
S
k
)
} be the sequence generated by the above AHTM, and (
u
∗
,
S
∗
) is an accumulation point of {
(
u
k
,
S
k
)
}. Then (
u
∗
,
S
∗
) is a local minimizer of (6).Proof.
Letf
(
u
,
S
)
=
E
u
1
T
+
S
-
d
2; then(25)
F
u
,
S
=
f
u
,
S
+
λ
T
S
0
=
E
u
1
T
+
S
-
d
2
+
λ
T
S
0
.Like the proof of Theorem 2, we derive that (
u
∗
,
S
∗
) is a saddle point of (6). We choose Δ
u to be a vector in R
m, and Δ
S satisfies T
S
∗
+
Δ
S
0
=
T
S
∗
0. Based on the convexity of f
(
u
,
S
), we have [
∇
u
f
(
u
∗
,
S
∗
)
]
T
Δ
u
≥
0 and [
∇
S
f
(
u
∗
,
S
∗
)
]
T
Δ
S
≥
0. Using this relation, we obtain that(26)
F
u
∗
+
Δ
u
,
S
∗
+
Δ
S
≥
f
u
∗
,
S
∗
+
∇
u
f
u
∗
,
S
∗
T
Δ
u
+
∇
S
f
u
∗
,
S
∗
T
Δ
S
+
λ
T
S
∗
+
Δ
S
0
≥
f
u
∗
,
S
∗
+
λ
T
S
∗
0
=
F
u
∗
,
S
∗
.Hence (
u
∗
,
S
∗
) is a local minimizer of (6).


## 3.1. Projected Alternative Direction Method for (5)

Let us look at the optimization problem (5). Inspired by ADM-based methods, it is easy to divide (5) into three subproblems:(7)
u
k
=
argmin
u
⁡
u
1
T
+
S
k
-
1
-
M
k
-
1
2
,
S
k
=
argmin
S
⁡
u
k
1
T
+
S
-
M
k
-
1
2
:
T
S
0
≤
s
,
M
k
=
u
k
1
T
+
S
k
-
E
∗
E
u
k
1
T
+
S
k
-
d
,where M
k is obtained by enforcing data consistency, where the aliasing artifacts corresponding to the residual in k
-
t space data E
∗
(
E
(
u
k
1
T
+
S
k
)
-
d
) are subtracted from u
k
1
T
+
S
k. Here, we denote these subproblems in (7) as u-subproblem, S-subproblem, and M-subproblem, respectively.The firstu-subproblem has the solution(8)
u
k
=
m
e
a
n
⁡
M
k
-
1
-
S
k
-
1
,where, for any Z
∈
C
m
×
n, mean(
·
)
∈
C
m is defined as (9)
m
e
a
n
Z
i
=
1
n
∑
j
=
1
n
Z
i
j
,
i
=
1,2
,
…
,
m
.Before establishing the solution of the second S-subproblem, we first give a proposition about sparse projection.Proposition 1.
LetΩ
=
S
:
T
S
0
≤
s. Then P
Ω
(
S
)
=
T
∗
(
T
S
)
max
⁡
(
s
), where P
Ω
(
S
) denotes all but s largest absolute value components of S.Proof.
LettingΩ
=
S
:
T
S
0
≤
s, we have(10)
P
Ω
S
=
argmin
TS
0
≤
s
⁡
X
-
S
2
=
argmin
TS
0
≤
s
⁡
T
X
-
T
S
2
=
argmin
Y
0
≤
s
⁡
T
X
-
Y
2
=
P
B
T
S
,where B
=
Y
:
Y
0
≤
s
=
Y
max
⁡
(
s
) and the third equation is satisfied by introducing Y
=
T
S. Then, it holds that (11)
P
Ω
S
=
T
∗
T
P
Ω
S
=
T
∗
P
B
T
S
=
T
∗
T
S
max
⁡
s
.The proof is completed.From Proposition1, we derive the solution of S-subproblem:(12)
S
k
=
P
Ω
M
k
-
1
-
u
k
1
T
.Hence the scheme is summarized as the following projected alternative direction method (PADM).In Algorithm1, M-subproblem is obtained by enforcing data consistency, so the convergence properties can be analyzed by considering the iterations of u-subproblem and S-subproblem. Now we will present the convergence result of PADM as follows.Algorithm 1:Projected alternative direction method for (5).
Input:
d: undersampled k
-
t space data
E: space-time encoding operator
T: sparsifying transform
Initialize: M
0
=
E
∗
d, S
0
=
0, s
While
u
k
=
m
e
a
n
(
M
k
-
1
-
S
k
-
1
)
S
k
=
P
Ω
(
M
k
-
1
-
u
k
1
T
)
M
k
=
u
k
1
T
+
S
k
-
E
∗
(
E
(
u
k
1
T
+
S
k
)
-
d
)
End
Output: L, STheorem 2.
Let{
(
u
k
,
S
k
)
} be the sequence generated by the above PADM, and (
u
∗
,
S
∗
) is an accumulation point of {
(
u
k
,
S
k
)
}. Then (
u
∗
,
S
∗
) is a local minimizer of (5).Proof.
LetF
(
u
,
S
)
=
E
u
1
T
+
S
-
d
2, for all S
∈
Ω. Indeed, one can observe that(13)
F
u
k
+
1
,
S
k
≤
F
u
,
S
k
,
F
u
k
,
S
k
+
1
≤
F
u
k
,
S
,
∀
S
∈
Ω
.It follows that (14)
F
u
k
+
1
,
S
k
+
1
≤
F
u
k
+
1
,
S
k
≤
F
u
k
,
S
k
.Hence, the sequence {
F
(
u
k
,
S
k
)
} is nonincreasing. Since (
u
∗
,
S
∗
) is an accumulation point of {
(
u
k
,
S
k
)
}, there exists a subsequence K such that l
i
m
k
∈
K
→
∞
(
u
k
,
S
k
)
=
(
u
∗
,
S
∗
). We then observe that {
F
(
u
k
,
S
k
)
}
k
∈
K is bounded, which together with the monotonicity of {
F
(
u
k
,
S
k
)
} implies that {
F
(
u
k
,
S
k
)
} is bounded below and hence l
i
m
k
→
∞
(
u
k
,
S
k
) exists. This observation and the continuity of {
F
(
u
k
,
S
k
)
} yield(15)
lim
k
→
∞
⁡
F
u
k
+
1
,
S
k
+
1
=
lim
k
→
∞
⁡
F
u
k
,
S
k
=
lim
k
∈
K
→
∞
⁡
F
u
k
,
S
k
=
F
u
∗
,
S
∗
.Using these relations and taking limits on both sides of (13) as k
∈
K
→
∞, we have(16)
F
u
∗
,
S
∗
≤
F
u
,
S
∗
,
F
u
∗
,
S
∗
≤
F
u
∗
,
S
,
∀
S
∈
Ω
.In addition, from the definition of Ω, we know that T
S
k
0
≤
s, which immediately implies T
S
∗
0
≤
s. Thus, (
u
∗
,
S
∗
) is a saddle point of (5). Let Δ
S be a variant such that S
∗
+
Δ
S
∈
Ω and let Δ
u be a vector in R
m. It then follows from (16) and the first-order optimality condition that [
∇
u
F
(
u
∗
,
S
∗
)
]
T
Δ
u
≥
0 and [
∇
S
F
(
u
∗
,
S
∗
)
]
T
Δ
S
≥
0. Henceforth, we obtain that(17)
F
u
∗
+
Δ
u
,
S
∗
+
Δ
S
≥
F
u
∗
,
S
∗
+
∇
u
F
u
∗
,
S
∗
T
Δ
u
+
∇
S
F
u
∗
,
S
∗
T
Δ
S
≥
F
u
∗
,
S
∗
,which implies that (
u
∗
,
S
∗
) is a local minimizer of (5).


## 3.2. Alternative Hard Thresholding Method for (6)

We can also divide the optimization problem (6) into the following u-subproblem, S-subproblem, and M-subproblem:(18)
u
k
=
argmin
u
⁡
u
1
T
+
S
k
-
1
-
M
k
-
1
2
,
S
k
=
argmin
S
⁡
u
k
1
T
+
S
-
M
k
-
1
2
+
λ
T
S
0
,
M
k
=
u
k
1
T
+
S
k
-
E
∗
E
u
k
1
T
+
S
k
-
d
.Compared with the above subsection, we find that the u-subproblem and M-subproblem are the same; hereafter we mainly explain how to solve the S-subproblem. By applying the separability of the objective function and the technique of operator splitting, the S-subproblem can be converted into m
n corresponding single variable minimization problems. Therefore, the following proposition will give the solution of the corresponding single variable minimization problem.Proposition 3.
Suppose thatx
∗ is a solution of problem (19)
min
⁡
φ
t
λ
x
=
x
2
-
2
t
x
+
λ
x
0
,where x
∈
R is variable and t
∈
R is a parameter. Then, (20)
x
∗
=
H
λ
0.5
t
=
t
,
i
f
t
>
λ
0.5
,
0
,
i
f
t
<
λ
0.5
,
t
o
r
0
,
o
t
h
e
r
w
i
s
e
.Proof.
Note thatφ
t
λ
(
0
)
=
0 and when x
≠
0, φ
t
λ
(
x
)
=
(
x
-
t
)
2
+
λ
-
t
2.(i)
Ift
>
λ
0.5, φ
t
λ
(
t
)
<
0, then x
∗
=
t.
(ii)
Ift
<
λ
0.5, φ
t
λ
(
t
)
>
0, then x
∗
=
0.
(iii)
Ift
=
λ
0.5, φ
t
λ
(
t
)
=
0, then x
∗
=
t or 0.Thus, we immediately obtain the proposition.It is straightforward to establish the following result.Proposition 4.
For the generall
0-regularized problem, (21)
min
⁡
A
X
-
B
2
+
λ
X
0
,where A
∈
R
p
×
m is some transform matrix and B is a given matrix in R
p
×
n. Then, the solution X
k is described as follows:(22)
X
k
=
H
λ
0.5
X
k
-
1
-
A
T
A
X
k
-
1
-
B
,where H
λ
0.5 is the hard thresholding operator defined as (23)
H
λ
0.5
X
i
j
=
X
i
j
,
i
f
X
i
j
>
λ
0.5
,
0
,
i
f
X
i
j
<
λ
0.5
,
X
i
j
o
r
0
,
o
t
h
e
r
w
i
s
e
.Proof.
The conclusion holds directly from Proposition3 by separating the objective function.From Proposition4, the solution of S-subproblem is characterized as(24)
S
k
=
T
∗
H
λ
0.5
T
M
k
-
1
-
u
k
1
T
.Based on the above argument, we have derived the following alternative hard thresholding method (AHTM) for (6).In the end of this subsection, we will establish the convergence result of Algorithm2.Algorithm 2:Alternative hard thresholding method for (6).
Input:
d: undersampled k
-
t space data
E: space-time encoding operator
T: sparsifying transform
Initialize: M
0
=
E
∗
d, S
0
=
0, λ
While
u
k
=
m
e
a
n
(
M
k
-
1
-
S
k
-
1
)
S
k
=
T
∗
(
H
λ
0.5
(
T
(
M
k
-
1
-
u
k
1
T
)
)
)
M
k
=
u
k
1
T
+
S
k
-
E
∗
(
E
(
u
k
1
T
+
S
k
)
-
d
)
End
Output: L, STheorem 5.
Let{
(
u
k
,
S
k
)
} be the sequence generated by the above AHTM, and (
u
∗
,
S
∗
) is an accumulation point of {
(
u
k
,
S
k
)
}. Then (
u
∗
,
S
∗
) is a local minimizer of (6).Proof.
Letf
(
u
,
S
)
=
E
u
1
T
+
S
-
d
2; then(25)
F
u
,
S
=
f
u
,
S
+
λ
T
S
0
=
E
u
1
T
+
S
-
d
2
+
λ
T
S
0
.Like the proof of Theorem 2, we derive that (
u
∗
,
S
∗
) is a saddle point of (6). We choose Δ
u to be a vector in R
m, and Δ
S satisfies T
S
∗
+
Δ
S
0
=
T
S
∗
0. Based on the convexity of f
(
u
,
S
), we have [
∇
u
f
(
u
∗
,
S
∗
)
]
T
Δ
u
≥
0 and [
∇
S
f
(
u
∗
,
S
∗
)
]
T
Δ
S
≥
0. Using this relation, we obtain that(26)
F
u
∗
+
Δ
u
,
S
∗
+
Δ
S
≥
f
u
∗
,
S
∗
+
∇
u
f
u
∗
,
S
∗
T
Δ
u
+
∇
S
f
u
∗
,
S
∗
T
Δ
S
+
λ
T
S
∗
+
Δ
S
0
≥
f
u
∗
,
S
∗
+
λ
T
S
∗
0
=
F
u
∗
,
S
∗
.Hence (
u
∗
,
S
∗
) is a local minimizer of (6).


## 4. Numerical Experiments

In this section, we conduct numerical experiments to compare the performance of IST [7], PADM, and AHTM for dynamic cardiac MRI data. We apply all these methods on MR images: cardiac perfusion and cardiac cine, which are available at the website http://cai2r.net/resources/software/ls-reconstruction-matlab-code. All experiments are implemented in MATLAB (MathWorks, Natick, MA) on a desktop computer with Intel Core I5 2.60 GHz CPU and 8 GB of RAM.We now address the initialization and the termination criteria for these methods. In particular, we set the stopping criterion asR
e
l
E
r
r
<
2.5
E
-
3, where(27)
R
e
l
E
r
r
≔
u
k
1
T
+
S
k
-
u
k
-
1
1
T
+
S
k
-
1
2
u
k
-
1
1
T
+
S
k
-
1
2
.We choose λ
=
0.01 in AHTM as suggested in [7]. However, we do not know the sparse parameter s beforehand. Hence, in this paper, we initialize s as one percent of pixels in each frame of T
S.


### 4.1. Cardiac Perfusion

We next conduct numerical experiments to test the performance of IST [7], PADM, and AHTM for dynamic cardiac perfusion MRI data. The computational results of the three methods are presented in Figures 1–3. In detail, we use L
+
S (the first row) to denote the original dynamic cardiac perfusion MRI frames, L (the second row) to denote the reconstructed background parts, and S (the third row) to denote the reconstructed moving parts. Comparing Figures 1–3, we can observe that PADM (see Figure 2) and AHTM (see Figure 3) are easier to recognize the pathological parts than IST (see Figure 1), which is due to the use of l
0 rather than its convex relaxation l
1. This capacity may be useful to identify lesions that are difficult to visualize in the original images.Figure 1
Performance of IST for dynamic cardiac perfusion.Figure 2
Performance of PADM for dynamic cardiac perfusion.Figure 3
Performance of AHTM for dynamic cardiac perfusion.In addition, Table1 reports the computational results for dynamic cardiac perfusion MRI data. Here, the number of iterations is given in the second column (i.e., Iteration). The third column (i.e., Time) lists the computing time in seconds. The relative error is reported in the last column (i.e., RelErr). It is evident to see that the “Time” and “Iteration” of PADM and AHTM are lower than these of IST. This is because we employ rank-one to avoid SVD.Table 1
Computational results for dynamic cardiac perfusion.
Algorithm
Iteration
Time
RelErr
IST
36
80.38
2.5
E
-
03
PADM
25
59.47
2.5
E
-
03
AHTM
31
63.25
2.4
E
-
03


### 4.2. Cardiac Cine

Similarly, Figures4–6 and Table 2 demonstrate the performance of IST [7], PADM, and AHTM for dynamic cardiac cine MRI data. From the third row, we find that Figures 5 and 6 can get clearer moving parts, because we can tune the sparsity of s and λ.Table 2
Computational results for dynamic cardiac cine.
Algorithm
Iteration
Time
RelErr
IST
26
140.13
2.4
E
-
03
PADM
18
81.49
2.5
E
-
03
AHTM
22
94.26
2.4
E
-
03Figure 4
Performance of IST for dynamic cardiac cine.Figure 5
Performance of PADM for dynamic cardiac cine.Figure 6
Performance of AHTM for dynamic cardiac cine.


## 4.1. Cardiac Perfusion

We next conduct numerical experiments to test the performance of IST [7], PADM, and AHTM for dynamic cardiac perfusion MRI data. The computational results of the three methods are presented in Figures 1–3. In detail, we use L
+
S (the first row) to denote the original dynamic cardiac perfusion MRI frames, L (the second row) to denote the reconstructed background parts, and S (the third row) to denote the reconstructed moving parts. Comparing Figures 1–3, we can observe that PADM (see Figure 2) and AHTM (see Figure 3) are easier to recognize the pathological parts than IST (see Figure 1), which is due to the use of l
0 rather than its convex relaxation l
1. This capacity may be useful to identify lesions that are difficult to visualize in the original images.Figure 1
Performance of IST for dynamic cardiac perfusion.Figure 2
Performance of PADM for dynamic cardiac perfusion.Figure 3
Performance of AHTM for dynamic cardiac perfusion.In addition, Table1 reports the computational results for dynamic cardiac perfusion MRI data. Here, the number of iterations is given in the second column (i.e., Iteration). The third column (i.e., Time) lists the computing time in seconds. The relative error is reported in the last column (i.e., RelErr). It is evident to see that the “Time” and “Iteration” of PADM and AHTM are lower than these of IST. This is because we employ rank-one to avoid SVD.Table 1
Computational results for dynamic cardiac perfusion.
Algorithm
Iteration
Time
RelErr
IST
36
80.38
2.5
E
-
03
PADM
25
59.47
2.5
E
-
03
AHTM
31
63.25
2.4
E
-
03


## 4.2. Cardiac Cine

Similarly, Figures4–6 and Table 2 demonstrate the performance of IST [7], PADM, and AHTM for dynamic cardiac cine MRI data. From the third row, we find that Figures 5 and 6 can get clearer moving parts, because we can tune the sparsity of s and λ.Table 2
Computational results for dynamic cardiac cine.
Algorithm
Iteration
Time
RelErr
IST
26
140.13
2.4
E
-
03
PADM
18
81.49
2.5
E
-
03
AHTM
22
94.26
2.4
E
-
03Figure 4
Performance of IST for dynamic cardiac cine.Figure 5
Performance of PADM for dynamic cardiac cine.Figure 6
Performance of AHTM for dynamic cardiac cine.


## 5. Conclusions

In this paper, we focused on the application of rank-one and transformed sparse decomposition for dynamic cardiac MRI. We established the projected alternative direction method and alternative hard thresholding method and showed their convergent results. Finally, numerical experiments demonstrated the efficiency and effectiveness of our proposed methods.



---

*Source: 169317-2015-07-12.xml*