Bi-Linear Modeling of Data Manifolds for Dynamic-MRI Recovery
Gaurav N. Shetty Konstantinos Slavakis Abhishek Bose Ukash Nakarmi Gesualdo Scutari Leslie Ying







= 5em



.pdf






24bp24bp
    

10bp14bp
    







centercolon=true




OT1pzcmit







=2500

















@=.6ex plus1pt
algorithmic1



algorithmt
















compat=newest















































































select coords between index/.style 2 args=
    x filter/.code=
        <#1>#2













 

[1]{}
  #1 
[2]⟨⟩#1
  | #2
[1]‖‖#1 · #1











































√()√()
#1#2
0=#1#20=0
0-0.20
2=height0 depth -0
00.4pt2

kostasstyle
2pt
2pt

0pt

.
1ex


kostasstyle

theoremTheorem

subalgoSubalgorithm
factFact
claimClaim
propProposition
lemmaLemma
corollaCorollary
remarkRemark
definitionDefinition
exampleExample
assumptionAssumption
problemProblem
resultResult


















  noitemerr
  latex@warningEmpty `thebibliography' environment
  






(<ref>)










  1]Gaurav N. Shetty
1]Kostas Slavakis
2]Abhishek Bose
3]Ukash Nakarmi
4]Gesualdo Scutari
5]Leslie Ying
[1]Department of Electrical Engineering, University at Buffalo. Emails: {gauravna,kslavaki}@buffalo.edu.
[2]Department of Computer Science, INRIA, Université Côte d’ Azur. Email: abhishek.bose@inria.fr.
[3]Department of Electrical Engineering and Radiology, Stanford University. Email: nakarmi@stanford.edu.
[4]Department of Industrial Engineering, Purdue University. Email: gscutari@purdue.edu.
[5]Department of Electrical Engineering and Biomedical Engineering, University at Buffalo. Email: leiying@buffalo.edu.


Bi-Linear Modeling of Data Manifolds for
  Dynamic-MRI Recovery
    [
    
===============================================================

 
  This paper puts forth a novel bi-linear modeling framework for data recovery
  via manifold-learning and sparse-approximation arguments and considers its
  application to dynamic magnetic-resonance imaging (dMRI). Each temporal-domain
  MR image is viewed as a point that lies onto or close to a smooth manifold,
  and landmark points are identified to describe the point cloud concisely. To
  facilitate computations, a dimensionality reduction module generates
  low-dimensional/compressed renditions of the landmark points. Recovery of
  high-fidelity MRI data is realized by solving a non-convex minimization task
  for the linear decompression operator and affine combinations of
  landmark points which locally approximate the latent manifold geometry. An
  algorithm with guaranteed convergence to stationary solutions of the
  non-convex minimization task is also provided. The aforementioned framework
  exploits the underlying spatio-temporal patterns and geometry of the acquired
  data without any prior training on external data or information. Extensive
  numerical results on simulated as well as real cardiac-cine
  MRI data illustrate noteworthy improvements of the advocated machine-learning
  framework over state-of-the-art reconstruction techniques.






§ INTRODUCTION


Magnetic-resonance imaging (MRI), a non-invasive, non-ionizing and high-fidelity
visualization technology, has found widespread applications in cardiac-cine,
dynamic contrast-enhanced and neuro-imaging, playing a key role in medical
research and diagnosis <cit.>. MRI's inherent limitations
and various physiological constraints often incur slow data acquisition, while
long scanning times make MRI an expensive process, may cause patient discomfort,
and hinder thus its usefulness.

Raw MRI data are observed in the k-space domain; the image (visual)
data are computed by the inverse Fourier transform of the k-space
ones <cit.>. A prominent way to speed up data acquisition
is to sample the k-space or frequency domain densely enough to ensure
reconstruction of a high-fidelity and artifact-free image via Fourier-transform
arguments <cit.>. In the case of dynamic (d)MRI, where an
extra-temporal dimension is added to the spatial domain, sampling is also
performed along the time axis. Due to MRI's slow scanning times, it becomes
difficult for the data acquisition process to keep up with the motion of the
organs or the fluid flow in the field of view
(FOV) <cit.>. It is a usual case to not be able to reach
the necessary sampling density, not only in the temporal direction but also in
the k-space domain, to guarantee artifact-free reconstructed images. This
“under-sampling” inflicts signal aliasing and
distortion <cit.>.

Naturally, a lot of the MRI-research effort has been focusing on developing
reconstruction algorithms that improve the spatio-temporal resolution of MR
images given the highly under-sampled k-space data. To this end,
artificial-intelligence (AI) approaches have been very recently placed at the
focal point of MRI research; examples are convolutional neural networks
(CNNs) <cit.>, deep
variational networks(VNs) <cit.>, and generative
adversarial neural networks (GANs) <cit.>. AI methods learn non-linear mappings via extensive offline
training on large-scale datasets, different from or in addition to the acquired
data, and use those learned non-linear mappings to map the observed
low-resolution, or, undersampled data to their high-fidelity counterparts. In
contrast to the data-driven AI methods, the present work, as well as
the following prior-art schemes, assume no offline training on
large-scale datasets and resort solely to the observed data.

A popular approach to exploit the underlying spatio-temporal patterns within
dMRI data is compressed sensing (CS) <cit.>. Low-rank structures <cit.> and total-variation-based
schemes <cit.> have also
been explored at length and found to produce promising results for slow varying
dynamic data. For instance, <cit.> proposes the estimation,
first, of a temporal basis of image time series via singular-value
decomposition, prior to formulating a sparsity inducing convex-recovery
task. Nevertheless, these schemes seem to be less effective when it comes to
dMRI with extensive inter-frame motion or with a low number of temporal frames,
as mentioned in <cit.>. Acceleration capabilities of CS combined with the motion-robustness of radial
  imaging techniques have also found their place in MRI reconstruction and have
  shown promising results in cardiac and respiratory motion correction. The
  motion resolved strategy <cit.> introduces extra motion dimensions
  by sorting the continuously acquired k-space data into distinct motion states
  followed by a sparsity and total variation based CS technique for
  reconstruction. There has also been a growing interest in MRI-recovery by
dictionary-learning (DL) schemes <cit.>. In <cit.>, for example, the dMRI data
are decomposed in two components: A low-rank one, that captures the temporal
(video) background, and a sparse one, described via spatio-temporal patch-based
DL that models the (dynamic) foreground.

Manifold-learning techniques have also been employed to recover dMRI data from
highly under-sampled observations <cit.>. In <cit.>, a graph-Laplacian matrix is
formed via the Euclidean distances between points of a data cloud and is
subsequently used as a regularizer in a convex-recovery
task. Improving on <cit.>, a bandlimited modeling
  of the data points, coupled with an improved estimation of the Laplacian
  matrix, was proposed in <cit.> to make the reconstruction
  resilient to noise and patient motion, computationally inexpensive and less
  demanding memory-wise. The work in <cit.> proposes a kernel
  variation to the graph Laplacian estimation in <cit.>.  A
popular path followed by manifold-learning schemes is to perform dimensionality
reduction of the collected high-dimensional data prior to applying a
reconstruction algorithm, , <cit.>. In the
non-MRI context, <cit.> capitalizes, also, on a
Euclidean-distance-based Laplacian matrix to perform dimensionality reduction
prior to reconstructing data by local principal component analysis. In the
previous schemes, all of the observed points participate in the
dimensionality-reduction task, raising thus computational burdens, especially in
cases where the number of data is excessively large. Methods that perform
dimensionality reduction on properly chosen small-cardinality subsets of the
observed data cloud have been introduced for clustering and classification, but
not for regression tasks <cit.>.

This paper follows the manifold-learning path and serves a two-fold objective:
[label=*)]

  * Present a machine-learning framework, the bi-linear modeling of data manifolds (BiLMDM), that contributes novelties to exploiting local and latent data structures via a sparsity-aware and bi-linear optimization task; and

  * apply BiLMDM to the dMRI-data recovery problem.

In a nutshell, BiLMDM can be described as follows. Each vector of data, observed
from an undersampled dMRI temporal frame, is modeled as a point
onto or close to an unknown manifold, embedded in a Euclidean space. The only
assumption imposed on the manifold is smoothness <cit.>. Landmark points are chosen to concisely describe the
  observed data-vector cloud, and are mapped to a lower-dimensional space, as
  in <cit.>, to effect compression and enable
  low-computational footprints in the proposed algorithmic solutions. Motivated
  by the smooth-manifold hypothesis, the proposed work approximates each data
  vector as an affine combinations of neighboring landmark points
  ( Fig. <ref>). A locally
  bi-linear factorization model, novel for data representations, is then
formed to model/fit the point cloud: One factor gathers the coefficients of the
previous affine combinations of the landmark points, while the other one serves
as the linear decompression operator that unfolds the points back to
the image dimensions. Improving on our preliminary
results <cit.>, a highly modular bi-linear optimization task is
tailored to the dMRI-data recovery problem, penalized by terms which account for
sparsity along the temporal axis and other modeling assumptions. A
successive-convex-approximation algorithm is proposed to guarantee convergence
to a stationary solution of the previous bi-linear optimization task.

 The proposed work is validated
  against the following state-of-the-art schemes: Partially separable sparsity
  aware model (PS-Sparse) <cit.>, joint manifold learning and
  sparsity aware (MLS) framework <cit.>, smoothness
  regularization on manifolds (SToRM) <cit.>, low rank and
  adaptive sparse signal model (LASSI) <cit.> and extra
  dimensional golden angle radial sparse parallel
  (XD-GRASP) <cit.>. In spite of sharing with the previous
  state-of-the-art schemes the principle of identifying low-dimensional
  structures for high-dimensional data and exploiting this structure to
  reconstruct data, BiLMDM's original contributions in data approximations can
  be summarized as follows: 
  [label=*)]

  
  * BiLMDM departs from the mainstream manifold-learning approach of
    identifying a data-graph Laplacian matrix that penalizes an optimization
    (data-recovery) task as in <cit.>, and uses instead the fundamental geometric principle
    of tangent spaces of smooth manifolds to search for data
    dependencies/patterns through local affine combinations of landmark points
    ( Fig. <ref>);

  
  * it offers a novel bi-linear factorization model, where one linear factor
    captures the local smoothness of the manifolds, via local affine
    combinations, while the other one accounts for the linear decompression
    operator which maps the low-dimensional affine approximations back to the
    high-dimensional raw-data/input space. The affine combinations and the
    decompression operator are jointly identified via a highly modular
    non-convex data recovery task, which may accommodate also additional prior
    information, , periodicity over the temporal axis of the dMRI data. Unlike the
    traditional bi-linear schemes, BiLMDM doesn't follow the classical dictionary learning approach
    of decomposing the image into two components as discussed earlier and model the dynamic foreground of the image as a sparse combination of adaptive dictionary atoms 
    (which act as image basis) of MR image signals. 

  
  A more detailed description of the differences of BiLMDM with the
  state-of-the-art dMRI recovery schemes is deferred to
  Sec. <ref>. The efficacy of BiLMDM was observed on both
  synthetically generated and experimentally acquired cardiac cine MR data,
  undersampled via both Cartesian and radial trajectories. BiLMDM has
  consistently outperformed all the competing techniques and produced
  reconstructions which were sharp, free of undesirable artifacts, deformations
  and temporal bleeding.

The rest of the paper is organized as follows: Section <ref> describes,
in short, the dMRI acquisition scheme in the k-space
domain. Section <ref> details the BiLMDM's modeling assumptions,
while Sec. <ref> describes the algorithm to solve the proposed
non-convex minimization task. The extensive numerical tests of
Section <ref> showcase that BiLMDM outperforms state-of-the-art dMRI
recovery schemes. The manuscript summarises the numerical results in
  Section <ref> and concludes in
Section <ref>. Finally, the appendix gathers basic mathematical
facts and expressions that are essential for the implementation of the algorithm
described in Sec. <ref>.



§ BI-LINEAR MODELING OF DATA MANIFOLDS




 §.§ dMRI data description




MRI data 𝒴∈^N_p× N_f
( denotes the set of all complex-valued numbers) are observed in
k-space (frequency domain), which spans an area of size
N_p× N_f
( Figs. <ref> and
<ref>), with N_p standing for the number
of phase-encoding lines and N_f for the number of frequency-encoding
ones <cit.>. Data 𝒴 can be considered as
the two-dimensional (discrete) Fourier transform ℱ(·) of the
image-domain data
𝒳∈^N_p× N_f, ,
𝒴 =
ℱ(𝒳) <cit.>. Without any loss of
generality, this study assumes that the “low-frequency” part of
𝒴 is located around the center of the
N_p× N_f area. Availability of the data over the whole
k-space is infeasible in practice; k-space is usually severely
under-sampled <cit.>. There exist several strategies to
sample the k-space; examples are the 1-D Cartesian
(Fig. <ref>) and the radial
(Fig. <ref>) ones, where the “white” lines in
Figs. <ref> and
<ref> denote the available/sampled data, while data
in the “black” areas are not observed. A general trend among sampling
strategies is to put more emphasis on low-frequency components, which carry
contrast information and with high SNR, and select few high-frequency
components, which comprise high-resolution image details. The 1-D Cartesian
sampling pattern emulates the acquisition of k-space pixels via the 1-D Gaussian
distribution, acquiring a large number of samples in the central k-space area
while sampling few ones from the “high-frequency” area
( Fig. <ref>). The radial-sampling pattern
consists of radial spokes which yield dense sampling at the center of k-space,
while the sampling density is decreased as the spokes move away from the center
( Fig. <ref>).

In dMRI, an additional dimension is added to the MRI k-space to accommodate time (the axis vertical
on the paper in Fig. <ref>), resulting in the augmented (k,t)-space. The dMRI
(k,t)-space can be viewed, in other words, as the N_fr-fold Cartesian product of the
(N_p× N_f)-sized MRI k-space, where N_fr represents the number
of observed MRI frames over time. In dMRI, k-space 𝒴_j and image-domain
𝒳_j data are connected via 𝒴_j = ℱ(𝒳_j),
j∈1, …, N_fr. The (k,t) space is usually highly under-sampled. To extract
reliable information from the (k,t)-space data, this work follows <cit.> and considers a small number ν (≪ N_p) of phase-encoding lines, coined "navigator (pilot) data" (the gray-colored area in Fig. <ref>), to learn the
intrinsic low-dimensional structure of the data.

To facilitate processing, the (k,t)-space data are vectorized. More specifically,
(𝒴_j) stacks one column of 𝒴_j below the other to yield the
complex-valued N_k× 1 vector y_j :=
(𝒴_j) where N_k := N_pN_f denotes the number of pixels in every frame. To avoid notation clutter, ℱ still denotes the
two-dimensional (discrete) Fourier transform even when applied to vectorized versions of image
frames:
ℱ[(𝒳_j)] := [ℱ(𝒳_j)]=
(𝒴_j). All vectorized k-space frames are gathered in the
N_k× N_fr matrix
Y := [y_1, y_2, …, y_N_fr] so that the vectorized
image-domain data are
X := ℱ^-1(Y) := [ℱ^-1(y_1),
ℱ^-1(y_2), …, ℱ^-1(y_N_fr)], where
ℱ^-1(·) denotes the inverse two-dimensional (discrete) Fourier transform. The
navigator data of the jth k-space frame ( Fig. <ref>),
j ∈{1, 2, …, N_fr}, are gathered into a ν N_f× 1 vector
y_j^nav. All navigator data comprise the ν N_f× N_fr
matrix
Y_nav := [y_1^nav, y_2^nav, …,
y_N_fr^nav].





 §.§ Modeling assumptions


The high-dimensional navigator data y_j^nav_j=1^N_fr carry
useful information about spatio-temporal dependencies in the (k,t)-space. To promote parsimonious
data representations, especially in cases where N_fr attains large values, it is
desirable to extract a subset
ℓ_k_k=1^N_ℓ⊂y_j^nav_j=1^N_fr
(N_ℓ≤ N_fr), called landmark points, which provide a “concise
description,” in a user-defined sense, of the data cloud
y_j^nav_j=1^N_fr. To this end, the following assumption, often
met in manifold-learning approaches <cit.>, imposes structure on
y_j^nav_j=1^N_fr.


  Data y_j^nav_j=1^N_fr lie on a smooth low-dimensional
  manifold M <cit.> embedded in the high-dimensional Euclidean space
  ^ν N_f ( Fig. <ref>).


For example, the most well-known case of M is a linear
subspace, which is the main
  hypothesis behind principal component analysis (PCA). PCA assumes that the
  data lie close to or onto a linear subspace in a low-dimensional space. Along
  the same lines, based on As. <ref> and the concept of the
tangent space of a smooth manifold, it is conceivable that neighbouring landmark
points cooperate affinely to describe vector y_j^nav (the
gray-colored area in Fig. <ref> depicts all
possible affine combinations of
ℓ_k_1, ℓ_k_2, ℓ_k_3). In other words, every temporal frame is approximated by a combination of landmark points, where these landmark points
  can be viewed as a representative of a particular group of frames sharing
  similar phase and motion characteristics.  Upon defining the
ν N_f× N_ℓ matrix
Λ := [ℓ_1, ℓ_2, …, ℓ_N_ℓ],
the above notion can be established
  by the existence of an N_ℓ× 1 vector b_j
that renders the approximation error
y_j^nav - Λb_j small, where
· denotes the standard Euclidean norm of space
^ν N_f. Since affine combinations are desirable,
b_j is constrained to satisfy
1_N_ℓ^⊺b_j = 1, where 1_N_ℓ
stands for the all-one N_ℓ× 1 vector and superscript ⊺
denotes vector/matrix transposition. Moreover, motivated by the low-dimensional
nature of M (As. <ref>), only a few landmark points
are considered to cooperate to represent y_j^nav in its closed
vicinity, , b_j is sparse. The previous arguments are summarized
into the following modeling hypothesis.


  There exist a sparse N_ℓ× N_fr matrix B, with
  1_N_ℓ^⊺B = 1_N_fr^⊺, and a
  ν N_f× N_fr matrix E_2, which gathers approximation errors, such that
  (s.t.) Y_nav = ΛB + E_2.


The previous assumption holds true in the prototypical case where
M is a linear subspace and Y_nav comprises column
vectors which lie close to or onto M. Any matrix which includes as
column vectors any basis of M may serve as Λ. Since the
affine hull <cit.> of the columns of Λ
coincides with M, there exists surely a coefficient matrix
B that satisfies the affine constraints
1_N_ℓ^⊺B =
1_N_fr^⊺ and
Y_nav = ΛB. It is not now difficult to see that
  As. <ref> covers also the case where M
  is a union ∪_q=1^Q M_q of linear subspaces
  M_q_q=1^Q, with Y_nav comprising columns
  vectors which lie close to or onto M_q_q=1^Q. Hence,
  As. <ref> lays the foundations for more general
  cases where data vectors lie close to or onto a union of smooth manifolds,
  with different dimensions. In such a way, As. <ref>
  offers theoretical support and justification in cases where data comprise
  disparate groups, where intra-group dependencies can be captured by a single
  smooth manifold, while inter-group relations can be only viewed via the union
  of those manifolds, , the union of data drawn from a static gray-colored
  background and data describing the movement of a beating heart.

Although several strategies may be implemented to identify the landmark points
Λ, a greedy optimization methodology, introduced in
<cit.>, is adopted here. In short, at every step of the
algorithm, a landmark point is selected from
y_j^nav_j=1^N_fr that maximizes, over all
un-selected y_j^nav_j=1^N_fr,
the minimum distance to the landmark points which have been already selected up
to the previous step of the algorithm. The algorithm of <cit.>
scores a computational complexity of order
𝒪(N_ℓ N_fr), which is naturally heavier than that of a
procedure that selects ℓ_k_k=1^N_ℓ randomly from
y_j^nav_j=1^N_fr.

Still, the landmark points (columns of Λ) are high dimensional. To
meet restrictions imposed by finite computational resources, it is desirable to
reduce the dimensionality of Λ. To this end, the methodology of
<cit.>, which is motivated by <cit.>, is employed. The approach comprises of two steps:


  

  
  * Given Λ and a user-defined λ_W>0, solve
    
    min_W∈^N_ℓ×
          N_ℓ   *Λ - ΛW_F^2 + λ_W W_1 
    s.to   1_N_ℓ^⊺W =
            1_N_ℓ^⊺ and (W)
            = 0 ,

    where ·_F stands for the Frobenius norm of a matrix. Since ℓ_k_k=1^N_ℓ lie on
    the manifold M, then according to As. <ref> and
    Fig. <ref>, any
    point taken from ℓ_k_k=1^N_ℓ may be faithfully approximated by an
    affine combination of the rest of the landmark points. In other words, there exists a matrix
    W s.t. Λ≈ΛW. With
    1_N_ℓ^⊺W = 1_N_ℓ^⊺ manifesting
    the previous desire for affine combinations, the constraint (W) = 0 is used
    to exclude the trivial solution of the identity matrix I_N_ℓ for
    W. Task (<ref>) is an affinely constrained composite convex minimization
    task, and, hence, the framework of <cit.> can be employed to solve it, due to
    the flexibility by which <cit.> deals with affine constraints when compared with
    state-of-the-art convex optimization techniques.

  
  * Once W is obtained from the previous step and for a user-defined integer number
    d ≤min{N_ℓ, ν N_f}, solve
    
    min_Λ̌∈^d×
          N_ℓ   Λ̌ -
            Λ̌W_F^2 s.toΛ̌Λ̌^ = I_d ,

    where the constraint Λ̌Λ̌^ = I_d is
    used to exclude the trivial solution of Λ̌ = 0, and the superscript
     denotes the Hermitian transpose of a matrix. The solution of the previous task is
    nothing but the complex conjugate transpose of the matrix which comprises the d minimal
    eigenvectors of
    (I_N_ℓ - W) (I_N_ℓ - W)^.
  


To this point, the identified landmark points were used to capture a latent structure of the MR data cloud and the structure was further used to compress the landmark points to a lower dimension, thus, facilitating the search for affine combinations in a lower dimensional space. However, there arises a need to unfold the points back to the dimensions of the navigator data which leads us to the following model.  


  There exist an N_k× d matrix G_3 and an
  N_k× N_ℓ matrix E_3, which gathers all approximation
  errors, s.t. Λ = G_3 Λ̌ + E_3.


Matrix G_3 can be viewed as the “decompression” operator
which reconstructs the “full” Λ from its low-dimensional
representation Λ̌.


The relationship between the
  undersampled k-space frames and navigator frames can be established, for the
  case of Fig. <ref>, as 
Y_nav = ΩY, where Ω is a matrix
with binary entries {0, 1} that select entries of Y. The inverse
problem of recovering Y from its partial Y_nav is
viable under errors, , an approximation of Y can be obtained via
Ω^†Y_nav, where the decompressor
Ω^† stands for the Moore-Penrose pseudo-inverse of
Ω. The following modeling hypothesis generalizes this argument.


  There exist an N_k×ν N_f matrix G_1 and an
  N_k× N_fr matrix E_1, which gathers all
  approximation errors, s.t. Y = G_1 Y_nav +
  E_1.


A similar modeling assumption can be found in the classical principal component analysis
(PCA) <cit.>, where G_1 serves as the “decompression” operator
(usually an orthogonal matrix) that maps the “compressed” Y_nav back to the
original data Y under the approximation error E_1. However, PCA identifies
best-linear-fit approximations to data, which might be an inappropriate modeling assumption for
several types of data <cit.>. The following discussion departs from PCA and
establishes modeling assumptions that allow non-linear data geometries.


Putting modeling assumptions <ref>,
<ref> and <ref> together, it can be verified that there exist matrices
G and E s.t. Y = GΛ̌B + E. Upon
defining U:= ℱ^-1(G), and since
GΛ̌B = ℱ(U) Λ̌B =
ℱ(UΛ̌B), by virtue of the linearity of ℱ,
the following bi-linear model between Y and the unknowns (U, B) is
established:

    Y = ℱ(UΛ̌B)
      + E .

Bi-linearity means that if U (or B) is fixed to a
specific value, then Y is linear with respect to B (or U), modulo the error E term. Interestingly, the linearity of ℱ^-1 suggests that
the previous modeling hypothesis holds true also in the image domain:
ℱ^-1(Y) = UΛ̌B + ℱ^-1(E).



§ THE BI-LINEAR RECOVERY TASK AND ITS ALGORITHMIC SOLUTION
 

In practice, only few (k,t)-space data are known. To explicitly take account of the limited number
of data, a (linear) sampling/masking operator 𝒮(·) is introduced, where
𝒮(Y) leaves the entries of Y as they are at sampled or observed
positions of the k-space domain while nullifying all the rest. The sampling operator is capable
of mimicking any sampling strategy, such as Cartesian, radial, spiral, . It is also often in
dMRI that image frames capture a periodic process, , heart movement, other than the static
background. In other words, it is reasonable to assume that in (<ref>), the
one-dimensional Fourier transform ℱ_t of the 1× N_fr time profile of
every one of the N_k pixels, , every row of the matrix
ℱ_t(UΛ̌B), is a sparse vector.

All of the previous modeling assumptions are incorporated in the following recovery task: Given
the positive real-valued parameters λ_1, λ_2, λ_3, C_U, solve

    min_(U, B, Z)    12*𝒮(Y) -
        𝒮ℱ(UΛ̌B)_F^2^T1 
        + λ_12*Z -
        ℱ_t(UΛ̌B)_F^2^T2
        + λ_2 *Z_1_T3 + λ_3
        *B_1_T4
    s.to    *Ue_i≤ C_U, ∀ i∈1, …,
        d_C1; 1_N_ℓ^⊺B = 
        1_N_fr^⊺_C2; 
       U∈^N_k× d; B∈^N_ℓ× N_fr; Z∈^N_k×
        N_fr ,

where e_i denotes the ith column of the identity matrix
I_d. Elaborating more on the recovery task (<ref>), T1
corresponds to the data-fit term, while T2 and T3 introduce the auxiliary
variable Z, used to impose a sparsity constraint on
ℱ_t(UΛ̌B). T4 imposes a sparsity
constraint on B, following the discussion on
As. <ref>. Bound C_U in C1 is used to prevent
unbounded solutions for U due to the scaling ambiguity in the bi-linear
term UΛ̌B. C2 adds the affine constraint
discussed in As. <ref>. 
  Moreover, it is worth noticing here that the reduction of dimensionality,
  achieved via Λ̌, reduces also the number of columns and
  unknowns, , degrees of freedom, of U.






The proposed
  scheme shares with PS-Sparse <cit.>,
  SToRM <cit.>, (and its other variants bandlimited-SToRM <cit.>, navigator-less SToRM <cit.>) and MLS <cit.> similarities in
  terms of using navigator lines to learn the underlying data manifold. Instead
  of using the eigen vectors of the co-variance (PS-Sparse) and Laplacian (SToRM and its variants) matrices as the temporal basis, BiLMDM relies on robust sparse
  embedding <cit.> to estimate a low dimensional latent structure. SToRM
  and BiLMDM both aim at learning the structure of “smooth” data
  manifolds. Smoothness in SToRM is established globally over the data cloud by
  penalizing the recovery task via several operator norms of a graph-Laplacian
  matrix. On the other hand, BiLMDM enforces smoothness locally over the
  manifold via neighborhoods and affine approximations (patches) of the tangent
  spaces (Fig. <ref>), and the data manifold can
  be viewed as an approximation of the union of all those patches (as discussed
  in As. <ref>). PCA-inspired schemes, such as
  PS-Sparse, search also for local affine/linear patches to
  manifolds. Nevertheless, BiLMDM does not impose any orthogonality constraints
  on any of its factors U and B, as PCA via the singular
  value decomposition (SVD) does. Identification of landmark points to learn an underlying manifold is a novel contribution of the proposed work, given that all of the prior-art methods rely on the entire data cloud for learning. The incorporation of the landmark points in the bi-linear factorization model UΛ̌B
  serves also as a natural encapsulation of the underlying data geometry into
  the optimization task, extending thus the usual way that PCA interprets its
  factors as an orthogonal/decorrelating basis and coefficients. Given that MLS
  employs the same approach of robust sparse embedding, the key feature of
  BiLMDM that sets itself apart from MLS, as well as other manifold-learning
  methods, is the bi-linear factorization model used to represent the MR data
  cloud. This bi-linear factorization model helps to capture the geometry of the
  point cloud locally by imposing a sparsity-promoting penalty on the affine
  combinations, thus, restricting the sharing of data among frames that show
  common phase and structural characteristics; local exploitation of
  dependencies is promoted over any global one. Such a bi-linear model differs also from LASSI <cit.>, which adopts the popular approach of viewing the data matrix as the superposition of a low-rank and a sparse component. The low-rank component describes the static image background, while the sparse component, modeled via a classical bi-linear dictionary-learning term, describes the dynamic foreground. In contrast to LASSI, BiLMDM uses bi-linear modeling to describe the dynamic foreground and the static background together. BiLMDM does not also require any
  pre-processing procedures, like sorting the acquired data into cardiac and
  respiratory phases, as in the motion resolved strategies of
  XD-GRASP <cit.>, before solving the reconstruction
  task.

The successive-convex-approximation framework of <cit.> is
employed to solve (<ref>) and is presented in a concise form in
steps <ref>–<ref> of
Alg. <ref>. Convergence to a stationary solution of
(<ref>) is
guaranteed <cit.>. Step <ref> of
Alg. <ref> comprises convex minimization sub-tasks. More
specifically, at every step of the algorithm, given
(U_n, B_n, Z_n), the following estimates are required
(for τ_U, τ_B>0):

  
    Û_n∈_U    12*𝒮(Y) -
          𝒮ℱ(UΛ̌B_n)_F^2 + τ_U2*U - U_n_F^2 
        + λ_12*Z_n - ℱ_t(UΛ̌B_n)_F^2 
    s.to     *Ue_i≤ C_U,∀ i∈1, …, d . 
    B̂_n∈_B    12*𝒮(Y) - 𝒮ℱ(U_n
          Λ̌B)_F^2 + τ_B2*B-
          B_n_F^2 
        + λ_12*Z_n - ℱ_t(U_n Λ̌B)_F^2 + λ_3*B_1
          
    s.to    1_N_ℓ^⊺B =
          1_N_fr^⊺ .


Both tasks in (<ref>) can be viewed as affinely constrained composite convex
minimization tasks, hence allowing the use of <cit.>, as described in
Alg. <ref>. From a computational complexity perspective, it is worth pointing
out that the proposed scheme relies on minimization sub-tasks, and computational complexities depend
thus on the solver adopted for solving those sub-tasks. Here, <cit.> employs only
first-order information (gradients) and proximal mappings, , soft-thresholding rules and
projection mappings. The implementation of <cit.> for the specific tasks
(<ref>) is presented in Alg. <ref>, and details are deferred
to the appendix section of this manuscript. Notice also that the previous minimization sub-tasks
can be solved in parallel. Furthermore, problems <ref> can be solved inexactly at each iteration (with increasing precision, as described in <cit.>), which contributes to reducing the computation cost of each iteration; we refer to <cit.> for more details.





§ NUMERICAL RESULTS





  x       bi          error       PS          UK          ST          LASSI     XDGRASP
  4x      0.0438      0.00034     0.0437      0.0437      0.0526      0.0539    0.049
  8x      0.0454      0.00031     0.0449      0.0449      0.064       0.0552    0.0505
  12x     0.0461      0.00084     0.0462      0.0463      0.0762      0.0630    0.0521
  16x     0.0478      0.00079     0.048       0.048       0.085       0.0716    0.0591
  20x     0.0488      0.00025     0.055       0.0512      0.0937      0.0732    0.0653
  24x     0.0499      0.00083     0.0752      0.072       0.0984      0.0960    0.0693


  x       bi          error          PS          UK          ST       LASSI
  5x      0.1033      0.00031        0.2939      0.2680      0.2224   0.2287
  9x      0.2283      0.00029        0.3118      0.3400      0.3447   0.2929
  13x     0.2525      0.00025        0.3274      0.3459      0.3877   0.3543
  18x     0.2931      0.00027        0.3450      0.3589      0.4464   0.3832
  24x     0.3235      0.00021        0.3725      0.3714      0.5012   0.4052


  x       bi          error         PS          UK          ST        LASSI     XDGRASP
  6x      0.0488      0.00011        0.0607      0.0632      0.0530    0.0721   0.0535
  12x     0.0507      0.00015        0.0752      0.0791      0.0558    0.0759   0.0561
  17x     0.0527      0.00011        0.0863      0.0892      0.0583    0.0791   0.0587
  21x     0.0592      0.00012        0.0889      0.0932      0.0606    0.0849   0.0631
  34x     0.0636      0.00013        0.0959      0.1029      0.0668    0.0913   0.0674


















The proposed framework is tested and validated on three datasets, used also
in <cit.>:
[label=*)]

  * Magnetic-resonance extended cardiac-torso (MRXCAT) cine phantom <cit.>;

  * cardiac phantom generated from real MR scans <cit.>

  * prospectively undersampled free breathing, real time, cardiac cine data.
 
All experiments were conducted on a 12-core Intel(R) 2.40 GHz Linux-based system with 48GB
RAM, with implementations realized in MATLAB <cit.>. The proposed method is compared
with: PS-Sparse <cit.>, MLS <cit.>, SToRM <cit.>, LASSI <cit.> and XD-GRASP <cit.>.

The publicly available MATLAB implementations of SToRM <cit.>,
LASSI <cit.> and XD-GRASP <cit.> were
employed. MATLAB code was also written to realize the algorithmic solutions
of <cit.> to the convex-optimization problems posed in PS-Sparse
and MLS. Parameters were tuned to produce the least NRMSE for each algorithm at
each sampling ratio and on every dataset. The proposed framework currently underperforms
  in terms of computation times versus the competing algorithms and their
  optimized software implementations. Since the present manuscript serves as a
  proof of concept of BiLMDM, an optimized and time-efficient C/C++ version of
  the developed MATLAB code falls beyond the scope of this work and it will be
  presented in the near future via an open-source software-distribution venue.

Tasks (<ref>),
  (<ref>) and (<ref>) require parameter tuning to
  produce the desirable MR images from the scanner acquired data. These
  parameters were chosen empirically to minimize the reconstruction error for
  the datasets where the gold standard was known. A range of values for every
  parameter was realized, and as long as these parameters are within that range,
  the NRMSEs were not so sensitive to the parameters. This empirical research
  into the behaviour of the parameter values revealed the following: for hyper
  parameters penalizing the periodicity along the time series (λ_1,
  λ_2), higher the ratio λ_2/λ_1, severe temporal bleeding
  was observed; hence, it should be low enough to avoid temporal bleeding and at
  the same time exploit the periodicity along the temporal axis.  The parameter
  which penalizes the sparsity on B, λ_3, controls the size of
  the neighborhood used for approximation of the time series image. Higher the
  value of λ_3, lower the size of the neighborhood, allowing thus frames
  with very similar characteristics to share data among themselves. However,
  decreasing the value excessively would result in dissimilar frames influencing each other. Parameter N_ℓ determines the number
  of landmark points which act as a representative of a collection of frames
  with similar phase and motion characteristics. The value for N_ℓ
  should be high enough to describe the data cloud concisely. For the datasets
  used in this work, desirable values range between 15%-20% of
  N_fr. Any increase N_ℓ beyond that didn't yield any
  significant changes in the image quality. Bounding constant C_U set
  to 1 usually yields adequate estimates for U, although it is worth
  pointing out that the estimates for U are not sensitive to a change
  in C_U. The step sizes for convergence of Alg. <ref>
  (γ_0=0.9, ζ_0=0.001) and the step size for
  Alg. <ref> (α = 0.5) ensures desirable
  results. These values were reached by observing the effects of these
  parameters on the NRMSE values for the reconstruction of the two synthetically
  generated cardiac cine dataset.

Quality of reconstruction is evaluated by the normalized-root-mean-square error (NRMSE), defined as

    NRMSE := X - X̂_FX_F ,

where X represents the fully sampled, high-fidelity and original
image-domain data, while X̂ represents an estimate of X
computed by the reconstruction schemes. In
  addition to the NRMSE, the effectiveness of the algorithms are also inspected
  along the edges using the high-frequency error norm (HFEN) and sharpness
  measures. HFEN is defined as

    HFEN = LoG(X) -
      LoG(X̂)_2LoG(X)_2 ,

where LoG is a rotationally symmetric Laplacian of Gaussian filter. As
in <cit.>, a filter kernel of size 15 × 15 with
standard deviation of 1.5 pixels is employed. Two sharpness measures are used:
M1 (intensity variance based) and M2 (energy of the image gradient based), as
described in <cit.>. Better insights into the structural
information are extracted using the Structural Similarity (SSIM)
index <cit.> which investigates the local patterns in the pixel
intensities after normalizing for luminance and contrast. The metrics (except
for sharpness measure) can be used for assessment only when the ground truth
images are present and hence are used for the two phantom datasets and not the
prospectively undersampled data.

The proposed framework is validated over a range of
undersampling/acceleration rates, defined by
N_k N_fr / (# of acquired voxels). 1-D Cartesian
(Fig. <ref>) as well as radial
(Fig. <ref>) sampling were applied to 
both the phantom datasets. To save space, only the 1-D
Cartesian-sampling results are demonstrated for the MRXCAT dataset, while
radial-sampling ones are shown for the phantom generated
using real MR scans. Nevertheless, BiLMDM's performance against the
  competing reconstruction algorithms follows a similar trend also for sampling
  strategies not included in the manuscript due to space limitations.










 §.§ MRXCAT phantom


The extended cardiac torso (XCAT) framework was used to generate the MRXCAT
phantom <cit.>; a breath-hold cardiac cine data of spatial size
(N_p,N_f) = (408, 408) corresponding to a spatial resolution of
1.56 × 1.56 mm^2 for a FOV of 400 × 400 mm^2. The
cardiac cine phantom generated is spread across N_fr=360 time frames in the temporal
direction, consisting of 15 cardiac cycles and 24 cardiac phases. The MRXCAT-phantom dataset is
characterized by its limited inter-frame variations and periodic nature along the temporal
direction, unlike the other two datasets.

On the quantitative front, Fig. <ref> provides quantitative comparisons in terms of NRMSE for a range of acceleration factors; while  Fig. <ref> provides the fluctuations in the NRMSE across the time frames in addition to Tab. <ref> which summarizes all the performance metrics for an acceleration factor of 20x. It is evident from Fig. <ref>
that the proposed BiLMDM consistently outperforms the state-of-the-art schemes over the entire range
of acceleration rates. It is worth noticing here that BiLMDM scores similar NRMSE values to
PS-Sparse and MLS, at specific acceleration rates, even though it uses a small subset of the
navigator data (landmark points) in this case N_ℓ = 50, d = 4, in contrast to PS-Sparse and MLS that utilize the whole set of navigator data. BiLMDM's performance is consistent for every frame in the time series as displayed in Fig. <ref>, and exhibits the least NRMSE and least NRMSE deviation for every frame, unlike XD-GRASP, LASSI and SToRM. Tab. <ref> provides evidence that BiLMDM not only produces result closest to the ground truth, as exhibited by both the NRMSE and SSIM values, but also produces the sharpest image of all indicated by the smallest HFEN value and the largest M1 measure supported by a high M2 measure.

Shifting focus towards the qualitative aspect, Fig. <ref> focuses on a dynamic region of interest (ROI) from the end-diastolic and
end-systolic phases. The proposed framework exhibits reconstruction of a high-resolution image with
contrast close to the gold standard, as opposed to the alias-infected reconstruction of XD-GRASP and SToRM or the
noisy reconstructions of PS-Sparse, MLS and LASSI. The error maps in Figs. <ref>
and <ref> provide visual proof of the improvements in reconstruction. BiLMDM
shows significant improvements at the image edges, as visualized in the error maps for both the
diastolic and systolic phases and further supported by HFEN, M1 and M2 values. Moreover, PS-Sparse, MLS and LASSI exhibit faint-right artifacts (marked by red arrows in Fig. <ref>),
unlike BiLMDM. The temporal cross sections in Fig. <ref> are consistent with the
spatial results in Fig. <ref>. Reconstructions from XD-GRASP and SToRM exhibit temporal blurring, while LASSI reconstructions show presence of noisy artifacts in Fig. <ref>. SToRM, LASSI and XD-GRASP also witness breaks in the periodicity of the cardiac phases as pointed by red arrows in Fig. <ref>.







 §.§ Cardiac phantom generated from real MR scan


Besides synthetically generated MR images, real human-cardiac cine MR data were also
acquired <cit.> with the following parameters: spatial size
(N_p,N_f) = (200,256), FOV = 273 × 350 mm^2,
resulting in a spatial resolution of 1.36 × 1.36 mm^2. The data were acquired during a
single breath hold. Multiple time-wraps (introducing temporal variations) and quasi-periodic spatial
deformation <cit.> (to model respiration) were then used to generate a temporal sequence
with N_fr = 256 frames. The results in Fig. <ref>, similar to
those on the previous cine dataset, show that BiLMDM outperforms the rest of the techniques with
NRMSE values ranging from 0.05 to 0.061 for N_ℓ=35, d=12. Fig. <ref>, further supports the ability of BiLMDM to create reconstructions with the least NRMSE fluctuations across the many frames. As seen in Tab. <ref>, BiLMDM exhibits the least HFEN, M1 and M2 measure values which are testament to the fact that reconstructions obtained are sharper and are precise around the edges. The high similarity index proves that the BiLMDM is capable of producing results closest to the ground truth.  

With regards to Figs. <ref> and <ref>, significant
blurring and deformations can be observed for PS-Sparse, MLS and LASSI, in contrast to the sharper
images produced by BiLMDM for both the end-diastolic and end-systolic phases. There are slight
improvements over SToRM and XD-GRASP, which are illustrated via the error maps. At the same time, it is worth noting that in the end-systolic phase, PS-Sparse, MLS and LASSI produce an enlargement (deformation) of the area marked with a red rectangle on the image. The temporal cross sections in Fig. <ref>
indicate significant motion blurring for PS-Sparse, MLS and LASSI. Moreover, it can be verified that
a large number of artifacts are present in the temporal cross section produced by LASSI. Significant temporal bleeding is observed in results from PS-Sparse and MLS which is consistent with the undesirable enlargement in the end-systolic phase discussed earlier. Even though SToRM and XD-GRASP provide reconstructions almost as good as BiLMDM, it can be seen that there are signs of temporal blurring in the marked areas in Fig. <ref> for the corresponding methods. On the
contrary, BiLMDM provides a temporal cross-section very close to the gold standard.



 §.§ Prospectively Undersampled Cardiac Cine Data


The prospectively undersampled real-time free breathing cardiac cine data was acquired using a FLASH sequence from a volunteer breathing normally under the following acquisition parameters: TR/TE = 5.8/4 ms, FOV 284×350 mm, spatial resolution = 1.8 mm. A 12 channel scanner was used to continuously acquire 4500 phase encoding lines under 1-D Cartesian trajectory. These phase encoded lines were divided in groups of 15 (including 5 navigator lines) to form 300 frames, resulting in a data matrix of size 156×192×300 for each channel. The above acquisition corresponds to a undersampling factor of nearly 10x. For fairness in comparison, the multi-channel data was reconstructed coil by coil for all methods and then combined with the sum-of-squares strategy.

Fig. <ref> shows the reconstruction results from the proposed and state-of-the methods. The validation skips most of the performance metrics except for the sharpness measure due to the absence of a ground truth. Hence the validation is solely dependent on the visual quality as shown in Fig. <ref> and the sharpness measures. The full FOV images exhibit that BiLMDM along with SToRM, unlike other methods show well-defined structures and result in a sharper image. LASSI exhibits some aliasing effects, while PS-Sparse and MLS observes some blurring in the spatial frame. This is further supported by sharpness measures (M1 followed by M2): PS-Sparse (8.0×10^-8, 81.7), MLS (8.2×10^-8, 79.9), SToRM (8.6×10^-7, 91.8), LASSI (8.2×10^-8, 83.6), XD-GRASP (8.2×10^-7, 80.6) and BiLMDM (8.8×10^-7, 91.7). Comparing the temporal cross sections for the reconstructed images, even though the FOV images look pretty sharp for SToRM and BiLMDM, SToRM suffers from more temporal blurring than BiLMDM, while LASSI and PS-Sparse show grainy artifacts in the temporal cross section.






















§ DISCUSSION


The proposed BiLMDM scheme relies on data
  sharing among different frames by approximating each frame by an affine
  combination of its neighboring frames on the manifold that describes locally
  the MR data cloud. This manifold and sparse-approximation based recovery model
  enables the reconstruction of desirable, good quality and artifact free MR
  images from severely undersampled k-space measurements. The quality of the
  reconstructed images and the corresponding performance in numbers discussed
  under Sec. <ref> is a testament to the fact that the BiLMDM
  outperforms the competing state-of-the-art methods. BiLMDM has consistently
  provided better results for both synthetically generated, retrospectively
  undersampled as well as experimentally acquired, prospectively undersampled
  cardiac cine data.

  Methods like PS-Sparse and MLS have reconstructed desirable MR images in the
  case of the MRXCAT phantom, under cartesian sampling trajectories; however,
  they suffered severe motion blurring and deformations (enlargements) in the
  case of the cardiac cine data retrospectively sampled via radial
  trajectories. In contrast, SToRM produced good quality images for the data
  using radial sampling but suffered severe spatial blurring and temporal
  bleeding in the cases of Cartesian sampling. However, BiLMDM reconstructions
  have been artifact free, without any temporal bleeding irrespective of what
  sampling strategy was used to acquire data. Similar results to the ones
  discussed here hold also for the sampling strategies not included in the
  manuscript, for all employed algorithms. Even in the case of the
  experimentally acquired data, where techniques like PS-Sparse, MLS, XD-GRASP
  and LASSI suffered from some missing structures and distortions, while SToRM
  suffered temporal blurring, BiLMDM produced images which seem to best preserve
  the features of the cardiac structure spatially and temporally. BiLMDM seems
  capable of reconstructing desirable MR images irrespective of the sampling
  strategy into consideration, which appears to be a dominant factor in the
  competing state-of-the-art methods.





§ CONCLUSIONS


This paper proposed the novel bi-linear modeling for data manifolds (BiLMDM); a
new framework for data reconstruction using manifold-learning and
sparse-approximation arguments. BiLMDM comprises several modules: Extracting a
set of landmark points from a data cloud helps in learning the latent manifold
geometry while identifying low-dimensional renditions of the landmark points
facilitates efficient means for data storage and computations. Finally, a
bi-linear optimization task is used to achieve data recovery. Quantitative and
qualitative analyses on dynamic MRI data, described in Sec. <ref>,
provided evidence that the proposed BiLMDM achieves improvements in dMRI image
reconstruction and artifact suppression over state-of-the-art approaches such as
PS-Sparse, MLS, LASSI, XD-GRASP and SToRM in cardiac cine
data. BiLMDM reconstruction is not limited to cine MR data but has
  also produced reliable reconstructions for perfusion MR data. The efficacy of
  BiLMDM on perfusion data and extensive comparison tests versus
  state-of-the-art schemes are reserved for a future publication. By
introducing BiLMDM, this work paves the way for devising efficient ways to
utilize fewer data points for reconstruction than state-of-the-art solutions and
opens the door for further advances in geometric data-approximation methods.







§ MATHEMATICAL PRELIMINARIES



For positive integers m, n, the space of matrices ^m× n is equipped with the
inner product AB := (A^B),
∀A, B∈^m× n, where the superscript  denotes the
Hermitian transpose of a matrix. It is worth noting that the inner product is not commutative:
BA = AB, where the overline symbol
denotes the complex conjugate of a number. The induced norm of ^m× n by the
previous inner product coincides with the Frobenius norm of a matrix:
·_F = ··^1/2. Moreover, the spectral norm
A_2 of matrix A∈^m× n is defined as
λ_max^1/2(A^A), where λ_max(·) stands for the
maximum eigenvalue of a symmetric matrix. In the case where m = n, then
A_2 = λ_max(A).

Given the positive integers m_1, m_2, n_1, n_2 and the linear mapping
ℒ: ^m_1× n_1→^m_2× n_2, the adjoint of ℒ is
the linear mapping ℒ^*: ^m_2× n_2→^m_1× n_1 defined as
Aℒ(B) = ℒ^*(A)B,
∀A∈^m_2× n_2, ∀B∈^m_1× n_1. For
example, with regards to the sampling mapping 𝒮(·) in Sec. <ref>, its
adjoint 𝒮^*(·) = 𝒮(·), , 𝒮(·) is self-adjoint, and
𝒮^2(·) = 𝒮(·). Moreover, for the MATLAB implementation of the Fourier
transform <cit.>, ℱ^* = N_kℱ^-1 and
ℱ_t^* = N_frℱ_t^-1. If matrix A∈^m× n
is viewed as a linear mapping A: ^n→^m, then
A^* = A^. Mapping ℒ(·) is called Lipschitz continuous,
with coefficient L>0, if
ℒ(A) - ℒ(B)_F≤ L A -
  B_F, ∀A, B.

Given a convex function g(·): ^m× n→ and a positive real number
λ, the proximal mapping
_λ g(·): ^m× n→^m× n is defined as
_λ g(A) := min_Bλ g(B) + (1/2)
A-B_F^2. For example, in the case where g becomes the indicator
function ι_𝒞 with respect to a closed convex set
𝒞⊂^m× n, , ι_𝒞(A) := 0, if
A∈𝒞, while ι_𝒞(A) := +∞, if
A∉𝒞, then _λι_𝒞 becomes the metric
projection mapping onto 𝒞:
_λι_𝒞(A) = min_B∈𝒞A-B_F. Moreover, in the case where g is the ℓ_1-norm
·_1, then the (i,j)th entry of _λ·_1(A) is given by
the soft-thresholding rule <cit.>

    [_λ·_1(A)]_ij = [A]_ij(1 -
      λmax{λ, |[A]_ij|})
       .





§ SOLVING FOR Û_N


This refers to the convex minimization sub-task (<ref>) and provides important
details essential to the implementation of Alg. <ref>. All derivations,
including those in Sec. <ref>, are performed on the basis of viewing the complex-valued
U as ((U), (U)), where (·) and (·) stand for the
real and imaginary parts, respectively, of a complex-valued matrix. Gradients are
not considered in the complex-differentiability sense (Cauchy-Riemann
conditions) <cit.>. Nevertheless, to save space and use compact mathematical
expressions, all subsequent results, including those in Sec. <ref>, are stated in their
complex-valued form.

Upon defining the convex constraint
𝒞_i := U Ue_i≤ C_U , ∀ i∈1, …,
  d, (<ref>) can be expressed as

    Û_n ∈min_U g_1(U) +
      g_2(U) ,

where

  g_1(U) 
       :=         12 𝒮(Y)    -
    𝒮ℱ(UΛ̌B_n)^2_F
    +  τ_U2U - U_n^2_F 

            + λ_12 Z_n-
      ℱ_t(UΛ̌B_n)_F^2 ,  

  g_2(U) 
      :=        ∑_i=1^d ι_𝒞_i(U) . 


It can be verified that ∀U,

  ∇g_1(U) 
       =         [ N_k ℱ^-1𝒮 ℱ 
    (UΛ̌B_n) + λ_1 N_fr
    UΛ̌B_n ] B_n^ Λ̌^
  

            + τ_U (U - U_n) 

            - [ N_k ℱ^-1𝒮(Y) + λ_1
  N_frℱ_t^-1(Z_n) ] B_n^
  Λ̌^ . 

By virtue of the fact ℱ^-1𝒮ℱ (U)_F≤U_F, ∀U, it can be also verified that ∀U_1,
U_2,

    ∇g_1(U_1) - ∇g_1(U_2)_F
       ≤*(N_k + λ_1 N_fr)
        Λ̌B_n B_n^Λ̌^ +
        τ_UI_d_2 ·U_1 - U_2_F ,

which yields the Lipschitz coefficient

    L := (N_k + λ_1 N_fr) ·λ_max( Λ̌B_n B_n^Λ̌^)
      + τ_U  .


The ith column of _λ g_2(U), ∀λ>0 and ∀U,
is computed by the (metric) projection mapping onto 𝒞_i:

    _λ g_2(U) e_i = C_Umax{C_U, Ue_i}Ue_i  .

Lastly, since (<ref>) is not affinely constrained, <cit.> suggests

    T(U) = U,   ∀U ,

for Alg. <ref>.



§ SOLVING FOR B̂_N


This refers to the convex minimization sub-task (<ref>) and provides important
details, essential to the implementation of Alg. <ref>. Upon defining the affine
constraint
𝒞_aff := B|1_N_ℓ^⊺B = 1_N_fr^⊺, (<ref>)
can be expressed as

    B̂_n ∈min_B∈𝒞_aff g_1(B) +
      g_2(B) ,

where

  g_1(B)
       =        
  12 𝒮(Y) -
    𝒮ℱ(U_nΛ̌B)^2_F 
  + τ_B2B - B_n^2_F
 
            + λ_12 Z_n -
      ℱ_t(U_nΛ̌B)_F^2  , 
 
  g_2(B)
      =        λ_3B_1  . 


It can be verified that ∀B,

  ∇g_1(B)
       =         Λ̌^U_n^
  [N_k ℱ^-1𝒮
  ℱ(U_nΛ̌B) + λ_1 N_fr
  U_nΛ̌B] 

            + τ_B (B - B_n) 

            -Λ̌^ U_n^
  [N_k ℱ^-1𝒮(Y) 
  + λ_1 N_fr ℱ_t^-1(Z_n)] , 

and by steps similar to those in Sec. <ref>, ∇g_1(·) is Lipschitz
continuous with coefficient

    L = (N_k+ λ_1 N_fr) ·λ_max(Λ̌^U_n^U_nΛ̌) + τ_B  .

According to (<ref>), _λ g_2(·) becomes:
∀B,

    [_λ g_2]_ij = [B]_ij( 1 - λλ_3max{λλ_3,|[B]_ij|})
       .

Moreover, according to <cit.>, mapping T(·), used in
Alg. <ref>, is set equal to the projection mapping onto
𝒞_aff: ∀B,

    T(B) = B - 1N_ℓ1_N_ℓ
      (1_N_ℓ^⊺B - 1_N_fr^⊺)  .







10

url@samestyle








Liang.Lauterbur.book
Z.-P. Liang and P. C. Lauterbur, Principles of Magnetic Resonance
  Imaging: A Signal Processing Perspective.1em plus 0.5em minus
  0.4emIEEE Press, 2000.

Petersen.Middleton.62
D. P. Petersen and D. Middleton, “Sampling and reconstruction of
  wave-number-limited functions in N-dimensional Euclidean spaces,”
  Information and Control, vol. 5, no. 4, pp. 279–323, 1962.

Liang.Lauterbur.dMRI.94
Z.-P. Liang and P. C. Lauterbur, “An efficient method for dynamic magnetic
  resonance imaging,” IEEE Trans. Medical Imag., vol. 13, no. 4, pp.
  677–686, 1994.

wang2016accelerating
S. Wang, Z. Su, L. Ying, X. Peng, S. Zhu, F. Liang, D. Feng, and D. Liang,
  “Accelerating magnetic resonance imaging via deep learning,” in
  Proc. ISBI, 2016, pp. 514–517.

chen2017low
H. Chen, Y. Zhang, M. K. Kalra, F. Lin, Y. Chen, P. Liao, J. Zhou, and G. Wang,
  “Low-dose CT with a residual encoder-decoder convolutional neural
  network,” IEEE Trans. Medical Imag., vol. 36, no. 12, pp.
  2524–2535, 2017.

zhu2018image
B. Zhu, J. Z. Liu, S. F. Cauley, B. R. Rosen, and M. S. Rosen, “Image
  reconstruction by domain-transform manifold learning,” Nature, vol.
  555, no. 7697, p. 487, 2018.

jin2017deep
K. H. Jin, M. T. McCann, E. Froustey, and M. Unser, “Deep convolutional neural
  network for inverse problems in imaging,” IEEE Transactions on Image
  Processing, vol. 26, no. 9, pp. 4509–4522, 2017.

aggarwal2018model
H. K. Aggarwal, M. P. Mani, and M. Jacob, “Model based image reconstruction
  using deep learned priors (MODL),” in 2018 IEEE 15th International
  Symposium on Biomedical Imaging (ISBI 2018).1em plus 0.5em minus
  0.4emIEEE, 2018, pp. 671–674.

schlemper2018deep
J. Schlemper, J. Caballero, J. V. Hajnal, A. N. Price, and D. Rueckert, “A
  deep cascade of convolutional neural networks for dynamic MR image
  reconstruction,” IEEE Trans. Medical Imag., vol. 37, no. 2, pp.
  491–503, 2018.

hammernik2018learning
K. Hammernik, T. Klatzer, E. Kobler, M. P. Recht, D. K. Sodickson, T. Pock, and
  F. Knoll, “Learning a variational network for reconstruction of accelerated
  MRI data,” Magnetic Resonance in Medicine, vol. 79, no. 6, pp.
  3055–3071, 2018.

mardani2017recurrent
M. Mardani, E. Gong, J. Y. Cheng, J. Pauly, and L. Xing, “Recurrent generative
  adversarial neural networks for compressive imaging,” in Proc. IEEE
  CAMSAP, 2017, pp. 1–5.

mardani2017deep
M. Mardani, E. Gong, J. Y. Cheng, S. Vasanawala, G. Zaharchuk, M. Alley,
  N. Thakur, S. Han, W. Dally, J. M. Pauly et al., “Deep generative
  adversarial networks for compressed sensing automates MRI,” arXiv
  e-print, 2017, 1706.00051.

lustig2006kt
M. Lustig, J. M. Santos, D. L. Donoho, and J. M. Pauly, “k-t SPARSE: High
  frame rate dynamic MRI exploiting spatio-temporal sparsity,” in
  Proc. ISMRM, vol. 2420, 2006.

Jung.07
H. Jung, J. C. Ye, and E. Y. Kim, “Improved k-t BLAST and k-t SENSE using
  FOCUSS,” Physics in Medicine and Biology, vol. 52, no. 11, pp.
  3201–3226, 2007.

otazo2010combination
R. Otazo, D. Kim, L. Axel, and D. K. Sodickson, “Combination of compressed
  sensing and parallel imaging for highly accelerated first-pass cardiac
  perfusion MRI,” Magnetic Resonance in Medicine, vol. 64, no. 3, pp.
  767–776, 2010.

liang2012k
D. Liang, E. V. R. DiBella, R.-R. Chen, and L. Ying, “k-t ISD: Dynamic
  cardiac MR imaging using compressed sensing with iterative support
  detection,” Magnetic Resonance in Medicine, vol. 68, no. 1, pp.
  41–53, 2012.

lingala2011ktslr
S. G. Lingala, Y. Hu, E. V. R. DiBella, and M. Jacob, “Accelerated dynamic
  MRI exploiting sparsity and low-rank structure: k-t SLR,” IEEE
  Trans. Medical Imag., vol. 30, no. 5, pp. 1042–1054, 2011.

zhao2012pssparse
B. Zhao, J. P. Haldar, A. G. Christodoulou, and Z.-P. Liang, “Image
  reconstruction from highly undersampled (k,t)-space data with joint partial
  separability and sparsity constraints,” IEEE Trans. Medical Imag.,
  vol. 31, no. 9, pp. 1809–1820, 2012.

block2007undersampled
K. T. Block, M. Uecker, and J. Frahm, “Undersampled radial MRI with multiple
  coils: Iterative image reconstruction using a total variation constraint,”
  Magnetic Resonance in Medicine, vol. 57, no. 6, pp. 1086–1098, 2007.

knoll2011second
F. Knoll, K. Bredies, T. Pock, and R. Stollberger, “Second order total
  generalized variation (TGV) for MRI,” Magnetic Resonance in
  Medicine, vol. 65, no. 2, pp. 480–491, 2011.

feng2014golden
L. Feng, R. Grimm, K. T. Block, H. Chandarana, S. Kim, J. Xu, L. Axel, D. K.
  Sodickson, and R. Otazo, “Golden-angle radial sparse parallel MRI:
  combination of compressed sensing, parallel imaging, and golden-angle radial
  sampling for fast and flexible dynamic volumetric MRI,” Magnetic
  Resonance in Medicine, vol. 72, no. 3, pp. 707–717, 2014.

poddar2016dynamic
S. Poddar and M. Jacob, “Dynamic MRI using smoothness regularization on
  manifolds (SToRM),” IEEE Trans. Medical Imag., vol. 35, no. 4, pp.
  1106–1115, 2016.

ravishankar2017.lassi
S. Ravishankar, B. E. Moore, R. R. Nadakuditi, and J. A. Fessler, “Low-rank
  and adaptive sparse signal (LASSI) models for highly accelerated dynamic
  imaging,” IEEE Trans. Medical Imag., vol. 36, no. 5, pp. 1116–1128,
  2017.

feng2016xd
L. Feng, L. Axel, H. Chandarana, K. T. Block, D. K. Sodickson, and R. Otazo,
  “XD-GRASP: Golden-angle radial MRI with reconstruction of extra
  motion-state dimensions using compressed sensing,” Magnetic resonance
  in medicine, vol. 75, no. 2, pp. 775–788, 2016.

awate2012spatiotemporal
S. P. Awate and E. V. R. DiBella, “Spatiotemporal dictionary learning for
  undersampled dynamic MRI reconstruction via joint frame-based and
  dictionary-based sparsity,” in Proc. ISBI, 2012, pp. 318–321.

wang2014compressed
Y. Wang and L. Ying, “Compressed sensing dynamic cardiac cine MRI using
  learned spatiotemporal dictionary.” IEEE Trans. Biomed. Engineering,
  vol. 61, no. 4, pp. 1109–1120, 2014.

caballero2014dictionary
J. Caballero, A. N. Price, D. Rueckert, and J. V. Hajnal, “Dictionary learning
  and time sparsity for dynamic MR data reconstruction,” IEEE Trans.  Medical Imag., vol. 33, no. 4, pp. 979–994, 2014.

Nakarmi2016accelerating
U. Nakarmi, Y. Zhou, J. Lyu, K. Slavakis, and L. Ying, “Accelerating dynamic
  magnetic resonance imaging by nonlinear sparse coding,” in Proc.  ISBI, 2016.

Wang2017parallel
Y. Wang, N. Cao, Z. Liu, and Y. Zhang, “Real-time dynamic MRI using parallel
  dictionary learning and dynamic total variation,” Neurocomputing,
  vol. 238, pp. 410–419, 2017.

nakarmi2017m
U. Nakarmi, K. Slavakis, J. Lyu, and L. Ying, “M-MRI: A manifold-based
  framework to highly accelerated dynamic magnetic resonance imaging,” in
  Proc. ISBI, 2017, pp. 19–22.

poddar2018free
S. Poddar, Y. Mohsin, D. Ansah, B. Thattaliyath, R. Ashwath, and M. Jacob,
  “Free-breathing cardiac MRI using bandlimited manifold modelling,”
  arXiv preprint arXiv:1802.08909, 2018.

poddar2018recovery
S. Poddar and M. Jacob, “Recovery of noisy points on bandlimited surfaces:
  Kernel methods re-explained,” in 2018 IEEE International Conference on
  Acoustics, Speech and Signal Processing (ICASSP).1em plus 0.5em minus
  0.4emIEEE, 2018, pp. 4024–4028.

Usman.Manifold.15
M. Usman, D. Atkinson, C. Kolbitsch, T. Schaeffter, and C. Prieto, “Manifold
  learning based ECG-free free-breathing cardiac CINE MRI,” J.  Magnetic Resonance Imag., vol. 41, no. 6, pp. 1521–1527, 2015.

nakarmi2018mls
U. Nakarmi, K. Slavakis, and L. Ying, “MLS: Joint manifold-learning and
  sparsity-aware framework for highly accelerated dynamic magnetic resonance
  imaging,” in Proc. ISBI, 2018, pp. 1213–1216.

ahmed2019free
A. H. Ahmed, Y. Mohsin, R. Zhou, Y. Yang, M. Salerno, P. Nagpal, and M. Jacob,
  “Free-breathing and ungated cardiac cine using navigator-less spiral
  SToRM,” arXiv preprint arXiv:1901.05542, 2019.

shen2017nonlinear
Y. Shen, P. A. Traganitis, and G. B. Giannakis, “Nonlinear dimensionality
  reduction on graphs,” in Proc. IEEE CAMSAP, 2017.

silva2006selecting
J. Silva, J. Marques, and J. Lemos, “Selecting landmark points for sparse
  manifold learning,” in Proc. NIPS, 2006, pp. 1241–1248.

chen2006improved
Y. Chen, M. Crawford, and J. Ghosh, “Improved nonlinear manifold learning for
  land cover classification via intelligent landmark selection,” in
  Proc. IEEE IGARSS, 2006, pp. 545–548.

Landmark.MinMax
V. De Silva and J. B. Tenenbaum, “Sparse multidimensional scaling using
  landmark points,” Stanford University, Tech. Rep., 2004.

Tu.book.08
L. W. Tu, An Introduction to Manifolds.1em plus 0.5em minus
  0.4emNew York: Springer, 2008.

RSE.13
K. Slavakis, G. B. Giannakis, and G. Leus, “Robust sparse embedding and
  reconstruction via dictionary learning,” in Proc. CISS, Baltimore:
  USA, Mar. 2013.

slavakis2017bi
K. Slavakis, G. N. Shetty, A. Bose, U. Nakarmi, and L. Ying, “Bi-linear
  modeling of manifold-data geometry for dynamic-MRI recovery,” in
  Proc. IEEE CAMSAP, 2017.

Saul.Roweis.03
L. K. Saul and S. T. Roweis, “Think globally, fit locally: Unsupervised
  learning of low dimensional manifolds,” J. Machine Learning
  Research, vol. 4, pp. 119–155, 2003.

Rockafellar.convex.analysis
R. T. Rockafellar, Convex Analysis.1em plus 0.5em minus
  0.4emPrinceton, NJ: Princeton University Press, 1970.

Elhamifar.Vidal.nips.11
E. Elhamifar and R. Vidal, “Sparse manifold clustering and embedding,” in
  Proc. NIPS, Granada: Spain, Dec. 2011.

slavakis.FMHSDM
K. Slavakis and I. Yamada, “Fejér-monotone hybrid steepest descent method
  for affinely constrained and composite convex minimization tasks,”
  Optimization, 2018, DOI: <10.1080/02331934.2018.1505885>.

friedman2001elements
J. Friedman, T. Hastie, and R. Tibshirani, The Elements of Statistical
  Learning, 2nd ed.1em plus 0.5em minus 0.4emSpringer, 2009.

facchinei2015parallel
F. Facchinei, G. Scutari, and S. Sagratella, “Parallel selective algorithms
  for nonconvex big data optimization,” IEEE Trans. Signal Process.,
  vol. 63, no. 7, pp. 1874–1889, 2015.

candes2013unbiased
E. J. Candes, C. A. Sing-Long, and J. D. Trzasko, “Unbiased risk estimates for
  singular value thresholding and spectral estimators,” IEEE Trans.  Signal Process., vol. 61, no. 19, pp. 4643–4657, 2013.

wissmann2014mrxcat
L. Wissmann, C. Santelli, W. P. Segars, and S. Kozerke, “MRXCAT: Realistic
  numerical phantoms for cardiovascular magnetic resonance,” Journal of
  Cardiovascular Magnetic Resonance, vol. 16, no. 1, p. 63, 2014.

MATLAB.2018b
“MATLAB (R2018b),” The Mathworks, Inc., Natick, Massachusetts,
  2018.

code.storm
S. Poddar and M. Jacob, “SToRM MATLAB Package,” accessed in 2018. [Online].
  Available: <https://github.com/sunrita-poddar/l2SToRM>


code.lassi
S. Ravishankar, B. E. Moore, R. R. Nadakuditi, and J. A. Fessler, “LASSI
  MATLAB Package,” accessed in 2018. [Online]. Available:
  <https://web.eecs.umich.edu/ fessler/irt/reproduce/17/ravishankar-17-lra/>


code.xdgrasp
L. Feng, L. Axel, H. Chandarana, K. T. Block, D. K. Sodickson, and R. Otazo,
  “XD-GRASP MATLAB Package,” accessed in 2019. [Online]. Available:
  <https://cai2r.net/resources/software/xd-grasp-matlab-code>


ravishankar2011mr
S. Ravishankar and Y. Bresler, “MR image reconstruction from highly
  undersampled k-space data by dictionary learning,” IEEE transactions
  on medical imaging, vol. 30, no. 5, pp. 1028–1041, 2011.

subbarao1993focusing
M. Subbarao, T.-S. Choi, and A. Nikzad, “Focusing techniques,” Optical
  Engineering, vol. 32, no. 11, pp. 2824–2837, 1993.

wang2004image
Z. Wang, A. C. Bovik, H. R. Sheikh, E. P. Simoncelli et al., “Image
  quality assessment: From error visibility to structural similarity,”
  IEEE transactions on image processing, vol. 13, no. 4, pp. 600–612,
  2004.

lin2000thin
Y.-C. Tsai, H.-D. Lin, Y.-C. Hu, C.-L. Yu, and K.-P. Lin, “Thin-plate spline
  technique for medical image deformation,” J. Medical and Biological
  Eng., vol. 20, no. 4, pp. 203–210, 2000.

maleki2013asymptotic
A. Maleki, L. Anitori, Z. Yang, and R. G. Baraniuk, “Asymptotic analysis of
  complex LASSO via complex approximate message passing (CAMP),”
  IEEE Trans. Information Theory, vol. 59, no. 7, pp. 4290–4308, 2013.

Lang.complex
S. Lang, Complex Analysis, 4th ed.1em plus 0.5em minus
  0.4emNew York: Springer, 1999.





