Multi-Contrast CSMRI Using Common
Edge Structures with LiGME Model
Daichi Kitahara, Rikako Kato, Hiroki Kuroda, and Akira Hirabayashi
College of Information Science and Engineering, Ritsumeikan University, Shiga, Japan

Abstract—CSMRI is the high-speed magnetic resonance imaging (MRI) technique using the compressed sensing (CS) theory.
Based on the fact that multiple MR images of different contrasts,e.g., T1-weighted and T2-weighted images, are scanned in clinicalpractice, Ehrhardt et al. proposed multi-contrast CSMRI utilizing
the edge information of a different contrast image obtained fromthe full-sampling k-space data. In this paper, we propose to extend
the method of Ehrhardt et al. to the linearly involved generalized
Moreau enhanced (LiGME) model. Since a directional total variation based on the edge information becomes closer to a group `0
pseudo-norm by introducing the LiGME model, we will be able toreconstruct large edges more accurately. Simulations using actual
MR images demonstrate the effectiveness of the proposed method.
Index Terms—Multi-contrast MRI, compressed sensing, directional total variation, LiGME model, convex optimization.

I. I NTRODUCTION
Magnetic resonance imaging (MRI) [1]–[3] scans the insideof the human body noninvasively using a strong magnetic field.
Compared to X-ray computed tomography (CT) imaging, MRI
has the advantage of not exposing patients to X-rays. On theother hand, the disadvantage is that the examination time is relatively long. In clinical practice, multiple types of MR images
are scanned by changing the setting of the repetition time (TR)and the echo time (TE) in a single examination. Typical types
are T1-weighted and T2-weighted images as shown in Fig. 1.
T1-weighted images are good at depicting anatomical structures, while T2-weighted images are suitable for depicting lesions. By comparing such images of different characteristics, aradiologist makes an accurate diagnosis. It takes approximately
2 to 3 minutes to acquire each type of MR image, while CT
imaging takes a few seconds. Thus, there is a need to shorten
the MRI scanning time for reduction of the burden on patients.
In MRI, the Fourier coefficients of a target slice are observed
sequentially in the k-space. Therefore, reducing the number ofthe observed Fourier coefficients shortens the scanning time. If
we apply the inverse Fourier transform, which is the simplest
MR image reconstruction method, to such undersampled coefficients, a blurred image containing artifacts is reconstructed.
To reconstruct an accurate MR image from the undersampledcoefficients, the compressed sensing (CS) theory is used [4]–[12]. CS [13]–[16] is a technique to reconstruct a target signal
from a small number of measurements under the assumptionthat the target signal is sparsely represented in an appropriate
space. Based on this theory, Lustig et al. proposed a fast MRI
method called CSMRI [4], [5]. In this method, a cost function
This work was supported by JSPS Grants-in-Aid for Early-Career Scientists
Grant Number JP19K20361 (e-mail: d-kita@media.ritsumei.ac.jp).

(a) T1-weighted image(b) T2-weighted image

Fig. 1. Structural similarity between T1-weighted and T2-weighted images.

is defined as the sum of a data fidelity term, the `1 norm of thediscrete wavelet coefficients, and the isotropic total variation.
By minimizing this convex cost function with an iterative algorithm, a clear image can be reconstructed from the undersampled coefficients whose compression rate is about 1/3.
To improve the image quality, Ehrhardt et al. focused on the
structural similarity that can be found between different typesof MR images on the same slice [10]. As we can see in Fig. 1,
although their contrasts are different, edge structures are very
similar to each other. Therefore, they proposed a multi-contrast
CSMRI method in which the full-sampling is performed only
for a single type of MR image, while the Fourier coefficientsof other types of images are undersampled so that the entire
examination time can be shortened. They defined a directional
total variation using edges of a reference image obtained fromthe full-sampling, and then reconstructed the images of different contrasts from the undersampled data by minimizing thesum of the data fidelity term and the directional total variation.
In this paper, to further improve the reconstruction accuracy,
we propose to extend the method of Ehrhardt et al. to the linearly involved generalized Moreau enhanced (LiGME) model[17]. In recent years, Selesnick has proposed the generalized
minimax concave (GMC) function as an alternative measure of
the sparsity to the `1 norm [18]. The GMC function is continuous but nonconvex, and its shape is similar to the `0 pseudonorm and the `p quasi-norm (s.t. p ∈ (0, 1)). Even though the
GMC function is nonconvex, the sum of the data fidelity term(squared errors) and the GMC function can become convex ifa matrix parameter of the GMC function is appropriately set.
Thanks to the overall convexity, we can reach the global minimum without getting stuck in local minima. Abe et al. extended the GMC function, by changing the original cost fromthe `1 norm to any symmetric convex cost, as the generalized
Moreau enhanced (GME) function [17]. In addition, they proposed the LiGME model that strictly minimizes the sum of the
fidelity term and the GME function including a linear mapping.
However, it is difficult to design a matrix parameter guaranteeing the overall convexity, and hence the LiGME model has not
yet been introduced to relatively large-scale inverse problems.
In the proposed image reconstruction, we first express the
directional total variation in [10] as a group `1 norm including
a linear mapping that is the composition of the difference matrix, a rotation matrix and a weight matrix. Second, we replace
the group `1 norm with the GME function while appropriately
setting its matrix parameter to guarantee the overall convexity.
Third, we modify a proximal splitting algorithm of the LiGME
model in [17] so that the constraint set whose projection is easy
to compute can be added. Since the directional total variation
gets closer to a group `0 pseudo-norm by the LiGME model,
we can reconstruct large edges more accurately. Simulations
using the same T1-weighted and T2-weighted images as in [10]demonstrate that the proposed method reconstructs MR images
more accurately and quickly than the original method [10].
II. P RELIMINARIESLet R, R+ , and C be the sets of all real, nonnegative real,and complex numbers, respectively. The imaginary unit is denoted by ı ∈ C. For any c ∈ C, c̄ stands for its complex conjugate. We write vectors with boldface small letters and matrices
with capital letters. We express the identity matrix of order N
as IN . The transpose and Hermitian transpose of a vector or amatrix are denoted by (·)T and (·)H . The positive definitenessand semidefiniteness of a Hermitian matrix A ∈ CN ×N are denoted by A  ON and A  ON . For p > 0 and a vector x =(x1 , x2 , . . . ,qxN )T ∈ CN , the `p (quasi-)norm1 of x is definedp PNpQ s.t. GQ = N ,as kxkp :=i=1 |xi | . For integers G and
PQ qPG−1
G,Q2
define a group `1 norm as kxk2,1 := i=1
j=0 |xi+jQ | .
A. Basic MRIAmong various MRI sequences [1], the two most basic arespin echo (SE) sequences and gradient echo (GRE) sequences.
Although MR images are generally complex-valued, the phases
greatly vary in the GRE sequences, while there are few phasefluctuations in the SE sequences. This paper uses the model ofthe ideal SE sequences in [10] for simplicity, where each MR
image can be modeled as a nonnegative real vector x ∈ RN+.
In the k-space, lattice sampling of the Fourier coefficients is1
y[k1 , k2 ] = √
N

N1
2 −1

N2
2 −1

X

X

x[n1 , n2 ]e−ı2π

k1 n1
k2 n2
N1 + N2,NN
n1 =− 21 n2 =− 22(1)where2 even numbers N1 and N2 s.t. N1 N2 = N denote the1 Note that the function k·k satisfies the condition of the norm in CN ifpp ≥ 1, and the `0 pseudo-norm of x is defined as kxk0 := limp→+0 kxkpp .
2 In the ideal SE model, y[0, 0], y[0, − N2 ], y[− N1 , 0], y[− N1 , − N2 ] ∈ R22
22
holds. Moreover, for k1 6= − N21 and k2 6= − N22 , y[k1 , k2 ] = y[−k1 , −k2 ],y[− N21 , k2 ] = y[− N21 , −k2 ], and y[k1 , − N22 ] = y[−k1 , − N22 ] also hold.

numbers of vertical and horizontal pixels of an image x, k1 =− N21 , − N21 + 1, . . . , N21 − 1, k2 = − N22 , − N22 + 1, . . . , N22 − 1,and x[n1 , n2 ] is the pixel value of x. The observation modelin (1) is the discretization of an integral [3], and the actual kspace data y[k1 , k2 ] ∈ C includes both model and observation
errors. By concatenating all the observed k-space data y[k1 , k2 ]as yfull ∈ CN , the observation model in (1) is expressed as
yfull = F x + εfull ,(2)where F ∈ CN ×N is the normalized discrete Fourier transform
matrix and a noise vector εfull ∈ CN includes the two errors.
When the energy kεfull k22 of noise is small, the MR image x
can be easily reconstructed with high accuracy by projecting
the result F H yfull of the inverse discrete Fourier transform, forthe observed data yfull , onto RN(F H yfull ).3 Here the+ as PRN+
N
N
N
projection PRN: C → R+ onto R+ can be computed for any+
input complex vector by truncating the imaginary parts and replacing negative values of the remaining real parts with zeros.
B. Compressed Sensing MRI (CSMRI)In MRI, the k-space data y[k1 , k2 ] is observed sequentially,
not simultaneously. Hence, reducing the number of the observations shortens the scanning time. We observe only M (< N )Fourier coefficients y[k1 , k2 ] and concatenate them as a vectory ∈ CM . The observation model in CSMRI is expressed as
y = Fu x + ε,(3)where Fu ∈ CM ×N is the undersampled Fourier transform matrix, ε ∈ CM is a noise vector, and M
N is the compression ratio.
If we represent the sampling points inP
the k-space as a binaryNmatrix S := (si,j ) ∈ {0, 1}M ×N s.t. j=1 si,j = 1 for all iPMand i=1 si,j ≤ 1 for all j, then Fu is expressed as Fu = SF .
Let ypad ∈ CN be the padded k-space data, where the unobserved coefficients y[k10 , k20 ] are filled with y[−k10 , −k20 ] or
0.4 Then, PRN(F H ypad ) is the simplest reconstruction result,5+but this cannot be used for diagnosis because of blurring andwrap-around artifacts. In [4], [5], to reconstruct a clear image
from the undersampled data y, Lustig et al. proposed to solve
1minimize kFu x − yk22 + λkΦxk1 + µTV(x). (4)2x∈RN+In (4), the first term is the data fidelity, i.e., squared errors fory, and the second and third terms mean the regularizations. In
the ideal SE sequences, since x is real-valued, we can use regularization terms similar to those in standard image processing
techniques. The matrix Φ ∈ RN ×N denotes a discrete wavelet
transform, and TV : RN → R+qis the isotropic total variation
PN
computed by TV(x) := n=1 ∆x2h,n + ∆x2v,n = kDxk2,N
2,1
with the two-dimensional difference matrix D := [DhT DvT ]T .6
3 In the GRE sequences, the reconstructed image is F H yfull . Moreover, in
the SE sequences, x can be reconstructed with high accuracy even from almost
half yhalf of y if we utilize the complex conjugate relations in Footnote 2.
4 This padding is called the partial Fourier method or half Fourier method.
5 In the GRE sequences, the simplest reconstruction is F H y
H
pad := Fu y.
6 Note that D ∈ {−1, 0, 1}N ×N also computes the differences betweenhthe left-end and the right-end pixel values. Similarly, Dv ∈ {−1, 0, 1}N ×N
computes the differences between the top-end and the bottom-end pixel values.

The parameters λ > 0 and µ > 0 are weights for the regularization terms. This method gives a good reconstruction result
from the observed data whose compression rate is about 1/3.
C. Multi-Contrast CSMRI Using Common Edge Structures
In clinical practice, a certain part of the body is scanned repeatedly with different settings. Thus, multiple MR images ofdifferent contrasts, e.g., T1-weighted and T2-weighted images,on the same slice are given, and a radiologist makes an accurate diagnosis by comparing them. Ehrhardt et al. focused on
this diagnostic process, and proposed to use the full-samplingin (2) for a single setting and the undersampling in (3) for othersettings, in order to shorten the overall examination time while
preserving the image quality [10]. For simplicity, here we acref
quire a T1-weighted image from the full-sampling yfull∈ CN
while a T2-weighted image from the undersampling y ∈ CM .
If noise is small, the T1-weighted image can be easily reconref
structed with high accuracy as u := PRN(F H yfull). By using+
the reconstructed T1-weighted image u as a reference, we estimate the T2-weighted image x from its undersampled data y.
Ehrhardt et al. proposed a multi-contrast CSMRI method as
1(5)minimize kFu x − yk22 + µdTV(x),2x∈RN+
where dTV : RN → R+ is the directional total variation that
evaluates the similarity of edges of x and u, and is defined by
"
#
N
X
∆xh,n
T
dTV(x) :=(I2 − ξn ξn )
∆xv,n
n=1
2
N q
X
∆x2h,n + ∆x2v,n − (2 − kξn k22 )(ξnT ∆xn )2 (6)=n=1with the horizontal and vertical normalized differences of u:
"
#
"
#
ξn,1
∆uh,n
1
ξn :=
:= q
. (7)
ξn,2
∆u2 + ∆u2 + η 2 ∆uv,n
h,n

v,n

In (7), a small positive value η > 0 is added to avoid numericalinstabilities. The vector ξn ∈ R2 is the nth edge information ofu, which is parallel to ∆un := (∆uh,n , ∆uv,n )T ∈ R2 andsatisfies 0 ≤ kξn k22 < 1. When k∆xn k22 = ∆x2h,n + ∆x2v,n is
kept constant, dTV(x) takes the minimum if each ∆xn is parallel to ξn and takes the maximum if ∆xn is orthogonal to ξn .
Hence, the T2-weighted image, having common edges with the
T1-weighted image, is reconstructed more accurately than (4).
III. M ULTI -C ONTRAST CSMRI WITH L I GME M ODEL
A. Signal Recovery Based on LiGME Model
In signal recovery based on convex optimization, we estimate a signal x ∈ RN from an observation y ∈ RM by solving
1
minimize kAx − yk22 + µΨ(Lx),
(8)
N
2
x∈R
where A ∈ RM ×N is an observation matrix, L ∈ RK×N is alinear mapping, and Ψ : RK → R is a lower semicontinuous
convex function. In standard CS-based signal recovery methods [13]–[16], the function Ψ is typically set to the `1 norm.

While CS aims to reconstruct an unknown sparse signal x
from a small number of observations y (s.t. M < N ), the `0
pseudo-norm and the `p quasi-norm (s.t. p ∈ (0, 1)) can more
precisely promote the sparsity than the `1 norm. However, theproblem in (8) becomes nonconvex when we substitute the `0
pseudo-norm or the `p quasi-norm to the regularization term Ψ,and thus finding a global minimizer is challenging.
For the case when the linear mapping L is the identity matrix IN and the function Ψ is the `1 norm, Selesnick proposed
to use the generalized minimax concave (GMC) function [18]:12ΨB (w) := Ψ(w) − min Ψ(v) + kB(w − v)k2 , (9)2v∈RK
as an alternative regularization to the `1 norm, in the form of(8), where a matrix parameter B ∈ RJ×K controls the shape of
the GMC function. The GMC function ΨB is continuous andbetter approximates the `0 pseudo-norm than the `1 norm [18].
Even though the GMC function in (9) is nonconvex, the overall
problem in the form of (8) is convex if the matrix B is appropriately set, and it is possible to find a global minimizer.
Abe et al. [17] generalized (9) for any coercive convex function Ψ satisfying Ψ(−w) = Ψ(w) for all w ∈ RK , and call it
the generalized Moreau enhanced (GME) function. Moreover,they proposed the linearly involved GME (LiGME) model:7
1(10)
minimize kAx − yk22 + µΨB (Lx).
N
2x∈R
Note that L is an arbitrary matrix in (10). Abe et al. clarified asufficient condition of B for the overall convexity of (10) as
AT A − µLT B T BL  ON ,(11)

and provided a proximal splitting algorithm that successfully
obtains a global optimal solution under the condition in (11).
B. Extension of Multi-Contrast CSMRI in (5) to LiGME ModelIn this paper, we propose to extend the optimization problemin (5) to the LiGME model for further improvement of the reconstruction accuracy. First, we show that the problem in (5)can be represented in the form of (8). When ξn 6= 0 holds forevery n, the directional total variation dTV in (6) is written by"#" ξn,1
#"
#
ξn,2
N
X
wn 0
∆xh,n
kξn k2kξn k2dTV(x) =. (12)
ξn,1
ξ∆x
0 1 − kξn,2
v,n
k
kξ
k
n=1
n 2
n 2
2

Note that (12) is a generalized form of dTV in [20] with small
positive weights wn ∈ (0, 1], and (6) is reproduced by letting
wn := 1 − kξn k22 . From (12), it can be seen that dTV decomposes the horizontal and vertical differences into components
parallel and orthogonal to ξn , and computes the `2 norm after
multiplying each parallel component by the small weight wn .
Furthermore, the condition ξn 6= 0 is not restrictive since weonly have to directly compute the `2 norm of ∆xn for n s.t.
ξn = 0. Therefore, with the horizontal and vertical difference
matrices Dh ∈ {−1, 0, 1}N ×N and Dv ∈ {−1, 0, 1}N ×N , an
1
2
7 In [19], (Ψ◦L) (x) := Ψ(Lx)−min
Bv∈RN [Ψ(Lv)+ 2 kB(x−v)k2 ],similar to ΨB (Lx) := Ψ(Lx) − minv∈RK [Ψ(v) + 12 kB(Lx − v)k22 ], isalso proposed. However, in [19], the authors consider only the denoising, i.e.,
A = IN , and Ψ ◦ L is also restricted to the one-dimensional total variation.

orthogonal matrix R ∈ R2N ×2N corresponding to rotations of∆xn , and a diagonal matrix W ∈ (0, 1]2N ×2N consisting of
wn , the directional total variation dTV can be expressed as
" #
" # !
2,N
Dh
Dh
=: Ψ W Rx
x ,dTV(x) = W R
Dv
Dv

L := W RD, Ψ := k·k2,N2,1 and g is the indicator function of
RN
,
the
image
reconstruction
algorithm9 is given as+ xk+1 = PRN [xk − 1 (AH A − µLT B H BL)xk
σ

+

µ T H

L
B Bvk − σµ LT wk + σ1 AH y]
−
σ

v
µ H
µ H
 k+1 = prox µτ Ψ [ τ B BL(2xk+1 − xk ) + vk − τ B Bvk ]which implies that the problem in (5) is represented by
" # !
Dh
1
2
minimize kFu x − yk2 + µΨ W R
x .
2
x∈RN
Dv
+

wk+1 = proxΨ∗ [L(2xk+1 − xk ) + wk ]2,1

Based on the above representation in the form of (8) and thediscussion in Section III.A, we propose to solve the following
optimization problem introducing the LiGME model:
" # !
Dh
1
2
x . (13)minimize kFu x − yk2 + µΨB W R
2
x∈RN
Dv
+
From (11), to guarantee the overall convexity of the problem
in (13), we have to design the matrix B in (9) satisfying
" #
Dh


F H S T SF − µ DhT DvT RT W B H BW R
 ON . (14)Dv
This condition is derived by substituting A := Fu = SF and
L := W RD = W R[DhT DvT ]T in (11). The matrices Dh and
Dv are block circulant with circulant blocks (BCCB) and thus
diagonalizable with the discrete Fourier transform F as Dh =b v F , where D
b h ∈ CN ×N and Dbv ∈b h F and Dv = F H D
F HD
N ×N
C
are complex diagonal matrices [6]. Although design of
B is difficult especially for a non-full row rank L [17], we candesign8 the matrix B satisfying (14), with these expressions, ass#"b † F OM ×N
θ SD
h
B=
RT W −1 ∈ C2M ×2N , (15)b v† F
2µ OM ×N S Db † ∈ CN ×N andwhere θ ∈ (0, 1], and two diagonal matrices D
h
N ×N
†
b
b v . Here B inb
are the pseudoinverses of Dh and D
Dv ∈ C(15) not only satisfies (14) but also can be computed quickly.
C. MR Image Reconstruction Algorithm
The proposed convex optimization problem in (13) is different from the original LiGME model in (10) because the nonnegative constraint x ∈ RN
+ is newly imposed. In fact, for sucha constrained LiGME model, we can also develop an iterativealgorithm that guarantees the convergence to a global optimal
solution by slightly modifying the splitting algorithm in [17].
Let g : CN → R ∪ {∞} be a convex function whose proximity operator is easily computed, and we consider a problem:
1minimize kAx − yk22 + µΨB (Lx) + g(x).
2x∈CN
From [17], [21], [22], we can prove that if ∀x Ψ(Lx) < ∞ or∀x g(x) < ∞ holds, a global optimal solution is obtained byonly adding the computation of prox σ1 g after the update of x
of the original algorithm in [17]. Since the proposed optimization problem in (13) corresponds to the case where A := Fu ,8 We can fortunately design the matrix B in (15), since not only F H F andu u
the difference matrices Dh and Dv can be diagonalized by the discrete Fouriertransform F , but also the inverse matrices of R and W are both easily given.

(16)with σ > 0 and τ > 0 which satisfy, for some κ > 0, σIN −
κ
2
κ H
T
H
2 A A − µL L  ON and τ I2N − µ( 2 + κ )B B  O2N .
IV. N UMERICAL S IMULATIONS
We show the effectiveness of the proposed method by numerical experiments using the T1-weighted and T2-weighted
images in Figs. 1(a) and 1(b), which were also used in the simulations of [10]. The number of pixels is 242×242 = 58,564 =N . We used the pixel values after division by 255 and used thesampling mask ‘cartesianX random 0 32’ in [10, Fig. 14] that
acquires M = 8,480 Fourier coefficients. Thus, the compression rate is M/N = 0.1448, approximately 1/7. Simulations
were conducted on an m-Book T510XN-M2SH5 (Windows 10,
Intel Core i7-8750H, 2.20 GHz, 32 GB) by MATLAB R2019a.
We compared the proposed method of (13) with the conventional method of (5) in terms of both reconstruction accuracy
and computation time. We exploited the code provided in [10]to reproduce the original results of the conventional method.
We first reconstructed the T2-weighted image from the observed data y by using the T1-weighed image as the reference.
An initial MR image for both the proposed and conventional
methods is generated by the inverse Fourier transform for thepadded k-space data ypad with the partial Fourier method, as
PRN(F H ypad ). Each parameter of the conventional method is+well tuned by the code of [10]. For the proposedp method, we
set the parameters to η = 0.035 in (7), wn := 1 − kξn k22 in(12), µ = 0.003 in (13), and θ = min(diag(W 2 )) in (15). The
step-size parameters σ and τ of the algorithm in (16) are set tob hH Db h + µDb vH Dbv 
σ = (κ − 1) + max diag κ2 S T S + µ D

,
T 


τ = (κ − 1) + 0.01 κ2 + κ2 max νhT νvT
where κ = 1.001, diag gives diagonal components of a squareb †H S T S Db † ) and νv = diag(Db v†H S T S D
b v† ).
matrix, νh = diag(DhhFigures 2(a), 2(b), 2(c), and 2(d) respectively show the target
image, the initial image, the result of the conventional method[10], and that of the proposed method, along with an enlarged
image of the red boxed area. We can find that the contrast of
the reconstructed image by the proposed method in Fig. 2(d) is
more clear than that by the conventional method in Fig. 2(c).
PSNR and SSIM values are also denoted in the captions. The
proposed method improved PSNR by about 1.3 dB and SSIM
by 0.02. The computation times for the proposed and conventional methods were 4.22 and 7.85 seconds, respectively. Thus,
the proposed method could reconstruct the better image about
twice as fast as the conventional method in this simulation.
9 For the non-ideal SE sequences where each MR image is a complex vectorx ∈ CN having few phase fluctuations, we only have to remove PRN in (16).
+ACKNOWLEDGMENT
We are grateful to Mr. Hiroki Okada and Mr. Yushi Sakakimoto in Kusatsu General Hospital for their valuable advice.
R EFERENCES(a) Target image(b) Initial image [22.81 / 0.6027](c) Conventional [26.78 / 0.8648](d) Proposed [28.09 / 0.8842]

Fig. 2. T1-weighted reference & T2-weighted target simulation [PSNR / SSIM].

We also reconstructed the T1-weighted image by using theT2-weighed image as the reference. We used the same values
for the parameters of the proposed method. PSNRs of the proposed and conventional methods were 24.43 dB and 23.63 dB,
and thus we confirmed that the proposed method was effective
in both CSMRI scenarios. Since T1-weighted images describe
anatomical structures better than T2-weighed images, the proposed method is suitable for the scenario in which T2-weighted
images are reconstructed with T1-weighed reference images.
V. C ONCLUSION AND F UTURE W ORKThis paper extended the multi-contrast CSMRI method of
Ehrhardt et al. to the LiGME model for improvement of thereconstruction accuracy. By the LiGME model, the directional
total variation based on the edge information of a different contrast image gets closer to a group `0 pseudo-norm, and the proposed method reconstructs edges more accurately in less computation time, which was shown by simulations. To the best of
the authors’ knowledge, this work is the first application of the
LiGME model to a large-scale underdetermined inverse problem. As future work, we plan to integrate the proposed methodwith parallel imaging techniques [23]–[25] and adapt the proposed method to MR images having large phase fluctuations.

[1] S. W. Young, Nuclear Magnetic Resonance Imaging: Basic Principles.
New York, NY, USA: Raven Press, 1984.
[2] P. G. Morris, Nuclear Magnetic Resonance Imaging in Medicine andBiology. Oxford, UK: Clarendon Press, 1986.
[3] J. P. Haldar, D. Hernando, S. K. Song, and Z. P. Liang, “Anatomicallyconstrained reconstruction from noisy data,” Magn. Reson. Med., vol. 59,
no. 4, pp. 810–818, 2008.
[4] M. Lustig, D. L. Donoho, and J. M. Pauly, “Sparse MRI: The applicationof compressed sensing for rapid MR imaging,” Magn. Reson. Imaging,vol. 58, no. 6, pp. 1182–1195, 2007.
[5] M. Lustig, D. L. Donoho, J. M. Santos, and J. M. Pauly, “Compressedsensing MRI,” IEEE Signal Process. Mag., vol. 25, pp. 72–82, 2008.
[6] Y. Chen, X. Ye, and F. Huang, “A novel method and fast algorithm for
MR image reconstruction with significantly under-sampled data,” Inv.
Prob. Imaging, vol. 4, no. 2, pp. 223–240, 2010.
[7] S. Ravishankar and Y. Bresler, “MR image reconstruction from highly
undersampled k-space data by dictionary learning,” IEEE Trans. Med.
Imaging, vol. 30, no. 5, pp. 1028–1041, 2011.
[8] C. Tang, N. Inamuro, T. Ijiri, and A. Hirabayashi, “Compressed sensing
MRI using double sparsity with additional training images,” in Proc.
ICASSP, New Orleans, LA, USA, 2017, pp. 801–805.
[9] B. Bilgic, V. K. Goyal, and E. Adalsteinsson, “Multi-contrast reconstruction with Bayesian compressed sensing,” Magn. Reson. Med., vol. 66, no.
6, pp. 1601–1615, 2011.
[10] M. J. Ehrhardt and M. M. Betcke, “Multicontrast MRI reconstruction
with structure-guided total variation,” SIAM J. Imaging Sci., vol. 9, no. 3,
pp. 1084–1106, 2016.
[11] A. Güngör, E. Kopanoğlu, T. Çukur, H. E. Güven, “Fast recovery of compressed multi-contrast magnetic resonance images,” in Proc. SPIE Med.
Imaging, Orlando, FL, USA, 2017, 6 pages.
[12] P. Song, L. Weizman, J. F. C. Mota, Y. C. Eldar, and M. R. D. Rodrigues,“Coupled dictionary learning for multi-contrast MRI reconstruction,”
IEEE Trans. Med. Imaging, vol. 39, no. 3, pp. 621–633, 2020.
[13] E. J. Candès, J. Romberg, and T. Tao, “Robust uncertainty principles:
Exact signal reconstruction from highly incomplete frequency information,” IEEE Trans. Inf. Theory, vol. 52, no. 2, pp. 489–509, 2006.
[14] D. L. Donoho, “Compressed sensing,” IEEE Trans. Inf. Theory, vol. 52,
no. 4, pp. 1289–1306, 2006.
[15] E. J. Candès, “The restricted isometry property and its implications forcompressed sensing,” Comp. Rendus Math., vol. 346, pp. 589–592, 2008.
[16] Y. C. Eldar and G. Kutyniok, Eds., Compressed Sensing: Theory andApplications. Cambridge, UK: Cambridge University Press, 2012.
[17] J. Abe, M. Yamagishi, and I. Yamada, “Linearly involved generalized
Moreau enhanced models and their proximal splitting algorithm underoverall convexity condition,” Inv. Prob., vol. 36, no. 3, 36 pages, 2020.
[18] I. Selesnick, “Sparse regularization via convex analysis,” IEEE Trans.
Signal Process., vol. 65, no. 17, pp. 4481–4494, 2017.
[19] I. Selesnick, A. Lanza, S. Morigi, and F. Sgallari, “Non-convex total
variation regularization for convex denoising of signals,” J. Math. Imaging Vis., vol. 62, no. 6–7, pp. 825–841, 2020.
[20] İ. Bayram and M. E. Kamasak, “Directional total variation,” IEEE Signal
Process. Lett., vol. 19, no. 12, pp. 781–784, 2012.
[21] L. Condat, “A primal–dual splitting method for convex optimization involving Lipschitzian, proximable and linear composite terms,” J. Optim.
Theory Appl., vol. 158, no. 2, pp. 460–479, 2013.
[22] L. Condat, D. Kitahara, A. Contreras, and A. Hirabayashi, “Proximal
splitting algorithms: A tour of recent advances, with new twists,” preprint arXiv:1912.00137, 46 pages, 2019.
[23] E. Gong, F. Huang, K. Ying, W. Wu, S. Wang, and C. Yuan, “PROMISE:
Parallel-imaging and compressed-sensing reconstruction of multicontrast
imaging using sharable information,” Magn. Reson. Med., vol. 73, no. 2,
pp. 523–535, 2015.
[24] I. Chatnuntawech, A. Martin, B. Bilgic, K. Setsompop, E. Adalsteinsson,and E. Schiavi, “Vectorial total generalized variation for accelerated
multi-channel multi-contrast MRI,” Magn. Reson. Imaging, vol. 34, no.
8, pp. 1161–1170, 2016.
[25] B. Bilgic, T. H. Kim, C. Liao, M. K. Manhard, L. L. Wald, J. P. Haldar,
and K. Setsompop, “Improving parallel imaging by jointly reconstructing multi-contrast data,” Magn. Reson. Med., vol. 80, pp. 619–632, 2018.