# Dynamic MRI reconstruction using low-rank plus sparse decomposition with
smoothness regularization

Chee-Ming Ting
School of Information Technology
Monash University Malaysia
&Fuad Noman
School of Information Technology
Monash University Malaysia
&Raphaël C.-W. Phan
School of Information Technology
Monash University Malaysia
&Hernando Ombao
Statistics Program
King Abdullah University of Science and Technology

###### Abstract

The low-rank plus sparse (L+S) decomposition model has enabled better
reconstruction of dynamic magnetic resonance imaging (dMRI) with separation
into background (L) and dynamic (S) component. However, use of low-rank prior
alone may not fully explain the slow variations or smoothness of the
background part at the local scale. In this paper, we propose a smoothness-
regularized L+S (SR-L+S) model for dMRI reconstruction from highly
undersampled k-t-space data. We exploit joint low-rank and smooth priors on
the background component of dMRI to better capture both its global and local
temporal correlated structures. Extending the L+S formulation, the low-rank
property is encoded by the nuclear norm, while the smoothness by a general
$\ell_{p}$-norm penalty on the local differences of the columns of L. The
additional smoothness regularizer can promote piecewise local consistency
between neighboring frames. By smoothing out the noise and dynamic activities,
it allows accurate recovery of the background part, and subsequently more
robust dMRI reconstruction. Extensive experiments on multi-coil cardiac and
synthetic data shows that the SR-L+S model outperforms several state-of-the-
art methods in terms of recovery accuracy.

_Keywords_ Dynamic MRI $\cdot$ low-rank $\cdot$ sparsity $\cdot$ smoothness
$\cdot$ proximal gradient.

##  1 Introduction

Dynamic magnetic resonance imaging (dMRI) is a powerful imaging modality
widely used in clinical applications such as cardiac and perfusion imaging. It
can reveal the spatial structure and temporal evolution of organs of interest.
However, traditional dMRI is limited by slow imaging speed and low signal-to-
noise ratio due to the additional temporal dimension and increased chance of
imaging artifacts. The application of compressed sensing (CS) Lustig et al.
(2007) combined with parallel MRI technique such as SENSE Pruessmann et al.
(1999) has enabled reconstruction of high-quality images from highly
undersampled k-space data, significantly accelerating the dMRI process Otazo
et al. (2010). To achieve undersampled MRI reconstruction, CS-based methods
leverage on sparsity representation of images in original or transform domain,
e.g., wavelets and temporal Fourier Lustig et al. (2007, 2006); Otazo et al.
(2010).

In addition to the sparse prior of CS, low-rank models have also been explored
for dMRI to characterize the temporal information correlation or redundancy
across dynamic MR frames arising from slow changes of the same tissues.
Earlier work on low-rank and sparse (L&S) like kt-SLR Lingala et al. (2011)
assume image sequence as both low-rank and sparse. Alternatively, the low-rank
plus sparse (L+S) decomposition is more natural for dynamic imaging because it
separates the space-time data matrix ${\bf X}={\bf L}+{\bf S}$ into a low-rank
component (${\bf L}$) and a sparse component (${\bf S}$) that can represent
the temporally-correlated background and dynamic foreground activities,
respectively Trémoulhéac et al. (2014); Otazo et al. (2015). By exploiting the
inherent low-rankness of the slowly-changing background and stronger sparsity
of images after the background subtracted, L+S has achieved better
reconstruction quality than the CS-based methods. Various efficient algorithms
have been developed to solve low-rank models for dMRI reconstruction. These
include fast algorithms based on proximal gradient and variable splitting
methods Lin and Fessler (2018), and alternating gradient descent and
minimization for MRI (altGDminMRI) Babu et al. (2023). However, the low-rank
regularization used in the above-mentioned studies can only describe the
common low-dimensional subspace of the whole image sequence, but may not fully
explain the slow variation or smoothness at the local range of the dynamic
frames (e.g., smooth contrast changes in perfusion MRI).

In this paper, we propose a novel smoothness-regularized L+S (SR-L+S)
framework for dynamic MRI reconstruction. Building on the L+S model, we
leverage on both the low-rank and smoothness priors to simultaneously capture
the global and local temporal correlations in the background components of the
dMRI sequence. The reconstruction problem is formulated as a convex
optimization, where nuclear norm and $\ell_{1}$ norm are used to,
respectively, promote low-rankness and and sparsity on ${\bf L}$ and ${\bf
S}$. We further incorporate a general smoothness regularizer that uses a
$\ell_{p}$-norm penalty on the successive differences of the columns of the
low-rank matrix ${\bf L}$ to enforce local similarity between neighbouring
frames along the temporal dimension. The additional regularizer can provide
further smoothing out of the random noise and foreground activities, thus
enabling better recovery of the background component ${\bf L}$, and
subsequently more robust image reconstruction. We consider both $\ell_{2}$ and
$\ell_{1}$ smoothness regularization penalties, where the $\ell_{1}$ special
case is an alternative formulation of the total variation (TV).

The main contributions of this work are as follows: (1) To our best knowledge,
this is among the first work to incorporate temporal smoothness regularization
into L+S reconstruction of dMRI. Although previous studies Lingala et al.
(2011); Yao et al. (2018) have combined low-rank and TV-based smoothness
regularization for dMRI, they applied them on the reconstructed image ${\bf
X}$ itself and not within the L+S framework as in the SR-L+S model that can
facilitate simultaneous robust reconstruction and separation of dMRI. TV has
also been used in L+S-based dMRI but only on the sparse component Wang et al.
(2020). Moreover, these studies applied the TV regularization mainly in the
spatial domain. In contrast, our model exploits the joint low-rank and
smoothness priors for better modeling of the background component in the
temporal domain to allow its accurate recovery from undersampled data. (2) We
develop an efficient algorithm based on proximal gradient method to solve a
convex optimization problem for SR-L+S model which consists of three
regularization terms: nuclear norm for background part ${\bf L}$,
$\ell_{p}$-norm for local smoothness in ${\bf L}$, and $\ell_{1}$-based
sparsity for ${\bf S}$, with closed-form updating equation for each involved
variable. (3) Experiments on cardiac perfusion and PINCAT synthetic data show
significant improvement on recovery accuracy of the proposed model
quantitatively and qualitatively over state-of-the-art low-rank methods.

##  2 METHODS

###  2.1 L+S decomposition for dMRI reconstruction

Consider a sequence of $N_{t}$ dynamic MR images of dimension $N_{x}\times
N_{y}$. It is usually formulated as a space-time matrix (Casorati matrix)
${\bf X}=[{\bf x}_{1},\ldots,{\bf x}_{N_{t}}]\in\mathbb{C}^{N_{x}N_{y}\times
N_{t}}$ in which each column represents a vectorized temporal image or frame.
A typical linear imaging model for dMRI can be written as

${\bf y}={\bf E}({\bf X})+{\bf n}$ (1)

where ${\bf y}\in\mathbb{C}^{P}$ is the undersampled ($k$-$t$)-space data
measured during the acquisition, ${\bf E}:\mathbb{C}^{N_{x}N_{y}\times
N_{t}}\rightarrow\mathbb{C}^{P}$ is the encoding or acquisition operator, and
${\bf n}$ is the noise vector. For parallel imaging with multiple receiver
coils, the dimension of the k-space data is $P=N_{s}N_{c}$ with $N_{s}$ the
total number of samples received from each receiver coil (across all frames)
and $N_{c}$ is the number of coils, and ${\bf E}={\bf F}_{u}{\bf S}_{c}$ where
${\bf F}_{u}$ is Fourier transform with undersampling and ${\bf S}_{c}$
denotes the coil sensitivities. The problem of image reconstruction is to
recover the clean image ${\bf X}$ from the undersampled data ${\bf y}$, which
is ill-posed under scenario $P\ll N_{x}N_{y}\times N_{t}$ and requires
regularization.

We consider the L+S framework Otazo et al. (2015); Lin and Fessler (2018) for
reconstruction and separation of dMRI data, modeled as a superposition of a
low-rank component ${\bf L}$ and a sparse component ${\bf S}$, ${\bf X}={\bf
L}+{\bf S}$. The low-rank matrix ${\bf L}$ represents the background component
of images, which is assumed to change slowly over time or exhibit high
correlation among frames. The ${\bf S}$ corresponds to the dynamic component
which is rapidly changing. It captures the innovation in each frame in ${\bf
X}$ with its background suppressed, which is already sparser that the original
image, and has a much sparser representation under a proper sparsifying
transform. The image reconstruction via L+S decomposition can be formulated as
a convex minimization problem

$\min_{{\bf L},{\bf S}}\frac{1}{2}\left\|{\bf E}({\bf L}+{\bf S})-{\bf
y}\right\|_{2}^{2}+\lambda_{L}\left\|{\bf
L}\right\|_{*}+\lambda_{S}\left\|\bf{T}({\bf S})\right\|_{1}$ (2)

where $\left\|{\bf L}\right\|_{*}$ is the nuclear norm or sum of singular
values of ${\bf L}$, and $\left\|\bf{T}({\bf S})\right\|_{1}$ is the
$\ell_{1}$ norm of transformed ${\bf S}$ under a sparsifying transform ${\bf
T}$. As in Otazo et al. (2015); Lin and Fessler (2018), we consider the
Fourier transform operator along the temporal dimension for ${\bf T}$. The
data fidelity term is captured by the $\ell_{2}$-norm term, the low-rankness
of ${\bf L}$ is induced by the nuclear norm, and the sparsity of the
transformed ${\bf S}$ by the $\ell_{1}$ penalty. The trade-off between these
terms is controlled by the regularization parameters $\lambda_{L}$ and
$\lambda_{S}$.

###  2.2 The proposed model: L+S with smoothness

The slowly-varying or possibly constant background components across frames
may not be fully explained only by the low-rank prior in (2). We exploit two
types of prior knowledge jointly on the background component ${\bf L}$ for
more robust image reconstruction: (1) Low-rankness which considers that the
entire temporal sequence in ${\bf L}=[{\bf l}_{1},\ldots,{\bf l}_{N_{t}}]$
resides in a common low-dimensional subspace, revealing highly-correlated
information at a global scale, and (2) Smoothness which considers that
neighbouring frames along the temporal dimension tend to vary smoothly,
reflecting local information correlation or similarity at a relatively local
scale.

By jointly applying the local and global correlated priors, we propose a novel
smoothness-regularized L+S (SR-L+S) model for dMRI reconstruction, which can
be formulated as

$\min_{{\bf L},{\bf S}}\frac{1}{2}\left\|{\bf E}({\bf L}+{\bf S})-{\bf
y}\right\|_{2}^{2}+\lambda_{L}\left\|{\bf
L}\right\|_{*}+\lambda_{S}\left\|\bf{T}({\bf
S})\right\|_{1}+\lambda_{D}\mathcal{R}({\bf L})$ (3)

where

$\mathcal{R}({\bf L})=\sum_{t=2}^{N_{t}}(\left\|{\bf l}_{t}-{\bf
l}_{t-1}\right\|_{p})^{p}$ (4)

is smoothness regularizer controlled by parameter $\lambda_{D}$. Minimizing
the $\ell_{p}$ ($p\leq 2$) norm on the differences of neighboring frames in
(4) can further enhance the temporal smoothness of the background component
${\bf L}$. We consider special cases $p=2$ and $p=1$. Let ${\bf D}$ be first-
order finite difference matrix

${\bf D}=\left[\begin{array}[]{ccccc}-1&1&0&\ldots&0\\\ 0&-1&1&\ldots&0\\\
\vdots&\vdots&\ddots&\ddots&\vdots\\\ 0&0&\ldots&-1&1\\\
\end{array}\right]\in\mathbb{R}^{(N_{t}-1)\times N_{t}}$

When $p=2$, (4) simplifies to $\mathcal{R}({\bf L})=\left\|{\bf L}{\bf
D}^{T}\right\|_{F}^{2}=\text{tr}({\bf L}{\bf D}^{T}{\bf D}{\bf L}^{T})$ where
$\left\|\cdot\right\|_{F}$ denotes the Frobenius norm. When $p=1$,
$\mathcal{R}({\bf L})=\left\|{\bf L}{\bf D}^{T}\right\|_{1}$ which is an
alternative formulation of the popular temporal total variation (TV).

###  2.3 Optimization

We develop a proximal gradient method (PGM) to solve the optimization problem
in (3) for the proposed RS-L+S decomposition. Here, we present the algorithm
for the $\ell_{1}$ smoothness prior, which can be applied similarly for the
$\ell_{2}$ case. By introducing an auxiliary variable $\boldsymbol{\Phi}={\bf
L}{\bf D}^{T}\in\mathbb{C}^{N_{x}N_{y}\times N_{t}-1}$, the problem is
equivalent to

$\displaystyle\min_{{\bf L},{\bf S},\boldsymbol{\Phi}}$
$\displaystyle\frac{1}{2}\left\|{\bf E}({\bf L}+{\bf S})-{\bf
y}\right\|_{2}^{2}+\lambda_{L}\left\|{\bf
L}\right\|_{*}+\lambda_{S}\left\|\bf{T}({\bf
S})\right\|_{1}+\lambda_{D}\left\|\boldsymbol{\Phi}\right\|_{1}$
$\displaystyle s.t.\ \ \boldsymbol{\Phi}={\bf L}{\bf D}^{T}$ (5)

We further relax the equality constraint by transforming (5) to

$\displaystyle\min_{{\bf L},{\bf S},\boldsymbol{\Phi}}$
$\displaystyle\frac{\mu}{2}\left\|{\bf E}({\bf L}+{\bf S})-{\bf
y}\right\|_{2}^{2}+\frac{1}{2}\left\|{\bf L}{\bf
D}^{T}-\boldsymbol{\Phi}\right\|_{F}^{2}$
$\displaystyle+\mu(\lambda_{L}\left\|{\bf
L}\right\|_{*}+\lambda_{S}\left\|\bf{T}({\bf
S})\right\|_{1}+\lambda_{D}\left\|\boldsymbol{\Phi}\right\|_{1})$ (6)

where $\mu>0$ is a relaxation parameter.

By defining $\boldsymbol{\Theta}=({\bf L};{\bf S};\boldsymbol{\Phi})$, the
problem (6) is a special case of a general convex problem of the form:

$\min_{\boldsymbol{\Theta}}F(\boldsymbol{\Theta})=g(\boldsymbol{\Theta})+f(\boldsymbol{\Theta})$
(7)

with

$\displaystyle g(\boldsymbol{\Theta})$
$\displaystyle=\mu(\lambda_{L}\left\|{\bf
L}\right\|_{*}+\lambda_{S}\left\|\bf{T}({\bf
S})\right\|_{1}+\lambda_{D}\left\|\boldsymbol{\Phi}\right\|_{1}),$
$\displaystyle f(\boldsymbol{\Theta})$ $\displaystyle=\frac{\mu}{2}\left\|{\bf
E}({\bf L}+{\bf S})-{\bf y}\right\|_{2}^{2}+\frac{1}{2}\left\|{\bf L}{\bf
D}^{T}-\boldsymbol{\Phi}\right\|_{F}^{2}$

where $g$ is convex but not necessarily smooth, $f$ is smooth, convex and
differentiable, and its gradient $\nabla f(\boldsymbol{\Theta})$ is Lipschitz
continuous with constant $L_{f}>0$. The PGM forms a local approximation to the
smooth term $f(\boldsymbol{\Theta})$ at a chosen point
${\boldsymbol{\Theta}}^{k}=({\bf L}^{k};{\bf S}^{k};{\boldsymbol{\Phi}}^{k})$,
and updates $\boldsymbol{\Theta}$ iteratively using proximal operator

$\displaystyle\boldsymbol{\Theta}_{k+1}$
$\displaystyle=\operatorname*{arg\,min}_{\boldsymbol{\Theta}}\
g(\boldsymbol{\Theta})+\frac{L_{f}}{2}\|\boldsymbol{\Theta}-{\bf
G}^{k}\|_{F}^{2}$ $\displaystyle\text{where}\ \ \ {\bf G}_{k}$
$\displaystyle={\boldsymbol{\Theta}}_{k}-\frac{1}{L_{f}}\nabla
f({\boldsymbol{\Theta}}_{k})$ (8)

The problem (8) can be separable into the following independent subproblems:

$\displaystyle{\bf L}_{k+1}$ $\displaystyle=\operatorname*{arg\,min}_{{\bf
X}}\mu\lambda_{L}\left\|{\bf L}\right\|_{*}+\frac{L_{f}}{2}\|{\bf L}-{\bf
G}^{L}_{k}\|_{F}^{2}$ (9) $\displaystyle{\bf S}_{k+1}$
$\displaystyle=\operatorname*{arg\,min}_{{\bf
S}}\mu\lambda_{S}\left\|\bf{T}({\bf S})\right\|_{1}+\frac{L_{f}}{2}\|{\bf
S}-{\bf G}^{S}_{k}\|_{F}^{2}$ (10) $\displaystyle\boldsymbol{\Phi}_{k+1}$
$\displaystyle=\operatorname*{arg\,min}_{\boldsymbol{\Phi}}\mu\lambda_{D}\left\|\boldsymbol{\Phi}\right\|_{1}+\frac{L_{f}}{2}\|\boldsymbol{\Phi}-{\bf
G}^{\Phi}_{k}\|_{F}^{2}$ (11)

where

$\displaystyle{\bf G}^{L}_{k}$ $\displaystyle={\bf
L}_{k}-\frac{\mu}{L_{f}}({\bf E}^{H}({\bf E}({\bf L}_{k}+{\bf S}_{k})-{\bf
y}))-\frac{1}{L_{f}}({\bf L}_{k}{\bf D}^{T}{\bf D}-{\boldsymbol{\Phi}}_{k}{\bf
D})$ $\displaystyle{\bf G}^{S}_{k}$ $\displaystyle={\bf
S}_{k}-\frac{\mu}{L_{f}}({\bf E}^{H}({\bf E}({\bf L}_{k}+{\bf S}_{k})-{\bf
y}))$ $\displaystyle{\bf G}^{\Phi}_{k}$
$\displaystyle={\boldsymbol{\Phi}}_{k}-\frac{1}{L_{f}}({\boldsymbol{\Phi}}_{k}-{\bf
L}_{k}{\bf D}^{T}).$

Here ${\bf E}^{H}$ denotes the adjoint operator of ${\bf E}$ which maps a
vector to a matrix. Each of these subproblems has closed-form solutions using
the proximal maps for nuclear norm or $\ell_{1}$ norm as follows:

$\displaystyle{\bf L}_{k+1}$
$\displaystyle={SVT}_{{\mu\lambda_{L}}/{L_{f}}}\left[{\bf G}^{L}_{k}\right]$
(12) $\displaystyle{\bf S}_{k+1}$ $\displaystyle={\bf
T}^{H}(\mathcal{S}_{{\mu\lambda_{S}}/L_{f}}\left[{\bf T}({\bf
G}^{S}_{k})\right])$ (13) $\displaystyle\boldsymbol{\Phi}_{k+1}$
$\displaystyle=\mathcal{S}_{{\mu\lambda_{D}}/{L_{f}}}\left[{\bf
G}^{\Phi}_{k}\right]$ (14)

where $SVT$ denotes the singular value thresholding operator for complex-
valued numbers given by ${SVT}_{\tau}[{\bf Z}]={\bf
U}\mathcal{S}_{\tau}[\boldsymbol{\Sigma}]{\bf V}^{H}$ where ${\bf Z}={\bf
U}\boldsymbol{\Sigma}{\bf V}^{H}$ is the singular value decomposition (SVD) of
a complex-valued matrix ${\bf Z}$, and $\mathcal{S}_{\tau}[{\bf Z}]$ denotes
the element-wise soft-thresholding or shrinkage operator at threshold
$\tau>0$, i.e., $\mathcal{S}_{\tau}[z]=\frac{z}{|z|}\max(|z|-\tau,0)$. Here
${\bf T}^{H}$ is adjoint operator of ${\bf T}$, defined by the corresponding
inverse Fourier transform. We terminate the iteration when $\left\|{\bf
L}_{k+1}+{\bf S}_{k+1}-({\bf L}_{k}+{\bf S}_{k})\right\|_{F}\leq
10^{-5}\left\|{\bf L}_{k}+{\bf S}_{k}\right\|_{F}$. The proposed PGM method
for solving the smoothness-regularized L+S decomposition and reconstruction of
dMRI is summarized in Algorithm 1.

Algorithm 1 Proximal Gradient L+S with Smoothness

0: Undersampled k-t data ${\bf y}$, data acquisition operator ${\bf E}$,
temporal Fourier transform ${\bf T}$, difference matrix ${\bf D}$,
regularization parameters $\lambda_{L}$, $\lambda_{S}$, $\lambda_{D}$, $\mu$,
$L_{f}$.

0: ${\bf M}_{0}={\bf L}_{0}={\bf E}^{H}{\bf y}$, ${\bf S}_{0}={\bf 0}$,
$\boldsymbol{\Phi}_{0}={\bf 0}$

1: while not converged do

2: ${\bf L}_{k+1}={SVT}_{{\mu\lambda_{L}}/L_{f}}\left[{\bf
L}_{k}-\frac{1}{L_{f}}(\mu{\bf M}_{k}+{\bf L}_{k}{\bf D}^{T}{\bf
D}-{\boldsymbol{\Phi}}_{k}{\bf D})\right]$

3: ${\bf S}_{k+1}={\bf T}^{H}(\mathcal{S}_{{\mu\lambda_{S}}/{L_{f}}}\left[{\bf
T}({\bf S}_{k}-\frac{\mu}{L_{f}}{\bf M}_{k})\right])$

4:
$\boldsymbol{\Phi}_{k+1}=\mathcal{S}_{{\mu\lambda_{D}}/{L_{f}}}\left[{\boldsymbol{\Phi}}_{k}-\frac{1}{L_{f}}({\boldsymbol{\Phi}}_{k}-{\bf
L}_{k}{\bf D}^{T})\right]$

5: ${\bf M}_{k+1}={\bf E}^{H}({\bf E}({\bf L}_{k+1}+{\bf S}_{k+1})-{\bf y}))$

6: end while

6: ${\bf L},{\bf S},\boldsymbol{\Phi}$

Table 1: Quantitative comparison of reconstruction using different methods on
PINCAT and cardiac perfusion data.

Dataset | Method | SER (dB) | PSNR (dB) | SSIM (n.u) | HFEN (n.u)
---|---|---|---|---|---
PINCAT | Zerofilled | 18.8056 | 31.2209 | 0.8069 | 0.5108
kt-SLR Lingala et al. (2011) | 18.8078 | 31.2230 | 0.8071 | 0.5108
L+S Otazo et al. (2015) | 24.5159 | 36.9767 | 0.9594 | 0.2513
L+S-FISTA Lin and Fessler (2018) | 24.6016 | 37.0687 | 0.9597 | 0.2465
L+S-POGM Lin and Fessler (2018) | 24.5765 | 37.0263 | 0.9597 | 0.2478
altGDminMRI1 Babu et al. (2023) | 25.4515 | 38.3994 | 0.9292 | 0.2378
SR-L+S-$\ell_{2}$ (ours) | 26.0206 | 38.5896 | 0.9437 | 0.2146
SR-L+S-$\ell_{1}$ (ours) | 27.7572 | 40.2958 | 0.9637 | 0.1711
Cardiac | Zerofilled | 6.1757 | 19.1928 | 0.5037 | 0.8448
kt-SLR Lingala et al. (2011) | -1.9176 | 10.9185 | 0.0977 | 1.3755
L+S Otazo et al. (2015) | 13.3158 | 26.2230 | 0.7874 | 0.4546
L+S-FISTA Lin and Fessler (2018) | 13.3129 | 26.2200 | 0.7868 | 0.4548
L+S-POGM Lin and Fessler (2018) | 15.3863 | 28.2842 | 0.8509 | 0.3243
altGDminMRI1 Babu et al. (2023) | 15.5309 | 28.3784 | 0.8678 | 0.2547
SR-L+S-$\ell_{2}$ (ours) | 16.3301 | 29.2059 | 0.8739 | 0.2715
SR-L+S-$\ell_{1}$ (ours) | 16.5811 | 29.4659 | 0.8787 | 0.2526

##  3 Experimental Results

###  3.1 Datasets

We evaluate the effectiveness of the proposed SR-L+S model for dMRI
reconstruction on two datasets. 1) the physiologically enhanced non-uniform
cardiac torso (PINCAT) numerical phantom data used in Lingala et al. (2011);
Lin and Fessler (2018) and 2) in vivo cardiac perfusion MRI data from Otazo et
al. (2015); Lin and Fessler (2018). The groud-truth PINCAT data has spatial
dimension $N_{x}\times N_{y}=128\times 128$ with $N_{t}=50$ temporal frames.
To emulate the multi-coil setting for the PINCAT data, we simulated
sensitivity maps of 32 coils (4 rings of 8 coils) compressed to $N_{c}=8$
coils, following setup in Lingala et al. (2011); Lin and Fessler (2018). The
cardiac data consist of images with size $N_{x}\times N_{y}=128\times 128$,
with $N_{t}=40$ temporal frames and $N_{c}=12$ coils with a retrospective
undersampling factor of 8 as in Lingala et al. (2011); Otazo et al. (2015);
Lin and Fessler (2018).

###  3.2 Performance metrics & comparison

To quantitatively assess the dMRI reconstruction performance of each method,
we utilized performance metrics commonly used in CS-based reconstructions.
These include signal error rate (SER), peak signal-to-noise ratio (PSNR),
structural similarity index (SSIM), and high-frequency error norm (HFEN)
Ravishankar and Bresler (2010). The reported PSNR and SSIM are averages over
all the frames. For HFEN, we employed a LoG filter with a kernel size of
$15\times 15$ pixels and a standard deviation of $1.5$ pixels. A higher SER,
PSNR and SSIM and lower HFEN indicate better reconstruction respective to the
original reference images.

We compared the proposed SR-L+S models with 5 state-of-the-art low-rank
methods for dMRI: kt-SLR Lingala et al. (2011), L+S with proximal gradient
method Otazo et al. (2015), L+S solved by fast iterative shrinkage-
thresholding (FISTA) algorithm Lin and Fessler (2018), L+S solved by proximal
optimized gradient method (POGM) Lin and Fessler (2018), and altGDminMRI1 Babu
et al. (2023) which used an alternating gradient descent combined with
minimization techniques assuming a hierarchical low-rank model. Zerofilled
reconstruction Bernstein et al. (2001) is also included as baseline.

###  3.3 Parameter settings

The regularization parameters were empirically selected on respective dataset
by choosing an parameter set that gives optimal reconstruction performance
over a range of values. For the PINCAT dataset, the SR-L+S-$\ell_{1}$ used
$\lambda_{L}$=0.01, $\lambda_{S}$=0.001, $L_{f}$=3, and $\lambda_{D}$=0.01;
the SR-L+S-$\ell_{2}$ used $\lambda_{L}$=0.01, $\lambda_{S}$=0.001, $L_{f}$=5,
and $\lambda_{D}$=0.01. For the Cardiac dataset, parameters
$\lambda_{L}$=0.01, $\lambda_{S}$=0.01, $L_{f}$=3, and $\lambda_{D}$=0.005 is
set for the SR-L+S-$\ell_{1}$, and $\lambda_{L}$=0.01, $\lambda_{S}$=0.01,
$L_{f}$=5, and $\lambda_{D}$=0.005 for the SR-L+S-$\ell_{2}$. We set $\mu=1$
for all experiments. For the competing methods, we followed the original
parameter settings provided by the authors for the implementation on the
PINCAT and cardiac data.

Figure 1: Comparison of reconstruction results of PINCAT synthetic data (14th
frame) with 8-fold undersampling using different methods. First row: Original
(fully-sampled) and reconstructed images. Second row: Enlarged views of yellow
boxed region. Third row: Error maps (display scale of [0, 0.2]) with
respective to the original image. Forth & fifth rows: y-t images at the
vertical cut line (slice $x=64$) in the reconstructed images, and the
corresponding error maps, respectively.

Figure 2: Comparison of reconstruction results of multi-coil cardiac perfusion
MRI with 8-fold undersampling using different methods. First row: Original
(fully-sampled) and reconstructed images. Second row: Enlarged views of
corresponding heart region with myocardial wall enhancement. Third row: Error
maps with respective to the original image. Forth & fifth rows: y-t images at
the vertical cut line (slice $x=64$) in the reconstructed images, and the
corresponding error maps, respectively.

###  3.4 Quantitative & qualitative results

####  3.4.1 Reconstruction accuracy

Table 1 shows the quantitative results of the proposed SR-L+S and comparison
methods measured by SER, PSNR, SSIM and HFEN. It is obvious that SR-L+S models
achieve the highest reconstruction accuracy compared to other methods, with
substantial improvements of all evaluation metrics on both datasets. The
highest SER, PSNR and SSIM, and lowest HFEN of our models indicate the minimal
level of noise and distortion in the reconstructed images, while preserving
the fine details and structural information of the original images. The
results imply that the incorporation of smoothness regularization of the
background components in SR-L+S can enhance dMRI reconstruction performance of
the original L+S formulations. Among the SR-L+S models, use of $\ell_{1}$
smoothness prior outperforms $\ell_{2}$ smoothness, suggesting the advantage
of piece-wise smoothness induced by the $\ell_{1}$-penalty. We can see that
L+S models perform better that the L&S model, with improved reconstruction by
the use of efficient algorithms such as FISTA, POGM and altGDminMRI1.

####  3.4.2 Visual results

Fig. 1 and Fig. 2 show the reconstruction results of different methods in the
spatial and time domains on the PINCAT and cardiac perfusion data,
respectively. The first and second rows show the reconstructed images (14th
frame) and magnified views or a particular region of interest (ROI) in the
yellow-box area. We can see that the proposed SR-L+S models provide visually
higher-quality reconstruction on both datasets. They reconstruct cleaner
images with lower aliasing artifacts and better preservation of sharp edges
and fine details of the images. Blurry reconstructions are evident
particularly for the kt-SLR and the original L+S models. The Kt-SLR also
suffers from noise enhancement on the cardiac data as also reported in Babu et
al. (2023). The lower residual artifacts of our method can be indicated
clearly in the error maps relative to the ground-truth images (third row). The
reconstruction errors of the other methods are more noticeable that that of
the proposed method, especially around the edges, e.g., of the papillary
muscle for the cardiac data. The y-t images (forth row) and its corresponding
error maps (fifth row) also show that our method can capture the dynamic
information more accurately. The reconstructed time profiles of competing
methods tend to over-smooth with more pronounced motion blurring especially on
the cardiac data. These qualitative findings consistently indicate the more
accurate and better reconstruction performance of our method , as corroborated
by the quantitative measures in Table 1.

(a)

(b)

Figure 3: The L and S component separation of cardiac perfusion dMRI sequence
of a patient with coronary artery disease. (a) The x-y view of the original
image (14th frame), the reconstructed image by ${\bf X}={\bf L}+{\bf S}$, and
the reconstructed ${\bf L}$ component and ${\bf S}$ component using different
L+S models. (b) The y-t views of reconstructed images by different models.

####  3.4.3 Separation of L and S components

Fig. 3 shows the reconstruction results on the cardiac data by ${\bf X}={\bf
L}+{\bf S}$, along with its recovered background ${\bf L}$ and dynamic ${\bf
S}$ components. From the x-y view (Fig. 3(a)), both the SR-L+S models provide
a better visualization of the perfusion defect (as indicated by an arrow) in
the sparse component ${\bf S}$ than the original L+S model, with the
background successfully suppressed. More contrast is observed between the
healthy part of the myocardium and the lesion. From the y-t view (Fig. 3(b)),
we can clearly see that the SR-L+S models present higher temporal fidelity
with respective to the original image sequence, with substantially reduced
temporal blurring artifacts. It is noteworthy that the our method can recover
${\bf L}$ that can better encapsulate the slow variation or temporal
correlations of the background structure, while revealing more pronounced
dynamic activities sparsely distributed in ${\bf S}$. In contrast, substantial
dynamic activities are still present in the background part of the images for
the original L+S model. This implies better separation of background and
foreground components of dMRI by our method, owing to the additional
smoothness regularization that can further smooth out the remaining noise and
dynamic activities from the background part. Nevertheless, we observe that the
use of $\ell_{2}$ smoothness prior in the SR-L+S model may in certain degree
over-smooth the background part compared to the $\ell_{1}$ penalty.

##  4 Conclusion

We developed a novel framework based on SR-L+R decomposition to reconstruct
dMRI from undersampled data. The SR-L+R is an extension of the L+S
formulation, that applies joint low-rank and smooth regularization to fully
exploit both the global and local temporal correlations in the slowly-varying
background component, enhancing reconstruction robustness to noise. An
efficient algorithm based on the proximal gradient method is developed to
solve a convex optimization for the SR-L+R decomposition. Experimental results
on dynamic cardiac and synthetic data shows that the proposed method can
provide better recovery of the dynamic MR images, and its background and
dynamic components, compared to existing methods. Future work could extend the
proposed model using tensor decomposition Roohi et al. (2017) and deep
unfolded network that utilizes deep networks to learn the regularization
parameters and the proximal mappings in the iterative solver as in Huang et
al. (2021); Liu et al. (2020).

## References

  * Lustig et al. [2007] Michael Lustig, David Donoho, and John M Pauly.  Sparse MRI: The application of compressed sensing for rapid mr imaging.  _Magn. Res. Med._ , 58(6):1182–1195, 2007.
  * Pruessmann et al. [1999] Klaas P Pruessmann, Markus Weiger, Markus B Scheidegger, and Peter Boesiger.  SENSE: Sensitivity encoding for fast MRI.  _Magn. Res. Med._ , 42(5):952–962, 1999.
  * Otazo et al. [2010] Ricardo Otazo, Daniel Kim, Leon Axel, and Daniel K Sodickson.  Combination of compressed sensing and parallel imaging for highly accelerated first-pass cardiac perfusion MRI.  _Magn. Res. Med._ , 64(3):767–776, 2010.
  * Lustig et al. [2006] Michael Lustig, Juan M Santos, David L Donoho, and John M Pauly.  kt SPARSE: High frame rate dynamic MRI exploiting spatio-temporal sparsity.  In _Proc. ISMRM_ , volume 2420, 2006.
  * Lingala et al. [2011] Sajan Goud Lingala, Yue Hu, Edward DiBella, and Mathews Jacob.  Accelerated dynamic MRI exploiting sparsity and low-rank structure: kt SLR.  _IEEE Trans. Med. Imag._ , 30(5):1042–1054, 2011\.
  * Trémoulhéac et al. [2014] Benjamin Trémoulhéac, Nikolaos Dikaios, David Atkinson, and Simon R Arridge.  Dynamic MR image reconstruction–separation from undersampled (k,t)-space via low-rank plus sparse prior.  _IEEE Trans. Med. Imag._ , 33(8):1689–1701, 2014\.
  * Otazo et al. [2015] Ricardo Otazo, Emmanuel Candes, and Daniel K Sodickson.  Low-rank plus sparse matrix decomposition for accelerated dynamic MRI with separation of background and dynamic components.  _Magn. Res. Med._ , 73(3):1125–1136, 2015.
  * Lin and Fessler [2018] Claire Yilin Lin and Jeffrey A Fessler.  Efficient dynamic parallel MRI reconstruction for the low-rank plus sparse model.  _IEEE Trans. Comput. Imag._ , 5(1):17–26, 2018\.
  * Babu et al. [2023] Silpa Babu, Sajan Goud Lingala, and Namrata Vaswani.  Fast low rank column-wise compressive sensing for accelerated dynamic MRI.  _IEEE Trans. Comput. Imag._ , 9:409–424, 2023.  doi:10.1109/TCI.2023.3263810.
  * Yao et al. [2018] Jiawen Yao, Zheng Xu, Xiaolei Huang, and Junzhou Huang.  An efficient algorithm for dynamic MRI using low-rank and total variation regularizations.  _Med. Imag. Anal._ , 44:14–27, 2018.
  * Wang et al. [2020] Dong Wang, David S Smith, and Xiaoping Yang.  Dynamic MR image reconstruction based on total generalized variation and low-rank decomposition.  _Magn. Res. Med._ , 83(6):2064–2076, 2020.
  * Ravishankar and Bresler [2010] Saiprasad Ravishankar and Yoram Bresler.  MR image reconstruction from highly undersampled k-space data by dictionary learning.  _IEEE Trans. Med. Imag._ , 30(5):1028–1041, 2010\.
  * Bernstein et al. [2001] Matt A Bernstein, Sean B Fain, and Stephen J Riederer.  Effect of windowing and zero-filled reconstruction of MRI data on spatial resolution and acquisition strategy.  _J. Magn. Res. Med._ , 14(3):270–280, 2001.
  * Roohi et al. [2017] Shahrooz F Roohi, Dornoosh Zonoobi, Ashraf A Kassim, and Jacob L Jaremko.  Multi-dimensional low rank plus sparse decomposition for reconstruction of under-sampled dynamic MRI.  _Pattern Recognition_ , 63:667–679, 2017.
  * Huang et al. [2021] Wenqi Huang, Ziwen Ke, Zhuo-Xu Cui, Jing Cheng, Zhilang Qiu, Sen Jia, Leslie Ying, Yanjie Zhu, and Dong Liang.  Deep low-rank plus sparse network for dynamic MR imaging.  _Med. Imag. Anal._ , 73:102190, 2021.
  * Liu et al. [2020] Risheng Liu, Yuxi Zhang, Shichao Cheng, Zhongxuan Luo, and Xin Fan.  A deep framework assembling principled modules for CS-MRI: unrolling perspective, convergence behaviors, and practical modeling.  _IEEE Trans. Med. Imag._ , 39(12):4150–4163, 2020\.