Also at ]Center for Advanced Materials Research, Research Institute of
Sciences and Engineering, University of Sharjah, PO Box 27272, Sharjah, United
Arab Emirates and ]Dept. of Mechanical and Materials Engineering, Florida
International University, Miami, FL33174, United States

Also at ]Nanoelectronics Integrated Systems Center (NISC), Nile University,
Cairo, Egypt and ]Dept. of Electrical and Computer Engineering, University of
Calgary, Calgary, Canada

# Tikhonov regularization for the deconvolution of capacitance from voltage-
charge response of electrochemical capacitors

Anis Allagui<EMAIL_ADDRESS>Dept. of Sustainable and Renewable Energy
Engineering, University of Sharjah, PO Box 27272, Sharjah, United Arab
Emirates  [  [  Ahmed S. Elwakil  Dept. of Electrical and Computer
Engineering, University of Sharjah, PO Box 27272, Sharjah, United Arab
Emirates  [  [

###### Abstract

The capacitance of capacitive energy storage devices can not be directly
measured, but can be estimated from the input and output signals expressed in
the time or frequency domains. Here the time-domain voltage-charge
relationship in non-ideal electrochemical capacitors is treated as an ill-
conditioned convolution integral equation where the unknown capacitance kernel
function is to be found. This comes from assuming _a priori_ that in the
frequency domain the charge is equal to the product of capacitance by voltage.
The computation of a stable solution to this problem particularly when dealing
with experimental data is highly sensitive to noise as it may lead to an
oscillating output even in the presence of small errors in the measurements.
In this work, the problem is treated using Tikhonov’s regularization method,
where a degree of damping is added to each singular value decomposition (SVD)
component of the solution, thus effectively filtering out the components
corresponding to the small singular values.

##  I Introduction

Electrochemical capacitors, also known as electric double-layer capacitors or
supercapacitors, can store large amounts of electrical energy in relatively
small volumes, and discharge it back in short times Simon, Gogotsi, and Dunn
(2014). They do so very efficiently with low charge-discharge hysteresis, with
extremely long cyclic lifetime, and with excellent safety features and cost
effectiveness, which renders them the option of choice for applications where
bursts of power are needed Shao  _et al._ (2018). The large energy storage
capability of these devices is based on the electrostatic double layer formed
at the interface between high-surface area porous electrodes in contact with a
liquid electrolyte Yang  _et al._ (2022). Specifically, the storage mechanism
in response to an external excitation involves ionic transports in opposite
directions (physical adsorption of counterions and repulsion of coions) which
can face different resistive paths depending on the local structure of the
electrodes (porosity, surface roughness, atomic-scale inhomogeneities Kerner
and Pajkossy (2000)), the dynamics of the eletric double-layer and the
formulation of the electrolytic bath Allagui  _et al._ (2020). As a result, a
frequency dispersion of capacitance, following a power-law decay of the form
$C(\omega)\propto(i\omega)^{\alpha-1}$ with $0<\alpha<1$, is commonly observed
in such devices and in other systems involving electrified
electrodes/electrolyte interfaces.

Now if the capacitance is a frequency-dependent function in the frequency
domain, it is natural that in the time domain it must be treated as a time-
dependent function. We recall that a frequency-domain representation of a
signal in terms of magnitude and phase as a function of frequency is an
indirect way of looking at time-domain data using time-to-frequency
transformations (e.g. Fourier or Laplace transforms). In the wide literature
we see that the same capacitive system is characterized on the one hand by a
constant phase element, and thus a capacitive function in the frequency
domain, and on the other hand by a constant capacitance from time-domain
measurements (see Allagui, Elwakil, and Eleuch (2021) and the references
therein). These inconsistencies have been raised and discussed in some of our
recent contributions on the topic Allagui, Elwakil, and Eleuch (2021); Allagui
and Fouda (2021); Allagui, Elwakil, and Fouda (2021); Allagui, Elwakil, and
Wang (2022), which led to renewed interest in establishing proper
characterization methods for capacitive devices Jeltsema (2023); Ortigueira
_et al._ (2023). It is needless to mention that erroneous calculation of
capacitance has direct implications on the computation of energy and power
performance of electrochemical capacitors needed for their proper integration
in larger circuit systems.

In this work we analyze the experimental time-domain charge-voltage response
of a non-ideal electric double-layer capacitor that shows a fractional
impedance behavior in the frequency domain of the constant phase element type.
The charge-voltage in the frequency domain are assumed to be related via a
multiplicative capacitive function, which is in line with the definition of
impedance function for linear systems. Thus, the corresponding input-output
time-domain problem is a convolution of a capacitive function with the applied
voltage resulting in the measured charge, from which we extract by using
Tikhonov regularization method the characteristic function of the device. This
is a continuation of our previous work on the same problem Allagui and Fouda
(2021), wherein we focused mostly on the transition between time- and
frequency-domain definitions of capacitance functions, but did not go far
enough with the proper analysis of real experimental data that are unavoidably
corrupted with noise.

Figure 1: Linear input-output convolution model of a capacitive system and the
effect of noise in the measured signal

##  II Methods

###  II.1 Theory

Electrochemical capacitors are assumed to be linear time-invariant (LTI),
causal and stable systems that map an input voltage signal $\mathbf{v[t]}$ to
an output charge $\mathbf{q[t]}$. Both input and output are real time-domain
signals. We rewrite the input voltage signal as a sequence of discrete values
such that $\mathbf{v}=\left[\bf v_{0},v_{1},\ldots,v_{N}\right]^{\text{T}}$ (T
for transpose) of finite length, and the same for the charge
$\mathbf{q}=\left[\bf
q_{0},q_{1},\ldots,q_{N}\right]^{\text{T}}=\left[f(v_{0}),f(v_{1}),\ldots,f(v_{N})\right]^{\text{T}}$.
Both voltage and charge vectors are assigned respectively to the equidistant
time instances $0^{+},t_{1},\ldots,t_{N}$ ($t_{i}=i\Delta t$, we assume
$\Delta t=1$). When a voltage $\bf v_{0}$ is applied on an capacitive device
at $t=0^{+}$, the corresponding measured charge is $\bf q_{0}=c_{0}v_{0}$,
where $\bf c_{0}$ is the value of the system impulse response at that instant.
The subsequent amount of charge $\bf q_{1}$ accumulated on the device after a
lapse of time $\Delta t$ becomes the sum of $\bf c_{0}v_{1}$ and $\bf
c_{1}v_{0}$ by considering the superposition principle in causal LTI systems.
This can be understood by realizing that the input voltage at time $\Delta t$
has now the value of $\bf v_{1}$ which, with the impulse response of the
device $\bf c_{0}$ results in a certain amount of charge $\bf c_{0}v_{1}$, but
the preceding value of voltage $\bf v_{0}$ adds also another amount of charge
$\bf c_{1}v_{0}$ that depends on the current response of the device $\bf
c_{1}$ at the instant $\Delta t$. Repeating this procedure for all discrete
times $t_{i}=i\Delta t$ until $N\Delta t$ (total time duration of the applied
signal) is the discrete convolution sum of the input voltage by the systems
impulse response which results in the output charge as:

$\bf q_{\tau}=(c*v)_{\tau}=\sum\limits_{i=0}^{N}c_{{\tau-i}}v_{i}$ (1)

A schematic illustration of this process is depicted in Fig. 1. Eq. 1 can be
expressed as a matrix multiplication representing a system of linear equations
describing the problem as:

$\bf q=v_{m}\,c$ (2)

Here $\bf v_{m}$ is a $(N+1)\times(N+1)$ lower triangular matrix of the
Toeplitz matrix, whose columns are shifted versions of the input signal vector
$\bf v$ ($\bf v_{i,j}=v_{i-j}$) Gray  _et al._ (2006), i.e.

$\bf{v_{m}}=\begin{bmatrix}\bf v_{0}&0&\dots&0\\\ \bf v_{1}&\bf
v_{0}&\dots&0\\\ \vdots&\vdots&\ddots&\vdots\\\ \bf v_{N}&\bf
v_{N-1}&\dots&\bf v_{0}\end{bmatrix}$ (3)

and the vector $\mathbf{c}=\left[\bf
c_{0},c_{1},\ldots,c_{N}\right]^{\text{T}}$ represents the unknown model
parameter values (capacitance function).

In continuous-time form, this forward problem is the convolution integral of
the input signal with the transfer function of the system, $c(t)$, which is
also referred to as the kernel of the convolution operation, such that:

$q(t)=(c\circledast v)(t):=\int_{0}^{t}c(t-\tau)v(\tau)d\tau$ (4)

The function $c(t)$ can be viewed as the system-level, macroscopic capacitance
of the device that takes into account (in a lumped form) all microscopic
events of charge transport and storage taking place inside the device in
response to an applied voltage.

The inverse problem of reconstructing the system response
$\mathbf{c}=\left[\bf c_{0},c_{1},\ldots,c_{N}\right]^{\text{T}}$ from the
convolution operation when the input and output signals are known as an ill-
conditioned deconvolution problem for system identification Mueller and
Siltanen (2012). Solutions to well-posed problems have the properties of
existence, uniqueness and stability of the solution (Hadamard criteria), but
this is not necessarily the case. When the possible sources of input signal
distortion have a linear effect on the output signal that can be lumped into
the identification of $\bf c$, one can obtain from Eq. 2 a naive solution as:

$\bf c=v_{m}^{-1}\,q$ (5)

if the matrix $\bf v_{m}$ is invertible. If some rows of $\bf v_{m}$ are close
to being linear combinations of each other, the matrix is nearly singular and
inversion cannot be done directly. Furthermore, when noise in the input and
output are uncorrelated, the problem given by Eq 1 takes instead the form:

$\bf q_{\tau}^{\text{meas}}=\sum\limits_{i=0}^{N}c_{{\tau-i}}v_{i}+{n_{\tau}}$
(6)

or in matrix notation as:

$\bf q^{\text{meas}}=v_{m}\,c+n=q+n$ (7)

to describe the measured data (Fig. 1). The vector $\bf n$ models errors
coming from all sources of noise, and can be considered as a random variable
with certain statistics: $E(\mathbf{n})=\mathbf{0}$ and
$E(\mathbf{n}\,\mathbf{n}^{\text{T}})=\mathbf{S^{2}}$, where $E(\cdot)$ is the
expectation operator and $\mathbf{S^{2}}$ is the positive definite variance
matrix. When the measurement errors are statistically uncorrelated, then the
variance matrix is diagonal, i.e.
$\mathbf{S^{2}}=\mathbf{\text{diag}}(s_{1}^{2},s_{2}^{2},\ldots,s_{N}^{2})$
where $s_{1}^{2},s_{2}^{2},\ldots,s_{N}^{2}$ are the standard deviations of
the errors, otherwise when the errors are correlated, the problem can still be
transformed to have a diagonal variance matrix by premultiplying Eq. 7 by the
inverse of the lower triangular Cholesky factor of $\mathbf{S^{2}}$ Rust and
O’Leary (2008). Even small perturbations in the signals can result in
amplified arbitrary perturbations in the solution, and in general at least one
of the conditions of existence, uniqueness and stability fails for the
obtained solution. Therefore, it is usually required to apply some form of a
computational regularization method to dampen out such instabilities and
arrive to a meaningful result.

We note that our problem is looked at from the time-domain perspective, but in
the frequency domain, it turns out to be more advantageous from a
computational point of view. Let $Q$, $C$ and $V$ be the discrete Fourier
transforms (DFT) of $\bf q$, $\bf c$ and $\bf v$ respectively, then the
discrete time series can be represented by the inverse DFTs (iDFT) as:

$\displaystyle\mathbf{v_{j}}$ $\displaystyle=$
$\displaystyle\frac{1}{N}\sum\limits_{m=0}^{N-1}{V_{m}}\exp\left[\frac{2i\pi
jm}{N}\right]$ (8) $\displaystyle\mathbf{c_{\tau-j}}$ $\displaystyle=$
$\displaystyle\frac{1}{N}\sum\limits_{n=0}^{N-1}{C_{n}}\exp\left[\frac{2i\pi(\tau-j)n}{N}\right]$
(9) $\displaystyle\mathbf{q_{\tau}}$ $\displaystyle=$
$\displaystyle\frac{1}{N}\sum\limits_{l=0}^{N-1}{C_{l}}\exp\left[\frac{2i\pi\tau
l}{N}\right]$ (10)

respectively. After a few algebraic manipulations, the original convolution
sum in the time domain given in Eq. 1 becomes a point-by-point multiplication
in the frequency domain of the DFTs of the signals, i.e.:

$Q=CV$ (11)

The spectrum or DFT of the solution to our problem is thus:

$C=\frac{Q}{V}$ (12)

which requires applying the iDFT to revert back to the time-domain solution
Allagui, Elwakil, and Wang (2022); Allagui, Elwakil, and Fouda (2021); Allagui
and Fouda (2021); Allagui, Elwakil, and Eleuch (2021). The DFT and the iDFT of
a signal can be computed efficiently by means of the fast Fourier transform
(FFT) algorithm in $O(N\ln N)$ operations Hansen (2002). However, frequency-
domain data are not always available, and usually one has to work with time-
domain data. Secondly, the solution can be unstable owing to the strong
response of high Fourier harmonics to arbitrary small distortion of data.

One way to recover $\bf c$ from Eq. 7, from which we will drop all subscripts
and superscripts for ease of notation, is to minimize the squares of the
errors, i.e.:

$r^{2}_{\min}(\mathbf{c})=\min_{\mathbf{c}}\\{(\bf q-v\,c)^{\text{T}}(\bf
q-v\,c)\\}$ (13)

The linear regression estimate Aster, Borchers, and Thurber (2018):

$\mathbf{\tilde{c}}=\bf(v^{\text{T}}v)^{-1}v^{\text{T}}q$ (14)

is the best linear unbiased estimate of the true value of $\mathbf{c}$ Rust
and O’Leary (2008). However, if the product $\bf(v^{\text{T}}v)$ is close to
singular, the problem has more solutions. The optimal solution giving the
least-squared error is found using the singular value decomposition (SVD)
method, which factorizes a matrix $\bf v$ (can be rectangular or square
matrix) as:

$\mathbf{v}=\mathbf{U\Sigma
V^{\text{T}}}=\sum\limits_{i=1}^{N}\mathbf{u}_{i}\,\sigma_{i}\,\mathbf{v}_{i}^{T}$
(15)

Here $\bf U$ and $\bf V^{\text{T}}$ are orthogonal matrices (i.e. $\bf
U^{\text{T}}U=V^{\text{T}}V=I$) with their columns ($\mathbf{u}_{i}$ and
$\mathbf{v}_{i}$ are the $i^{\textit{th}}$ columns of $\mathbf{U}$ and
$\mathbf{V}$, respectively) being the left and right singular vectors of $\bf
v$, respectively. The matrix $\bf\Sigma$ in between is a diagonal (scaling)
matrix whose diagonal elements $\sigma_{i}$ are the (nonnegative) singular
values of $\bf v$ ordered in a decreasing order in the diagonal as
$\sigma_{1}\geqslant\sigma_{2}\geqslant\ldots\geqslant\sigma_{N}\geqslant 0$
(note that some of the singular values may be zero). It can be shown that the
number of nonzero singular values is equal to the rank of $\mathbf{v}$, and
that the condition number of $\mathbf{v}$ is
$\text{cond}(\mathbf{v})=\sigma_{1}/\sigma_{N}$ Hansen (2002). The naive
solution to the problem is given by Hansen (2010):

$\mathbf{c}_{n}={\mathbf{V{\Sigma}^{-1}U^{\text{T}}q}}=\sum_{i=1}^{N}\frac{\mathbf{u}_{i}^{\text{T}}\mathbf{q}}{\sigma_{i}}\mathbf{v}_{i}$
(16)

Otherwise, if noise is involved in the measurements, one can compute an
estimate of the solution by simply truncating the sum at some $k<N$ such that:

$\mathbf{\tilde{c}}_{k}=\mathbf{V\tilde{\Sigma}^{-1}U^{\text{T}}q}=\sum_{i=1}^{k}\frac{\mathbf{u}_{i}^{\text{T}}\mathbf{q}}{\sigma_{i}}\mathbf{v}_{i}$
(17)

where $\bf\tilde{\Sigma}^{-1}$ is the inverse of $\bf{\Sigma}$ with some
elements of $\bf\tilde{\Sigma}^{-1}$ set to zero if the singular value is
below a certain threshold, i.e.
$\mathbf{\tilde{\Sigma}}=\text{diag}(\sigma_{1},\sigma_{2},\ldots,\sigma_{k},0,\ldots,0)$.
This means that some equations corrupted by numerical instabilities will be
discarded from the overall system of equations describing the device under
test. This transforms the ill-conditioned problem to a well-conditioned one,
but rank-deficient. As a result one should obtain a more robust and less
unstable estimate of the solution. The threshold or the filtering method to be
applied, which helps with the regularization of the solution such that
components with high frequencies are excluded, is usually related to the
quality of the data and the extent of the noise level. This is the so-called
truncated SVD (TSVD) technique, which is probably the simplest direct
regularization method that can be applied.

Figure 2: (a) Input voltage and accumulated charge (derived from the current
signal) time series on the Samxon EDLC, (b) sparsity pattern of the
248$\times$248 lower triangular Toeplitz matrix $\mathbf{v}$ constructed from
the voltage time series data in (a), (c) Picard plots of the singular values
$\sigma_{i}$ and the SVD-components $|u_{i}^{\text{T}}\mathbf{q}|$ and
$|u_{i}^{\text{T}}\mathbf{q}/\sigma_{i}|$ for the matrix $\mathbf{v}$ as a
function of index $i$ ($1\leqslant i\leqslant 248$), (d) plot of gradients of
$\sigma_{i}$ and $|u_{i}^{\text{T}}\mathbf{q}|$ for $1\leqslant i\leqslant 20$
(e) Tikhnov regularized capacitance function $\mathbf{\tilde{c}}_{\lambda}$
for different values of $\lambda$ and (f) their associated filter factors,
$f_{i}={\sigma_{i}^{2}}/{(\sigma_{i}^{2}+\lambda^{2})}$ in Eq. 18. (g) Plot of
Tikhnov-regularized capacitance function $\mathbf{\tilde{c}}_{\lambda}$ for
$\lambda=0.0403$

An alternative way to just setting the small singular values of $\mathbf{v}$
to zero is to instead increase them by small increments so that they lead to a
smaller contribution to the solution estimate. This is the Tikhonov
regularization method which gives an estimate to the solution as Golub,
Hansen, and O’Leary (1999); Hansen (2002):

$\mathbf{\tilde{c}}_{\lambda}=(\mathbf{v}^{\text{T}}\mathbf{v}+\lambda^{2}\mathbf{I})^{-1}\mathbf{v}^{\text{T}}\mathbf{q}=\sum\limits_{i=1}^{N}\left(\frac{\sigma_{i}^{2}}{\sigma_{i}^{2}+\lambda^{2}}\right)\frac{\mathbf{u}_{i}^{\text{T}}\mathbf{q}}{\sigma_{i}}\mathbf{v}_{i}$
(18)

being the result of the minimization problem:

$\min_{\mathbf{c}}\\{||\mathbf{q-vc}||^{2}+\lambda^{2}||\mathbf{c}||^{2}\\}$
(19)

Here $\lambda$ is a parameter to tradeoff between closeness to the true
solution and distance from noise. It is clear that as $\lambda\to 0$, the
estimate $\mathbf{\tilde{c}}_{\lambda}$ approaches the naive solution to the
problem, whereas as $\lambda\to\infty$, $\mathbf{\tilde{c}}_{\lambda}\to 0$.
In practical applications, $\lambda$ is bounded by the minimum and maximum of
the singular values, i.e. $\sigma_{N}$ and $\sigma_{1}$. The Tikhonov
solutions resemble the TSVD solutions when the truncation parameter $k$ and
the regularization parameter $\lambda$ are chosen such that
$\sigma_{k}\approx\lambda$ Hansen (2010). We will apply the Tikhonov
regularization method in what follows.

###  II.2 Experimental

A Samxon EDLC (DRL series, part No. DRL105S0TF12RR, rated 2.7 V, 1.0 F) is
taken as a test device. The electrical measurements are carried out on a
Biologic VSP-300 potentiostat equipped with an impedance spectroscopy (EIS)
module. The instrument’s control voltage is $\pm$10 V with a resolution of 1
$\mu$V on 60 mV range, and its current ranges are 500 mA to 10 nA with a
resolution of 760 fA. The EIS module frequency range is 7 MHz (3%, 3∘) down to
10 $\mu$Hz; 3 MHz (1%, 1∘).

##  III Results and discussion

Experimental analysis of electrochemical capacitors is usually carried out
based on their cyclic voltammetry response, i.e. current or accumulated charge
vs. triangular voltage input. In Fig. 2(a) we show the accumulated charge on
the device (computed from the time-integral of the current, also shown in the
figure) in response to a single linear voltage ramp from 0 to 2.5 V at the
voltage rate of 0.1 V s-1. The data points are downsampled for computational
purposes, and limited to the ones collected every 0.1 s only. Here we analyze
one half of a cycle (charging sequence), but the same can be done for the
subsequent discharging ramp. From a first glance we can see that the charge is
not perfectly proportional to the applied voltage Fouda  _et al._ (2020),
indicating that the capacitance cannot be treated simply as a constant, but
rather implicitly, as a time-dependent function.

The steps for extracting the capacitance function from the charge and voltage
vectors are as follows. First, in Fig. 2(b), we show for visualization
purposes the sparsity pattern of the matrix $\mathbf{v}$ (Eq. 3, of size
248$\times$248) constructed from the voltage time series data. Nonzero-valued
elements are colored in blue while zero-valued elements are displayed in
white. As mentioned above this is the lower triangular matrix of the Toeplitz
matrix whose columns are shifted versions of the time series voltage signal.

Second, we analyze the voltage matrix by SVD method Hansen (2007). Plots of
singular values $\sigma_{i}$ and the SVD-components
$|u_{i}^{\text{T}}\mathbf{q}|$ and $|u_{i}^{\text{T}}\mathbf{q}/\sigma_{i}|$
Hansen (1990) of the matrix $\mathbf{v}$ are shown in Fig. 2(c). The singular
values decay gradually and smoothly from $\sigma_{1}=\text{177.89}$ to
$\sigma_{248}=\text{0.0063}$ with $\text{cond}(\mathbf{v})=2.8246\times
10^{4}$. The absolute values of the SVD coefficients of the right-hand side,
i.e. $|u_{i}^{\text{T}}\mathbf{q}|$, show an initial decay for $1\leqslant
i\leqslant 20$ faster than that of $\sigma_{i}$ (gradients are plotted in Fig.
2(d)) indicating that our problem satisfies the discrete Picard condition in
this relatively narrow region Hansen (1990, 2007). However, if $\sigma_{i}$
decayed faster than $|u_{i}^{\text{T}}\mathbf{q}|$, then with an even small
perturbation in the measurements one should expect a highly fluctuating
solution dominated by the noise. We also note that for the same region of
$i\leqslant 20$ the terms $|u_{i}^{\text{T}}\mathbf{q}/\sigma_{i}|$ are
decaying, but at a slower rate than $|u_{i}^{\text{T}}\mathbf{q}|$. For
$i>20$, $|u_{i}^{\text{T}}\mathbf{q}|$ show a highly fluctuating behavior but
the overall trend is somehow leveling off at close to the noise level of the
measurements, whereas the terms $|u_{i}^{\text{T}}\mathbf{q}/\sigma_{i}|$ show
first a slight increase followed by fluctuations at around 10-3–10-2 in value.

Plots of the Tikhnov-regularized solutions $\mathbf{\tilde{c}}_{\lambda}$ for
different values of $\lambda$ (Eq. 18) logarithmically-spaced from the lowest
to the highest singular value (i.e. from 0.0063 to 177.89) are shown in Fig.
2(e). Plots of the associated filter factor terms,
$f_{i}={\sigma_{i}^{2}}/{(\sigma_{i}^{2}+\lambda^{2})}$ in Eq. 18, are shown
in Fig. 2(f); for $\sigma_{i}\gg\lambda$, $f_{i}\to 1$ and for
$\sigma_{i}\ll\lambda$, $f_{i}\to\sigma_{i}^{2}/\lambda^{2}$. The overall
trend of the capacitance function follows a power-law decay with time that is
clearly observed for small values of $\lambda$, whereas for large values of
$\lambda$, $\mathbf{\tilde{c}}_{\lambda}$ is greatly damped, and tends to
zero. Finally, in Fig. 2(g) we show the regularized solution for the
particular value of $\lambda=0.0403$. This value of $\lambda$ is obtained via
the L-curve method (not shown here) Hansen (2007), and seems to provide an
acceptable balance between perturbation and regularization errors in the
solution, knowing that other methods such as the generalized cross-validation
criterion may provide different results Jacquelin, Bennani, and Hamelin
(2003).

Figure 3: Plot of capacitance of the Samxon EDLC computed from the ratio of
time-domain (i) charge by voltage and (ii) differential charge by differential
voltage Figure 4: Plot of spectral impedance data of the Samxon EDLC device in
the Nyquist form of real vs. imaginary parts of impedance for different
voltage dc biases; inset shows plots of real and imaginary part of the
frequency-domain capacitance function given by Eq. 21

In comparison with the capacitance as commonly computed from the simple point-
by-point division of charge by voltage or the division of differential change
in charge by change in voltage, as shown in Fig. 3, it is clear that there are
discrepancies between the two resulting capacitance functions Allagui and
Fouda (2021); Allagui, Elwakil, and Fouda (2021); Allagui, Elwakil, and Eleuch
(2021). To further look into this, we collected the spectral impedance of the
device at different dc voltage biases from 0 to 2.5 V with stepping sine
excitations of 14.14 mV amplitude over the frequency range 10 mHz to 100 kHz.
Some of the data are plotted in Fig. 4 in terms of real vs. imaginary parts of
impedance, showing little effect of applied dc biases for this device. The
low-frequency capacitive branch from _ca._ 1 Hz down to 10 mHz can
satisfactorily be fitted with a constant phase impedance defined as:

$Z(s)=\frac{V(s)}{I(s)}=\frac{V(s)}{s\,Q(s)}=\frac{1}{s^{\alpha}C_{\alpha}}$
(20)

where $\alpha$ is a real coefficient that takes values between 0 and 1, and
$s=j2\pi f$. For example, for the case of 2.0 V dc bias, the fitting
parameters (using complex nonlinear least-squares data fitting) were found to
be $C_{\alpha}=0.999$ F sα-1 and $\alpha=0.974$ with $\chi^{2}=0.068$. The
effect of dc bias on these model parameters is also relatively small. From Eq.
20 a frequency-dependent, not constant, capacitance function is computed from
the ratio of charge by voltage as:

$C(s)=\frac{Q(s)}{V(s)}=C_{\alpha}s^{\alpha-1}$ (21)

given that we were dealing with convolution operations in the time domain.
Plots of the real and imaginary parts of capacitance as a function of
frequency for the three different dc biases are shown in the inset of Fig. 4.
The inverse Laplace transform (defined as $(2\pi
i)^{-1}\int_{\gamma-j\infty}^{\gamma+j\infty}F(s)e^{st}ds$) of $C(s)$ giving
the corresponding time-domain capacitance is a decaying power law function,
i.e.:

$c(t)=\frac{C_{\alpha}t^{-a}}{\Gamma(1-\alpha)}$ (22)

which is in accordance with the Tikhonov-regularized solution depicted in Fig.
2(e) for example. In Fig. 5, we replot $\mathbf{\tilde{c}}_{\lambda}$ with
$\lambda=0.0403$ along with the fitting model function $c(t)$ (Eq. 22) with
the best fitting parameters ($C_{\alpha}=0.160$ F sα-1, $\alpha=0.863$, using
nonlinear least-squares data fitting), and with ($C_{\alpha}=0.999$ F sα-1,
$\alpha=0.974$) from impedance spectroscopy analysis, noting that the model
parameters extracted from time-domain and frequency-domain data modeling
should not be necessarily the same. Besides the effect of the regularization
procedure, this explains the quantitative difference between the time-domain
capacitance functions presented in Fig. 5.

Figure 5: Plot of time-domain capacitance functions obtained by (i) Tikhnov
regularization of deconvolution problem of time-domain charge-voltage data,
and from inverse Laplace transform of frequency-domain fractional-order
impedance model

Thus, the discrepancies between the time-domain capacitances computed by
division of charge by voltage as shown in Fig. 3 vs. deconvolution as
presented in this work (Fig. 5) is due to confusion between definitions
Allagui, Elwakil, and Eleuch (2021). In this study, we assumed that the
recommended formula for computing the capacitance in the frequency domain is
$Q(s)=C(s)V(s)$ given that it is a direct extension of the very definition of
impedance. From there, the time-domain relationship between charge and voltage
should be the convolution $q(t)=(c\circledast v)(t)$, and not $q(t)=c(t)v(t)$.

Finally, we note that the same analysis shown above in Fig. 2 for the case of
charging linear voltage ramp has been repeated for the same device with a
constant charging voltage. We obtained comparable results that are not shown
here to avoid repetition and redundancy.

##  IV Conclusion

The capacitance function of a capacitive energy storage device cannot be
directly measured, but can be estimated from input and output signals. In this
study we showed how to treat the ill-posed convolution problem of capacitance-
by-voltage giving the accumulated charge on a non-ideal capacitive device that
can be corrupted with noise. We used the Tikhonov regularization method, and
verified the concordance between the regularized solution obtained from time-
domain data and the one obtained from frequency-to-time transformation of the
frequency-domain capacitance. The latter is derived from the definition of
impedance being the ratio of frequency-domain voltage by current.

## References

## References

  * Simon, Gogotsi, and Dunn (2014) P. Simon, Y. Gogotsi, and B. Dunn, “Where do batteries end and supercapacitors begin?” Science 343, 1210–1211 (2014).
  * Shao  _et al._ (2018) Y. Shao, M. F. El-Kady, J. Sun, Y. Li, Q. Zhang, M. Zhu, H. Wang, B. Dunn, and R. B. Kaner, “Design and mechanisms of asymmetric supercapacitors,” Chemical reviews 118, 9233–9280 (2018).
  * Yang  _et al._ (2022) Q. Yang, X. Jia, X. Li, Q. Yang, T. Zhang, X. Huang, Q. Zheng, C. Li, and J. Shao, “Understanding the capacitive charge in bulk porous electrodes by mathematical modeling,” Physical Review Applied 17, 044045 (2022).
  * Kerner and Pajkossy (2000) Z. Kerner and T. Pajkossy, “On the origin of capacitance dispersion of rough electrodes,” Electrochim. Acta 46, 207–211 (2000).
  * Allagui  _et al._ (2020) A. Allagui, H. Alnaqbi, A. S. Elwakil, Z. Said, A. Hachicha, C. Wang, and M. A. Abdelkareem, “Fractional-order electric double-layer capacitors with tunable low-frequency impedance phase angle and energy storage capabilities,” Appl. Phys. Lett. 116, 013902 (2020).
  * Allagui, Elwakil, and Eleuch (2021) A. Allagui, A. S. Elwakil, and H. Eleuch, “Highlighting a common confusion in the computation of capacitance of electrochemical energy storage devices,” J. Phys. Chem. C 125, 9591–9592 (2021).
  * Allagui and Fouda (2021) A. Allagui and M. E. Fouda, “Inverse problem of reconstructing the capacitance of electric double-layer capacitors,” Electrochim. Acta , 138848 (2021).
  * Allagui, Elwakil, and Fouda (2021) A. Allagui, A. S. Elwakil, and M. E. Fouda, “Revisiting the time-domain and frequency-domain definitions of capacitance,” IEEE Trans. Electron Devices 68 (2021).
  * Allagui, Elwakil, and Wang (2022) A. Allagui, A. S. Elwakil, and C. Wang, “Time-domain and frequency-domain mappings of voltage-to-charge and charge-to-voltage in capacitive devices,” IEEE Trans. Circuits Syst. II Express Briefs (2022).
  * Jeltsema (2023) D. Jeltsema, “Further comments on “origin of the curie-von schweidler law and the fractional capacitor from time-varying capacitance”[j. pow. sources 532 (2022) 231309],” J. Pow. Sources 558, 232600 (2023).
  * Ortigueira  _et al._ (2023) M. D. Ortigueira, V. Martynyuk, V. Kosenkov, and A. G. Batista, “A new look at the capacitor theory,” Fractal and Fractional 7, 86 (2023).
  * Gray  _et al._ (2006) R. M. Gray _et al._ , “Toeplitz and circulant matrices: A review,” Foundations and Trends in Communications and Information Theory 2, 155–239 (2006).
  * Mueller and Siltanen (2012) J. L. Mueller and S. Siltanen,  _Linear and nonlinear inverse problems with practical applications_ (SIAM, 2012).
  * Rust and O’Leary (2008) B. W. Rust and D. P. O’Leary, “Residual periodograms for choosing regularization parameters for ill-posed problems,” Inverse Problems 24, 034005 (2008).
  * Hansen (2002) P. C. Hansen, “Deconvolution and regularization with toeplitz matrices,” Numerical algorithms 29, 323–378 (2002).
  * Aster, Borchers, and Thurber (2018) R. C. Aster, B. Borchers, and C. H. Thurber,  _Parameter estimation and inverse problems_ (Elsevier, 2018).
  * Hansen (2010) P. C. Hansen,  _Discrete inverse problems: insight and algorithms_ (SIAM, 2010).
  * Golub, Hansen, and O’Leary (1999) G. H. Golub, P. C. Hansen, and D. P. O’Leary, “Tikhonov regularization and total least squares,” SIAM journal on matrix analysis and applications 21, 185–194 (1999).
  * Fouda  _et al._ (2020) M. E. Fouda, A. Allagui, A. S. Elwakil, S. Das, C. Psychalinos, and A. G. Radwan, “Nonlinear charge-voltage relationship in constant phase element,” AEU Int. J. Electron. Commun. 117 (2020).
  * Hansen (2007) P. C. Hansen, “Regularization tools version 4.0 for matlab 7.3,” Numerical algorithms 46, 189–194 (2007).
  * Hansen (1990) P. C. Hansen, “The discrete picard condition for discrete ill-posed problems,” BIT Numerical Mathematics 30, 658–672 (1990).
  * Jacquelin, Bennani, and Hamelin (2003) E. Jacquelin, A. Bennani, and P. Hamelin, “Force reconstruction: analysis and regularization of a deconvolution problem,” Journal of sound and vibration 265, 81–107 (2003).


# Magnetic Resonance Parameter Mapping using Self-supervised Deep Learning
with Model Reinforcement

Wanyu Bian  Albert Jang  and Fang Liu This work was supported in part by the
grants R21EB031185; R01AR081344; R01AR079442.” W. Bian, A. Jang and F. Liu are
with the Athinoula A. Martinos Center for Biomedical Imaging, Massachusetts
General Hospital and Harvard Medical School, Charlestown, MA 02129 USA.
(corresponding author: F. Liu<EMAIL_ADDRESS>
###### Abstract

This paper proposes a novel self-supervised learning method, RELAX-MORE, for
quantitative MRI (qMRI) reconstruction. The proposed method uses an
optimization algorithm to unroll a model-based qMRI reconstruction into a deep
learning framework, enabling the generation of highly accurate and robust MR
parameter maps at imaging acceleration. Unlike conventional deep learning
methods requiring a large amount of training data, RELAX-MORE is a subject-
specific method that can be trained on single-subject data through self-
supervised learning, making it accessible and practically applicable to many
qMRI studies. Using the quantitative $T_{1}$ mapping as an example at
different brain, knee and phantom experiments, the proposed method
demonstrates excellent performance in reconstructing MR parameters, correcting
imaging artifacts, removing noises, and recovering image features at imperfect
imaging conditions. Compared with other state-of-the-art conventional and deep
learning methods, RELAX-MORE significantly improves efficiency, accuracy,
robustness, and generalizability for rapid MR parameter mapping. This work
demonstrates the feasibility of a new self-supervised learning method for
rapid MR parameter mapping, with great potential to enhance the clinical
translation of qMRI.

Quantitative MRI, Self-supervised learning, Model reinforcement, Optimization

##  1 Introduction

Deep learning has had a profound impact on medical imaging, including MRI. The
breadth of its impact includes disease diagnosis, prognosis, image analysis,
and processing [1, 2]. In particular, deep learning methods for MRI
reconstruction are becoming increasingly popular due to their capability to
learn important image features directly from large datasets [3]. One popular
method is estimating unacquired k-space data through training end-to-end
convolutional neural networks (CNNs), using undersampled data as input and
corresponding fully sampled data as a reference [4, 5, 6, 7, 8, 9, 10]. This
enables the networks to learn the mapping function between the undersampled
and fully sampled k-space data pairs. A well-trained network can then
reconstruct undersampled k-space data for accelerated MRI.

Despite recent advances in deep learning MRI reconstruction [4, 5, 6, 7, 8, 9,
3, 10, 11, 12], several challenges remain, including limited training data
access. Acquiring fully sampled k-space data can be time-consuming and
expensive, and the data collection challenge for quantitative MRI is even more
pronounced. Quantitative MRI is an advanced method to quantify tissue MR
parameters by modeling the acquired MR signal. For example, a popular qMRI
method called variable Flip Angle (vFA) [13] to quantify tissue spin-lattice
relaxation time ($T_{1}$) requires acquisition for multiple repeated scans at
different flip angles. Since each acquisition can take several minutes, the
repeated scan time can be non-trivial. There has been some recent progress
using deep learning-based methods to accelerate qMRI [14, 15, 16, 17, 18, 19,
20, 21, 22, 23], including MANTIS by Liu et al. [18, 19], DeepDTI by Tian et
al. [20], MoDL-QSM by Feng et al. [21], and DOPAMINE by Jun et al. [22], all
of which use supervised learning to enable rapid MR parameter mapping from
undersampled k-space data.

Unsupervised or self-supervised learning methods have recently gained
increasing attention due to their nature of reduced data demand. In
unsupervised learning methods, the algorithm can identify patterns and
structures of the input data without explicit reference pair. Self-supervised
learning methods involve training a model to learn the properties and features
of the input data through self-instruction. Liu et al. recently [23] proposed
a self-supervised learning framework that incorporates an end-to-end CNN
mapping and an MRI physics model to guide the generation of MR parameter maps.
This method, referred to as REference-free LAtent map eXtraction (RELAX), has
shown excellent performance in reconstructing undersampled qMRI data for
$T_{1}$ and spin-spin relaxation time ($T_{2}$) mapping in the brain and knee.
Notably, RELAX only trains on undersampled k-space datasets without the need
for fully sampled k-space for reference, greatly alleviating the data
limitation. The end-to-end CNN mapping in RELAX is critical in converting the
undersampled k-space data into MR parameter maps through cross-domain
learning, despite the end-to-end mapping (e.g., deep U-Net) requires high
computing resources such as GPU memory, needs diverse training datasets to
provide sufficient image features and lacks convergence guarantees due to its
nonlinear network operations and deep structures.

Optimization algorithms can be used to unroll deep neural networks and have
been applied for MRI reconstruction [24, 25, 26]. Those algorithms stem from
classic mathematical methods to solve constrained optimization problems and
often have well-defined solutions. For example, Variational Network (VN) [4]
unrolls a variational model by iteratively updating along the gradient descent
direction. ADMM-net [27] learns to unroll the alternating direction method of
multipliers (ADMM). ISTA-net [28] unrolls iterative shrinkage thresholding
algorithm (ISTA) to solve the reconstruction problem with $\ell_{1}$
regularizations. PD-net [29] iterates the learnable primal-dual hybrid
gradient algorithm (PDHG) and learns the primal and dual variables in an
alternative fashion.

This paper proposes a new deep learning method for qMRI reconstruction by
marrying self-supervised learning and an optimization algorithm to unroll the
learning process. This proposed method is an extension of RELAX [23], thus
referred to as RELAX with MOdel REinforcement (RELAX-MORE), which jointly
enforces the optimization and MR physics models for rapid MR parameter
mapping. RELAX-MORE aims to enable efficient, accurate, and robust MR
parameter estimation while achieving subject-specific learning, featuring
training on single-subject data for further addressing the data limitation.

The rest of the paper is organized as follows: Section 2 introduces the theory
and a bi-level optimization scheme of RELAX-MORE. Section 3 details the
implementation of lower level and upper level optimization and describes the
experiment setup. Section 4 provides the results and discussion. Section 5
concludes the paper.

##  2 Theory

###  2.1 Supervised Learning and Self-supervised Learning for Accelerated MRI
Reconstruction

Let $\mathrm{x}\in\ \mathbb{C}^{N_{x}\times N_{y}\times N_{z}\times N_{c}}$ be
the measured MR image with a resolution $N_{x}\times N_{y}\times N_{z}$, where
$x,y,z$ are spatial dimensions and $N_{c}$ denotes the number of coil elements
in a typical multi-coil MRI reconstruction scenario. The corresponding
undersampled k-space measurement ${\mathrm{f}\in\mathbb{C}}^{N_{kx}\times
N_{ky}\times N_{kz}\times N_{c}}$ and
$\mathrm{f}=\mathbf{U}\mathcal{F}\mathbf{C}\mathrm{x}+n$ incoproates the MRI
encoding process, where $\mathbf{U}$ is an undersampling mask for accelerated
acquisition, $\mathcal{F}$ is the discrete Fourier transform, $\mathbf{C}$
denotes the coil sensitivity map and $n\in\ \mathbb{C}^{N_{kx}\times
N_{ky}\times N_{kz}\times N_{c}}$ is the signal noise during data acquisition,
and $kx,ky,kz$ are k-space dimensions.

Supervised learning uses the fully sampled image $\mathrm{x}^{\ast}\ \in\
\mathbb{C}^{N_{x}\times N_{y}\times N_{z}\times N_{c}}$ as reference and the
goal is to restore the image content from undersampled k-space data by
learning network parameters $\Theta$ through minimizing the following loss
function:

$\min_{\Theta}{\ell\left(\mathrm{x}\left(\mathrm{f}\middle|\Theta\right),\mathrm{x}^{\ast}\right)}.$
(1)

The loss function $\ell$ assesses the difference between the fully sampled
$\mathrm{x}^{\ast}$ and the network output as the reconstructed image
$\mathrm{x}\left(\mathrm{f}\middle|\Theta\right)$. Once trained with a large
amount of data pairs $(\mathrm{f},\mathrm{x}^{\ast})$, the network can learn
the mapping between undersampled measurements $\mathrm{f}$ and fully sampled
data $\mathrm{x}^{\ast}$ for further inference.

In self-supervised learning, the network is trained to learn the mapping from
undersampled k-space data $\mathrm{f}$ to itself. However, the MR physics
model for image formation must be incorporated into the learning pipeline to
guide the training process. The loss function in Equation (1) can be
reformulated as:

$\min_{\Theta}{\ell\left(\mathbf{U}\mathcal{F}\mathbf{C}\mathrm{x}\left(\mathrm{f}\middle|\Theta\right),\mathrm{f}\right)}.$
(2)

Notably, as minimizing Equation (2) involves no fully sampled data and solely
relies on learning the intrinsic structure and properties of the undersampled
data itself, making it a powerful concept for deep learning methods with
limited data [30, 31, 11].

###  2.2 Self-supervised Learning for Accelerated qMRI Reconstruction

To extend from MRI reconstruction to qMRI reconstruction, the MR parameter
maps can be denoted as $\Delta=\left\\{\delta_{i}\right\\}_{i=1}^{N}$ where
$\delta_{i}$ represents each MR parameter and $N$ is the total number of MR
parameters to be estimated. Given the MR signal model $\mathbf{M}$, which is a
function of $\Delta$, and the fully sampled image $\mathrm{x}^{\ast}$, the MR
parameters $\Delta$ can be estimated by minimizing the following problem:

$\min_{\Delta}{\ell\left(\mathbf{M}(\Delta),\mathrm{x}^{\ast}\right)}.$ (3)

This is also referred to as model-based qMRI reconstruction. Here, the MR
image $\mathrm{x}\in\mathbb{C}^{N_{x}\times N_{y}\times N_{z}\times
N_{c}\times N_{k}}$ and its undersampled k-space data $\mathrm{f}\in\
\mathbb{C}^{N_{kx}\times N_{ky}\times N_{kz}\times N_{c}\times N_{k}}$ need to
include $N_{k}$ measurements at varying imaging parameters dependent on
individual qMRI method. The Equation (3) can be formulated as:

$\min_{\Delta}{\ell\left(\mathbf{U}\mathcal{F}\mathbf{C}\mathbf{M}\left(\Delta\right),\mathrm{f}\right)},$
(4)

which incorporates the MRI encoding operators.

The above minimization problem (4) can be solved in self-supervised learning
by minimizing the following loss function:

$\min_{\Theta}{\ell\left(\mathbf{U}\mathcal{F}\mathbf{C}\mathbf{M}\left(\Delta\left(\mathrm{f}\middle|\Theta\right)\right),\mathrm{f}\right)},$
(5)

where the quantitative MR parameter
$\Delta\left(\mathrm{f}\middle|\Theta\right)$ are parametrized from deep
neural networks with learnable network parameters. In RELAX, minimizing
Equation (5) with respect to trainable network parameters $\Theta$ is to train
a deep neural network via network backpropagation for an end-to-end mapping
network.

###  2.3 The RELAX-MORE Mechanism

The proposed RELAX-MORE extends the RELAX framework by further considering
unrolling the mapping network using an optimization algorithm. This new method
is designed to solve a bi-level optimization problem, where the mapping
network learning and quantitative MR parameter estimation can be jointly
reformulated as:

$\displaystyle\min_{\Theta}\,\,$
$\displaystyle\ell(\mathbf{U}\mathcal{F}\mathbf{C}\mathbf{M}(\Delta\left(\mathrm{f}\middle|\Theta\right)),\mathrm{f})\quad\text{s.t.}$
(6a)
$\displaystyle\Delta(\mathrm{f}|\Theta)=\operatorname*{arg\,min}_{\Delta}\phi(\Delta),$
(6b) where
$\displaystyle\phi\left(\Delta\right):=\tfrac{1}{2}\parallel\mathbf{U}\mathcal{F}\mathbf{C}\mathbf{M}\left(\Delta\right)-\mathrm{f}\parallel_{2}^{2}+\lambda\mathcal{R}\left(\Delta\right).$
(6c)

Like RELAX, the upper level (6a) is a backward training process to optimize
network parameters $\Theta$ by minimizing loss function $\ell$ via normal
network backpropagation. However, unlike RELAX, the lower level (6b) can be
implemented as a forward process to solve for the optimal solution of
$\min\limits_{\Delta}{\phi\left(\Delta\right)}$ via a gradient descent based
method. The objective function $\phi$ is a regularized least square model for
quantitative MR parameter optimization. The weighted regularization
$\lambda\mathcal{R}$ provides prior information pertinent to the desired MR
parameter maps. Hand-crafted regularizations can be employed, such as total
variation (TV) ( $\parallel\cdot\parallel_{TV}$ ) or $\ell_{1}$ norm
($\parallel\cdot\parallel_{1}$) to promote different aspects of image
features. Recent studies show that regularization $\mathcal{R}_{\Theta}$ using
CNNs can improve feature characterization, suppressing the hand-crafted priors
[32, 33]. The weight $\lambda$ can be learned and integrated into
$\mathcal{R}_{\Theta}$. The objective function $\phi_{\Theta}$ depends on the
network parameters $\Theta$ that learned from regularizers, which can be
parametrized as the summation of CNNs with $\ell_{2,1}$ norm:

$\mathcal{R}_{\Theta}(\Delta):=\sum^{N}_{i=1}\|\mathcal{D}_{\Theta_{i}}(\delta_{i})\|_{2,1}.$
(7)

Each ${\mathcal{D}_{\Theta}}_{i}$ learns to extract the prior features of
different MR parameters $\delta_{i}$, for each $i=1,\ldots,N$. Therefore, the
RELAX-MORE for qMRI reconstruction can be mathematically described as solving
a CNN regularized bi-level optimization problem:

$\displaystyle\min_{\Theta}\,\,$
$\displaystyle\ell(\mathbf{U}\mathcal{F}\mathbf{C}\mathbf{M}(\Delta\left(\mathrm{f}\middle|\Theta\right)),\mathrm{f})\quad\text{s.t.}$
(8a)
$\displaystyle\Delta(\mathrm{f}|\Theta)=\operatorname*{arg\,min}_{\Delta}\phi_{\Theta}(\Delta),$
$\displaystyle=$
$\displaystyle\operatorname*{arg\,min}_{\Delta}\mathcal{H}(\Delta)+\mathcal{R}_{\Theta}(\Delta),$
(8b)

where data fidelity term
$\mathcal{H}(\Delta)=\tfrac{1}{2}\parallel\mathbf{U}\mathcal{F}\mathbf{C}\mathbf{M}\left(\Delta\right)-\mathrm{f}\parallel_{2}^{2}$.

##  3 Methods

###  3.1 The Lower Level Optimization: Forward Unrolling of A Learnable
Descent Algorithm

Proximal gradient descent, a widely used iterative algorithm, is implemented
to unroll the lower level optimization (8b) in RELAX-MORE. For each unrolled
iterative phase $t=1,\ldots,T$, the following steps is used to update MR
parameter $\delta_{i}$ for each $i=1,\ldots,N$ :

$\displaystyle\bar{\delta}_{i}^{(t)}$
$\displaystyle=\delta_{i}^{(t-1)}-\alpha^{(t)}_{i}\nabla\tfrac{1}{2}\parallel\mathbf{U}\mathcal{F}\mathbf{C}\mathbf{M}(\\{{\delta}_{i}^{(t-1)}\\}_{i=1}^{N})-\mathrm{f}\parallel_{2}^{2}$
(9a) $\displaystyle\delta_{i}^{(t)}$ $\displaystyle={\rm
prox}_{\mathcal{R}_{\Theta},\rho_{t}}({\bar{\delta}}_{i}^{(t)}),$ (9b)

where the proximal operator for regularization $\lambda\mathcal{R}(z)$ is
defined as:

${\rm
prox}_{\lambda\mathcal{R},\rho}(y):=\operatorname*{arg\,min}_{z}\tfrac{\rho}{2}\parallel
z-y\parallel+\lambda\mathcal{R}(z).$ (10)

This step involves implementing the proximal operation for regularization
$\mathcal{R}$, which is equivalent to finding the maximum-a-posteriori
solution for the Gaussian denoising problem at a noise level
$\sqrt{\lambda/\rho}$ [34, 35], thus the proximal operator can be interpreted
as a Gaussian denoiser. However, because the proximal operator ${\rm
prox}_{\mathcal{R}_{\Theta},\rho_{t}}$ in the objective function (8b) does not
admit closed form solution, a CNN is used to substitute ${\rm
prox}_{\mathcal{R}_{\Theta},\rho_{t}}$, where the network is constructed as a
residual learning network $\mathcal{G}_{\Theta}$. Network
$\mathcal{G}_{\Theta}$ is composed with CNN $\mathcal{D}_{\Theta}$, which
reflects Equation (7), and its adjoint $\widetilde{\mathcal{D}}_{\Theta}$ with
symmetric but reversed architecture to increase the network capability, and a
soft shrinkage operator $\mathcal{S}_{\beta}$ is applied in between. All the
substitutions are mathematically proven to be effective by [28, 36]. Then step
(9b) becomes

$\displaystyle\delta_{i}^{(t)}$
$\displaystyle=\mathcal{G}_{\Theta_{i}}(\bar{\delta}_{i}^{(t)})+\bar{\delta}_{i}^{(t)},$
(11a)
$\displaystyle=\widetilde{\mathcal{D}}_{\Theta_{i}}\circ\mathcal{S}_{\beta^{(t)}_{i}}\circ\mathcal{D}_{\Theta_{i}}(\bar{\delta}_{i}^{(t)})+\bar{\delta}_{i}^{(t)},\,\,\forall
i=1,\ldots,N,$ (11b)

where the soft thresholding operator is $\mathcal{S}_{\beta^{(t)}}={\rm
prox}_{\parallel\cdot\parallel_{2,1},\beta^{(t)}}(y)=[{\rm
sign}(y_{i})\max(|y_{i}|-\beta^{(t)},0)]$ for vector
$y=(y_{1},\ldots,y_{n})\in\mathbb{R}^{n}$ and $\beta^{(t)}$ is a soft
thresholding parameter that is updated at each phase $t=1,\ldots,T$. In
summary, the forward learnable proximal gradient descent algorithm to achieve
unrolled lower level optimization in RELAX-MORE can be detailed in Algorithm:

1 Input: $\delta_{i}^{(0)},\alpha_{i}^{(0)},\beta_{i}^{(0)},\,i=1,\cdots,N$

2 for  _$t=1$ to $T$_ do

3 for  _$i=1$ to $N$_ do

4 $\bar{\delta}_{i}^{(t)}=\delta_{i}^{(t-1)}-$
$\alpha_{i}^{(t)}\nabla\frac{1}{2}\|\mathbf{U}\mathcal{F}\mathbf{C}\mathbf{M}(\\{\delta_{i}^{(t-1)}\\}_{i=1}^{N})-\mathrm{f}\|_{2}^{2},$

5
$\delta_{i}^{(t)}=\widetilde{\mathcal{D}}_{\Theta_{i}}\circ\mathcal{S}_{\beta^{(t)}_{i}}\circ\mathcal{D}_{\Theta_{i}}(\bar{\delta}_{i}^{(t)})+\bar{\delta}_{i}^{(t)}$

6

7 $\mathrm{x}^{(t)}=\mathbf{M}(\\{\delta_{i}^{(t)}\\}_{i=1}^{N}),\,\,$

8

output $\\{\delta_{i}^{(T)}\\}_{i=1}^{N}$ and $\mathrm{x}^{(t)},\forall
t\in\\{1,...,T\\}$.

Algorithm 1 Forward Learnable Descent Algorithm for solving (8b)

The soft thresholding parameter $\beta^{(t)}$ and step size $\alpha^{(t)}$ are
learnable and updating during each phase. The final outputs of the algorithm
are estimated intermediate MR images $\mathbf{x}^{(t)},\forall t=1,\ldots,T$
for all iterative phases and the reconstructed MR parameter maps
$\delta_{i}^{(T)},\forall i=1,\ldots,N$ of the last phase.

###  3.2 Initial Input of the Algorithm

Using an initial MR parameter input $\delta_{i}^{(0)}$ that is closer to the
optimal solution of the Algorithm can lead to better results and eliminate the
need for excessive iterations in gradient descent. To achieve a favorable
input, initialization networks are used to learn $\delta_{i}^{(0)},\forall
i=1,…,N$ from undersampled k-space data. First, the undersampled k-space data
is converted into the image domain, followed by a Coil Combination Network to
combine all coil-wise images. $N_{k}$ coil-combined images are produced to
represent $N_{k}$ measurements during qMRI acquisition. Next, all images are
concatenated together as input into the $\delta$-initialization Network to
obtain $\delta_{i}^{(0)},\forall i=1,…,N$. The initial MR images can also be
obtained from the signal model using
$\mathrm{x}^{(0)}=\mathbf{M}(\\{\delta_{i}^{(0)}\\}_{i=1}^{N})$, which will be
subsequently used in the loss function $\ell$ for training the overall network
framework in upper level optimization in Equation (8a).

###  3.3 The Upper Level Optimization: Backward Training of Overall Network

Figure 1: Schematic framework for implementing RELAX-MORE.

The loss function (8a) for RELAX-MORE is a summation of the weighted mean
squared error (MSE) between undersampled k-space measurement $\mathrm{f}$ and
the k-space representation of estimated MR images $\mathrm{x}^{(t)},\forall
t=0,\ldots,T$ for all unrolled phases. Therefore, the overall network training
objective in the upper level optimization is to minimize the loss function
$\ell$ explicitly expressed as

$\ell(\mathbf{U}\mathcal{F}\mathbf{C}\mathrm{x}^{(t)},\mathrm{f})=\sum^{T}_{t=0}\gamma_{t}\parallel\mathbf{U}\mathcal{F}\mathbf{C}\mathrm{x}^{(t)}-\mathrm{f}\parallel^{2}_{2},$
(12)

where the $\gamma_{t}$ are the weighting factors for the MSE in order of the
importance of the output from each phase and the initial $\mathrm{x}^{(0)}$.

The detailed framework of the proposed RELAX-MORE is illustrated in Fig. 1 The
Coil combination Network and $\delta$-initialization Network applied 4 repeats
of Conv($3\times 3$ kernel size, 32 kernels) followed by a ReLU activation
function. The last convolution layer has a single kernel Conv($3\times 3$, 1).
The CNN $\mathcal{D}$, applied 8 Conv($3\times 3$, 64) with a ReLU and the
last convolution also has a single kernel Conv($3\times 3$, 1). The CNN
$\widetilde{\mathcal{D}}$ applied a reversed network structure.

###  3.4 Experiment Setup

To investigate the feasibility of RELAX-MORE method for reconstructing
accelerated qMRI, we used the widely studied $T_{1}$ mapping through the vFA
method [13] as an example. Herein, the MR signal model $\mathbf{M}$ can be
expressed as:

$\mathbf{M}_{k}(\mathbf{T}_{1},\mathbf{I}_{0})=\mathbf{I}_{0}\cdot\frac{(1-e^{-TR/\mathbf{T}_{1}})\sin{\eta_{k}}}{1-e^{-TR/\mathbf{T}_{1}}\cos{\eta_{k}}},$
(13)

where several MR images are acquired at multiple flip angles $\eta_{k}$ for
$k=1,\ldots,N_{k}$. $\mathbf{T}_{1}\in\mathbb{R}^{N_{x}\times N_{y}\times
N_{z}},\mathbf{I}_{0}\in\mathbb{C}^{N_{x}\times N_{y}\times N_{z}}$ are the
spin-lattice relaxation time map and proton density map, respectively,
reflecting the imaged tissue relaxation properties, with $T_{1}$ value
sensitive to many brain and knee joint diseases. The set of MR parameters
estimated in this model are $\Delta=\\{\mathbf{T}_{1},\mathbf{I}_{0}\\}$. The
flip angles and other imaging parameters, such as repetition time TR are pre-
determined.

The experiments include in-vivo studies on the brain and knee of healthy
volunteers and ex-vivo phantom studies, all of which were carried out on a
Siemens 3T Prisma scanner. For the brain study, the purpose is to investigate
the efficiency and performance of RELAX-MORE and compare it with other state-
of-the-art reconstruction methods. The vFA on the brain of five subjects was
performed in the sagittal plane using a spoiled gradient echo sequence at
imaging parameters TE/TR = $12/40$ ms, FA =
$5^{\circ},10^{\circ},20^{\circ},40^{\circ}$, FOV = $230\times 230\times 160$
mm and matrix size = $176\times 176\times 48$ with a dedicated 20-channel head
coil. For the knee study, the purpose is to investigate the robustness of
RELAX-MORE against severe noise contamination and image imperfection.
Therefore, the vFA on one knee was performed using a 4-channel receiving-only
flex coil at parameters TE/TR = $12/70$ ms, FA =
$10^{\circ},20^{\circ},40^{\circ}$, FOV = $160\times 137\times 108$ mm and
matrix size $224\times 192\times 36$. For the phantom study, the purpose is to
investigate the reconstruction accuracy of RELAX-MORE. The vFA phantom data
was acquired along the coronal plane with sequence parameters TE/TR = $12/80$
ms, FA = $5^{\circ},10^{\circ},20^{\circ},40^{\circ},60^{\circ}$, FOV =
$170\times 170\times 60$ mm and matrix size $128\times 128\times 12$ using the
20-channel head coil. In all 3 studies, to minimize bias on $T_{1}$ estimation
due to $B_{1}^{+}$ inhomogeneity and imperfect spoiling, $B_{1}^{+}$ maps were
separately acquired for compensation [37], and $169^{\circ}$ linear RF phase
increments between subsequent repetition cycles and strong gradient spoilers
were applied to minimize the impact of imperfect spoiling [38].

Figure 2: Exemplified masks: (a) 3x 1D Cartesian variable density
undersampling mask and (b) 4x 2D Poisson disk undersampling mask.

The fully acquired MRI k-space data were undersampled retrospectively using
two undersampling schemes: 1) 1D Cartesian variable density undersampling at
acceleration factor AF = $3\times$ with the 16 central k-space lines fully
sampled (Fig. 2(a)) and 2) 2D Poisson disk undersampling at AF = $4\times$
with the central $51\times 51$ k-space portion fully sampled (Fig. 2(b)). The
undersampling patterns were varied for each flip angle, like in previous
studies [18, 23]. The coil sensitivity maps were estimated from ESPIRiT [39].

The selection of network hyperparameters was non-trivial and optimized by
empirically searching on a reasonably large range. More specifically, we used
initial $\alpha_{i}^{(0)}=0.1,\beta_{i}^{(0)}=1e^{-5}$ for brain images;
initial $\alpha_{i}^{(0)}=0.5,\beta_{i}^{(0)}=1e^{-7}$ for knee images; and
$\alpha_{i}^{(0)}=0.9,\beta_{i}^{(0)}=1e^{-3}$ for phantom images. The loss
function applied $\gamma_{t}=1$ at $t=T$ and $\gamma_{t}=1e^{-4}$ for other
phase steps. The network parameters $\Theta$ were initialized using Xavier
initialization [40]. Because RELAX-MORE trains on single-subject data, the
batch size was set to include the entire image volume. The network was trained
for $10,000$ epochs using the Adam optimizer [41]. The learning rate, which
controls the step size taken in each iteration of optimizing the loss
function, was set to $1e^{-4}$ for brain and knee images, and $5e^{-4}$ for
the phantom image. It should be noted that, unlike conventional deep learning
reconstruction methods where training and testing are two separate steps, in
RELAX-MORE with subject-specific self-supervised learning, the reconstruction
has readily concluded once the training has converged for one subject. All the
programming in this study was implemented using Python language and PyTorch
package, and experiments were conducted on one NVIDIA A100 80GB GPU and an
Intel Xeon 6338 CPU at Centos Linux system.

##  4 Results

###  4.1 Impact of Unrolling Gradient Descent Algorithm

Figure 3: $T_{1}$ maps generated from RELAX-MORE at different phases of the
forward learnable descent algorithm $(3\leq T\leq 9)$ and $T_{1}$ map obtained
from fully sampled data using the standard pixel-wise fitting. While
performance gain is significant with increasing phases at lower phase numbers,
incremental performance gain is observed after phase 7.

The reconstruction performance of RELAX-MORE is affected by the degree of
unrolling the gradient descent algorithm. Fig. 3 demonstrates the evolution of
the estimated $T_{1}$ maps from fully sampled data at different total
unrolling phase numbers from $T=3$ to $T=9$ of the gradient descent algorithm.
The larger $T$ reflects deeper unrolling operation and thus requires more
computing resources. Referring to the zoom-in views of Fig. 3, one can observe
that with increasing total phase number, the algorithm starts to better
distinguish the $T_{1}$ values from white matter (WM), grey matter (GM), and
cerebrospinal fluid (CSF) regions. Initially blurry at a lower $T$, the tissue
details start to sharpen with increasing $T$, resembling better the fully
sampled $T_{1}$ map obtained using the standard pixel-wise fitting. However,
the differences between phase $T=7,8$ and $9$ are negligible, indicating that
performance gain can reach a plateau with a sufficient depth of unrolling
steps. Therefore, to ensure consistent experiment setup and balance the trade-
off between algorithm complexity and reconstruction performance, the number of
$T$ was set to $8$ for all the experiments thereafter. This result illustrates
the effectiveness of unrolling gradient descent algorithm in the RELAX-MORE
framework.

###  4.2 Comparison with State-of-the-Art Methods

RELAX-MORE was compared with two state-of-the-art non-deep learning qMRI
reconstruction methods and one self-supervised deep learning method. These
methods include 1) Locally Low Rank (LLR) [42], where image reconstruction was
first performed, followed by pixel-wise parameter fitting. LLR exploits the
rank-deficiency of local image regions along the acquisition parameter
dimension to accelerate parameter mapping; 2) Model-TGV [43], a model-based
qMRI reconstruction method that improves the piece-wise constant region
restriction of total variation through a generalization of the total variation
theory and 3) RELAX [23], an end-to-end self-supervised deep learning method
for rapid MR parameter mapping. LLR and Model-TGV were implemented using the
recommended optimization parameters in their original papers and code.
Contrary to the original RELAX implementation, which was trained using many
subjects, in the interest of a fair comparison, our implementation of RELAX
was carried out on a single-subject training using an Adam optimizer with
$10,000$ epochs.

Figure 4: (a) Qualitative comparison among different methods using $3\times$
1D Cartesian variable density undersampling mask; (b) Zoom-in view with
details of the region outside the brain; (c) Zoom-in view with in-brain
detail; (d) pixel-wise error maps of the in-brain region details.

Figure 5: (a) Qualitative comparison among different methods using $4\times$
2D Poisson disk undersampling mask; (b) Zoom-in view with details of the
cerebellum; (c) Zoom-in view with in-brain detail; (d) pixel-wise error maps
of the in-brain region details.

$T_{1}$ maps estimated from $3\times$ 1D Cartesian variable density
undersampling using the different methods are presented in Fig. 4. As can be
observed in Fig. 4(a), the zero-filled $T_{1}$ map obtained by pixel-wise
fitting the undersampled images is noisy, displaying ripple artifacts due to
aliasing. Although LLR can partially remove these artifacts, it retains the
noisy signature. Model-TGV averages out the noise, producing a cleaner tissue
appearance and better $T_{1}$ contrast, but it is over-smoothed, resulting in
blurry maps. On the other hand, RELAX removes noises and artifacts and
produces sharpened maps but remains somewhat blocky. This is hypothesized due
to the difficulty of converging the end-to-end network using single-subject
data. RELAX-MORE produces artifact-free $T_{1}$ maps with anticipated
excellent performance for noise removal, resulting in a good performance in
appearance and contrast. This is further witnessed in the zoom-in view of Fig.
4(b), where the zero-filled, LLR, and even fully sampled images cannot
reliably estimate $T_{1}$ in the cervical spine where the signal is low due to
insufficient head coil coverage. Those pixel-wise fitting methods are
typically prone to image noise contamination. Although Model-TGV and RELAX can
better estimate $T_{1}$ in these low SNR regions, the blurry Model-TGV results
make it difficult to distinguish between the spinal cord and subarachnoid
space, whereas RELAX shows a disconnect in the subarachnoid space which may be
due to over-sharpening. However, RELAX-MORE shows remarkably good performance
in maintaining contrast and tissue details at those regions, enabling clear
distinction between the spinal cord and subarachnoid space regions. In Fig.
4(c) of the cerebrum zoom-in view, RELAX-MORE consistently exhibits the best
reconstruction performance in correcting artifacts, removing noises, and
preserving tissue contrast and details. Referring to the error maps (Fig.
4(d)), which is the absolute difference between the estimated $T_{1}$ maps in
Fig. 4(c) and the fully sampled $T_{1}$ map, overall, zero-filled shows the
most significant error followed by LLR. Model-TGV and RELAX exhibit similar
error maps, whereas RELAX-MORE produces the least error.

$T_{1}$ maps estimated from $4\times$ 2D Poisson disk undersampling are shown
in Fig. 5. Comparing the images in Fig. 5(a), the overall signature
differences are like the $3\times$ 1D undersampling case except for the
undersampling artifacts being noise-like due to 2D undersampling. Referring to
the zoom-in view and comparing the estimated $T_{1}$ maps again show similar
differences to the $3\times$ 1D undersampling case. However, compared with
other methods, RELAX-MORE can clearly distinguish between WM and GM $T_{1}$ of
the cerebellum in Fig. 5(b), particularly in the posterior part. The absolute
difference maps taken between the estimated $T_{1}$ maps in Fig. 5(c) and the
fully sampled map, shown in Fig. 5(d), exhibit similar results to the
$3\times$ 1D undersampling case with RELAX-MORE showing the least error. This
is likely achieved through integrating deep learning capability and unrolling
gradient descent algorithm with proximal operator prioritizing noise
suppression without compromising the fidelity and clarity of the underlying
tissue structure.

Further performance comparison was carried out using peak-signal-to-noise-
ratio (PSNR), structural similarity index measure (SSIM), and normalized mean
squared error (NMSE) as evaluation metrics for $T_{1}$ maps. PSNR is a
measbure of the quality of the reconstruction, while SSIM is a measure of the
similarity between two maps. PSNR, SSIM, and NMSE are defined in the
following:

$PSNR(\mathrm{v},\mathrm{v}^{*})=20\log_{10}(\max(|\mathrm{v}^{*}|)/\tfrac{1}{N}\parallel\mathrm{v}-\mathrm{v}^{*}\parallel^{2}),$
(14)
$SSIM(\mathrm{v},\mathrm{v}^{*})=\frac{(2\mu_{\mathrm{v}}\mu_{\mathrm{v}^{*}}+c_{1})(2\sigma_{\mathrm{v}\mathrm{v}^{*}}+c_{2})}{(\mu_{\mathrm{v}}^{2}+\mu_{\mathrm{v}^{*}}^{2}+c_{1})(\sigma_{\mathrm{v}}^{2}+\sigma_{\mathrm{v}^{*}}^{2}+c_{2})},$
(15)
$NMSE(\mathrm{v},\mathrm{v}^{*})=\parallel\mathrm{v}^{*}-\mathrm{v}\parallel^{2}_{2}/\parallel\mathrm{v}^{*}\parallel^{2}_{2},$
(16)

where $\mathrm{v},\mathrm{v}^{*}$ represent the estimated map and reference
map. $\mu_{\mathrm{v}},\mu_{\mathrm{v}^{*}}$ are local means of pixel
intensity, $\sigma_{\mathrm{v}},\sigma_{\mathrm{v}^{*}}$ are the standard
deviations and $\sigma_{\mathrm{v}\mathrm{v}^{*}}$ is covariance between
$\mathrm{v}$ and $\mathrm{v}^{*}$, $C_{1}=(K_{1}L)^{2},C_{2}=(K_{2}L)^{2}$ are
two constants that avoid denominator to be zero, and $K_{1}=0.01,K_{2}=0.03$.
$L$ is the largest pixel value of the magnitude of image.

$T_{1}$ maps obtained from fully sampled data were used for $\mathrm{v}^{*}$,
as a reference to compare the performance among different methods. The metric
calculations were carried out on results from $3\times$ 1D Cartesian variable
density undersampling, using the brain regions from all five subjects.

Figure 6: Quantitative comparison between the estimated $T_{1}$ maps with
$3\times$ 1D Cartesian mask and the reference $T_{1}$ maps estimated from
fully sampled data. The comparison was made based on three different measures,
namely PSNR, SSIM, and NMSE with median line and mean markers $\times$.

The analysis results are presented in Fig. 6, where the average values of
PSNR, SSIM and NMSE along with corresponding standard deviations among five
subjects are presented using box-whisker plots. For both PSNR and SSIM, RELAX-
MORE produces the highest mean value, followed by RELAX, Model-TGV and LLR.
Both PSNR and SSIM of RELAX-MORE show the highest consistency, as evidenced by
its small standard deviation compared to other methods. Zero-filled shows the
lowest mean values for PSNR and SSIM with the largest standard deviation,
which can be attributed to aliasing artifacts. RELAX-MORE shows the lowest
mean NMSE, followed by RELAX, Model-TGV and LLR. Zero-filled shows the highest
mean NMSE. Overall, in agreement with the qualitative observation in Figs. 4
and 5, RELAX-MORE performs superiorly in terms of reconstruction fidelity,
structure and texture preservation, and noise suppression, outperforming other
state-of-the-art methods.

###  4.3 Ablation Study: Accuracy and Robustness

Figure 7: (a) Phantom and (b) corresponding $T_{1}$ maps generated from Zero-
filled, RELAX-MORE, and fully sampled phantom data. Zero-filled and RELAX-MORE
$T_{1}$ maps are estimated from data undersampled using a $3\times$ 1D
Cartesian undersampling mask. (c) Line profile, indicated in (b), along with
zoom-in views showing the difference in $T_{1}$ values.

An ablation study was carried out further to evaluate the parameter estimation
accuracy of the proposed method. To analyze RELAX-MORE in discrete tissue
environment, it was applied to undersampled data from a phantom composed of
five 20 mL vials consisting of peanut oil and $2\%,4\%,8\%$ Agar dissolved in
DI water and boiled egg white, all immersed in a 200 $\mu$M $\rm{MnCl_{2}}$
water bath (Fig. 7(a)). From left to right, Fig. 7(b) shows the $T_{1}$ maps
obtained from zero-filled data, RELAX-MORE and fully sampled data. Zero-filled
$T_{1}$ map exhibits ripple artifacts due to the undersampling, whereas these
artifacts are eliminated in RELAX-MORE. Comparing the line profile indicated
by the line in Fig. 7(b) and shown in Fig. 7(c), zero-filled $T_{1}$ profile
resembles amplitude modulation, where the periodic amplitude fluctuation stems
from aliasing effects and the higher frequency fluctuation reflects noise. In
good agreement with the fully sampled $T_{1}$ profile, RELAX-MORE’s $T_{1}$
profile not only removes the periodic amplitude fluctuation but also smoothens
the oscillations, demonstrating its ability to average out noise, remove
aliasing artifacts while maintaining high accuracy for parameter estimation in
a wide range.

Figure 8: Estimated $T_{1}$ maps for knee images using $3\times$ 1D Cartesian
undersampling mask.

The ablation study was also carried out to evaluate the parameter estimation
robustness of RELAX-MORE against severe imaging imperfection and the
generalizability of RELAX-MORE on different anatomies such as knee joints.
Fig. 8 shows the $T_{1}$ maps from a representative sagittal slice of the knee
obtained from zero-filled, RELAX-MORE, and fully sampled data. The MRI data
was purposely acquired using a 4-channel receiving-only flex coil, which
provides insufficient coil coverage and a poor signal-to-noise ratio for knee
MRI. As evident in Fig. 8, the fully sampled $T_{1}$ map presents enormous
noises at bone, cartilage, ligament, meniscus, and muscle and the fine
structure, such as the epiphyseal plate (highlighted by a white arrow) is
contaminated by noises. The zero-filled $T_{1}$ map from $3\times$ 1D
Cartesian acceleration also presents severe undersampling artifacts combined
with the noise rendering erroneous $T_{1}$ quantification. However, RELAX-MORE
successfully removes all the artifacts, suppresses the unwanted image noises
in $T_{1}$ quantification, and provides a surprisingly favorable
quantification of different knee joint structures. RELAX-MORE demonstrates its
high robustness and generalizability for reconstructing the in-vivo brain and
in-vivo knee joint, making it a widely applicable method for different
anatomies.

###  4.4 Transfer Learning: Computing Time Efficiency

RELAX-MORE is a subject-specific method that utilizes self-supervised learning
for efficient qMRI reconstruction. While RELAX-MORE can perform well using
single-subject data, as shown in all other experiments, the reconstruction
process (e.g., network training) needs to be conducted for each subject data.
In this experiment, we investigated transfer learning to improve the computing
time efficiency for the training/reconstruction of RELAX-MORE. Using transfer
learning, the network weights after training on one brain data were then
applied as the starting point to train the brain data of another subject.

Figure 9: (a) SSIM calculated as a function of epoch for training with and
without transfer learning. (b) The corresponding estimated $T_{1}$ maps with
(top row) and without (bottom row) transfer learning at epochs 500, 1000, 5000
and 10000, indicated by vertical green lines in (a).

In Fig. 9, the performance of RELAX-MORE with and without transfer learning is
presented by plotting the SSIMs for $T_{1}$ estimations as a function of epoch
number at $3\times$ 1D Cartesian acceleration. Comparing the two SSIM curves
in Fig. 9(a), it is observed that transfer learning starts at a higher initial
SSIM value during the initial training phase and increases faster compared to
the case without transfer learning. This implies that transfer learning
benefits the model to converge more efficiently due to commencing training
with pre-trained weights. This is further confirmed by comparing the $T_{1}$
maps generated with (Fig. 9(b) top row) and without (Fig. 9(b) bottom row)
transfer learning, where at low epoch values (500, 1000), transfer learning
generated $T_{1}$ maps show better $T_{1}$ clarity and contrast between
different tissue types compared to its non-transfer learning counterpart. As
the training epoch increases, both SSIM curves flatten, indicating that the
model has reached stable performance. The results suggest that transfer
learning can improve the computing time efficiency of RELAX-MORE with reduced
training/reconstruction time. At approximately 0.5 seconds per epoch using our
GPU device for training the 3D brain data, transfer learning can reduce the
reconstruction time to less than 10 min to reach good performance. With the
advance of new techniques [44], transfer learning can be a practical approach
to further improve the reconstruction timing efficiency for RELAX-MORE.

##  5 Conclusion

This paper proposes a novel self-supervised learning method, RELAX-MORE, for
qMRI reconstruction. This proposed method uses an optimization algorithm to
unroll a model-based qMRI reconstruction into a deep learning framework,
enabling the generation of highly accurate and robust MR parameter maps at
imaging acceleration. Unlike conventional deep learning methods requiring a
large amount of training data, RELAX-MORE is a subject-specific method that
can be trained on single-subject data through self-supervised learning.
Furthermore, in our experiments, RELAX-MORE outperforms several state-of-the-
art conventional and deep learning methods for accelerated $T_{1}$ maps. This
work also demonstrates several superior aspects of RELAX-MORE in overall
performance, accuracy, robustness, generalizability, and computing timing
efficiency, making this method a promising candidate for advancing accelerated
qMRI for many clinical applications.

## References

  * [1] G. Litjens, T. Kooi, B. E. Bejnordi, A. A. A. Setio, F. Ciompi, M. Ghafoorian, J. A. Van Der Laak, B. Van Ginneken, and C. I. Sánchez, “A survey on deep learning in medical image analysis,” _Medical image analysis_ , vol. 42, pp. 60–88, 2017.
  * [2] D. Shen, G. Wu, and H.-I. Suk, “Deep learning in medical image analysis,” _Annual review of biomedical engineering_ , vol. 19, pp. 221–248, 2017.
  * [3] F. Knoll, K. Hammernik, C. Zhang, S. Moeller, T. Pock, D. K. Sodickson, and M. Akcakaya, “Deep-learning methods for parallel magnetic resonance imaging reconstruction: A survey of the current approaches, trends, and issues,” _IEEE signal processing magazine_ , vol. 37, no. 1, pp. 128–140, 2020.
  * [4] 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.
  * [5] J. Schlemper, J. Caballero, J. V. Hajnal, A. Price, and D. Rueckert, “A deep cascade of convolutional neural networks for mr image reconstruction,” in _Information Processing in Medical Imaging: 25th International Conference, IPMI 2017, Boone, NC, USA, June 25-30, 2017, Proceedings 25_. Springer, 2017, pp. 647–658.
  * [6] Y. Han, J. Yoo, H. H. Kim, H. J. Shin, K. Sung, and J. C. Ye, “Deep learning with domain adaptation for accelerated projection-reconstruction mr,” _Magnetic resonance in medicine_ , vol. 80, no. 3, pp. 1189–1205, 2018.
  * [7] 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, pp. 487–492, 2018.
  * [8] G. Yang, S. Yu, H. Dong, G. Slabaugh, P. L. Dragotti, X. Ye, F. Liu, S. Arridge, J. Keegan, Y. Guo _et al._ , “Dagan: deep de-aliasing generative adversarial networks for fast compressed sensing mri reconstruction,” _IEEE transactions on medical imaging_ , vol. 37, no. 6, pp. 1310–1321, 2017.
  * [9] D. Lee, J. Yoo, S. Tak, and J. C. Ye, “Deep residual learning for accelerated mri using magnitude and phase networks,” _IEEE Transactions on Biomedical Engineering_ , vol. 65, no. 9, pp. 1985–1995, 2018.
  * [10] F. Liu, A. Samsonov, L. Chen, R. Kijowski, and L. Feng, “Santis: sampling-augmented neural network with incoherent structure for mr image reconstruction,” _Magnetic resonance in medicine_ , vol. 82, no. 5, pp. 1890–1904, 2019.
  * [11] B. Yaman, S. A. H. Hosseini, S. Moeller, J. Ellermann, K. Uğurbil, and M. Akçakaya, “Self-supervised learning of physics-guided reconstruction neural networks without fully sampled reference data,” _Magnetic resonance in medicine_ , vol. 84, no. 6, pp. 3172–3191, 2020.
  * [12] M. Blumenthal, G. Luo, M. Schilling, M. Haltmeier, and M. Uecker, “Nlinv-net: Self-supervised end-2-end learning for reconstructing undersampled radial cardiac real-time data,” in _ISMRM annual meeting_ , 2022.
  * [13] H. Z. Wang, S. J. Riederer, and J. N. Lee, “Optimizing the precision in t1 relaxation estimation using limited flip angles,” _Magnetic resonance in medicine_ , vol. 5, no. 5, pp. 399–416, 1987.
  * [14] Y. Zhu, J. Cheng, Z.-X. Cui, Q. Zhu, L. Ying, and D. Liang, “Physics-driven deep learning methods for fast quantitative magnetic resonance imaging: Performance improvements through integration with deep neural networks,” _IEEE Signal Processing Magazine_ , vol. 40, no. 2, pp. 116–128, 2023.
  * [15] C. Cai, C. Wang, Y. Zeng, S. Cai, D. Liang, Y. Wu, Z. Chen, X. Ding, and J. Zhong, “Single-shot t2 mapping using overlapping-echo detachment planar imaging and a deep convolutional neural network,” _Magnetic resonance in medicine_ , vol. 80, no. 5, pp. 2202–2214, 2018.
  * [16] O. Cohen, B. Zhu, and M. S. Rosen, “Mr fingerprinting deep reconstruction network (drone),” _Magnetic resonance in medicine_ , vol. 80, no. 3, pp. 885–894, 2018.
  * [17] H. M. Luu, D.-H. Kim, J.-W. Kim, S.-H. Choi, and S.-H. Park, “qmtnet: Accelerated quantitative magnetization transfer imaging with artificial neural networks,” _Magnetic Resonance in Medicine_ , vol. 85, no. 1, pp. 298–308, 2021.
  * [18] F. Liu, L. Feng, and R. Kijowski, “Mantis: model-augmented neural network with incoherent k-space sampling for efficient mr parameter mapping,” _Magnetic resonance in medicine_ , vol. 82, no. 1, pp. 174–188, 2019.
  * [19] F. Liu, R. Kijowski, L. Feng, and G. El Fakhri, “High-performance rapid mr parameter mapping using model-based deep adversarial learning,” _Magnetic resonance imaging_ , vol. 74, pp. 152–160, 2020.
  * [20] Q. Tian, B. Bilgic, Q. Fan, C. Liao, C. Ngamsombat, Y. Hu, T. Witzel, K. Setsompop, J. R. Polimeni, and S. Y. Huang, “Deepdti: High-fidelity six-direction diffusion tensor imaging using deep learning,” _NeuroImage_ , vol. 219, p. 117017, 2020.
  * [21] R. Feng, J. Zhao, H. Wang, B. Yang, J. Feng, Y. Shi, M. Zhang, C. Liu, Y. Zhang, J. Zhuang _et al._ , “Modl-qsm: Model-based deep learning for quantitative susceptibility mapping,” _NeuroImage_ , vol. 240, p. 118376, 2021.
  * [22] Y. Jun, H. Shin, T. Eo, T. Kim, and D. Hwang, “Deep model-based magnetic resonance parameter mapping network (dopamine) for fast t1 mapping using variable flip angle method,” _Medical Image Analysis_ , vol. 70, p. 102017, 2021.
  * [23] F. Liu, R. Kijowski, G. El Fakhri, and L. Feng, “Magnetic resonance parameter mapping using model-guided self-supervised deep learning,” _Magnetic resonance in medicine_ , vol. 85, no. 6, pp. 3211–3226, 2021.
  * [24] D. Liang, J. Cheng, Z. Ke, and L. Ying, “Deep magnetic resonance image reconstruction: Inverse problems meet neural networks,” _IEEE Signal Processing Magazine_ , vol. 37, no. 1, pp. 141–151, 2020.
  * [25] V. Monga, Y. Li, and Y. C. Eldar, “Algorithm unrolling: Interpretable, efficient deep learning for signal and image processing,” _IEEE Signal Processing Magazine_ , vol. 38, no. 2, pp. 18–44, 2021.
  * [26] M. T. McCann, K. H. Jin, and M. Unser, “Convolutional neural networks for inverse problems in imaging: A review,” _IEEE Signal Processing Magazine_ , vol. 34, no. 6, pp. 85–95, 2017.
  * [27] Y. Yang, J. Sun, H. Li, and Z. Xu, “Deep admm-net for compressive sensing mri,” in _Advances in Neural Information Processing Systems_ , vol. 29. Curran Associates, Inc., 2016\.
  * [28] J. Zhang and B. Ghanem, “Ista-net: Interpretable optimization-inspired deep network for image compressive sensing,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2018, pp. 1828–1837.
  * [29] J. Adler and O. Öktem, “Learned primal-dual reconstruction,” _IEEE transactions on medical imaging_ , vol. 37, no. 6, pp. 1322–1332, 2018.
  * [30] B. Yaman, S. A. H. Hosseini, and M. Akcakaya, “Zero-shot self-supervised learning for MRI reconstruction,” in _International Conference on Learning Representations_ , 2022.
  * [31] M. Akçakaya, S. Moeller, S. Weingärtner, and K. Uğurbil, “Scan-specific robust artificial-neural-networks for k-space interpolation (raki) reconstruction: Database-free deep learning for fast imaging,” _Magnetic resonance in medicine_ , vol. 81, no. 1, pp. 439–453, 2019.
  * [32] H. K. Aggarwal, M. P. Mani, and M. Jacob, “Modl: Model-based deep learning architecture for inverse problems,” _IEEE transactions on medical imaging_ , vol. 38, no. 2, pp. 394–405, 2018.
  * [33] S. Diamond, V. Sitzmann, F. Heide, and G. Wetzstein, “Unrolled optimization with deep priors,” _arXiv preprint arXiv:1705.08041_ , 2017.
  * [34] F. Heide, M. Steinberger, Y.-T. Tsai, M. Rouf, D. Pajak, D. Reddy, O. Gallo, J. Liu, W. Heidrich, K. Egiazarian _et al._ , “Flexisp: A flexible camera image processing framework,” _ACM Transactions on Graphics (ToG)_ , vol. 33, no. 6, pp. 1–13, 2014.
  * [35] S. V. Venkatakrishnan, C. A. Bouman, and B. Wohlberg, “Plug-and-play priors for model based reconstruction,” in _2013 IEEE Global Conference on Signal and Information Processing_. IEEE, 2013, pp. 945–948.
  * [36] W. Bian, Y. Chen, and X. Ye, “Deep parallel mri reconstruction network without coil sensitivities,” in _Machine Learning for Medical Image Reconstruction: Third International Workshop, MLMIR 2020, Held in Conjunction with MICCAI 2020, Lima, Peru, October 8, 2020, Proceedings 3_. Springer, 2020, pp. 17–26.
  * [37] L. I. Sacolick, F. Wiesinger, I. Hancu, and M. W. Vogel, “B1 mapping by bloch-siegert shift,” _Magnetic resonance in medicine_ , vol. 63, no. 5, pp. 1315–1322, 2010.
  * [38] C. Preibisch and R. Deichmann, “Influence of rf spoiling on the stability and accuracy of t1 mapping based on spoiled flash with varying flip angles,” _Magnetic Resonance in Medicine: An Official Journal of the International Society for Magnetic Resonance in Medicine_ , vol. 61, no. 1, pp. 125–135, 2009.
  * [39] M. Uecker, P. Lai, M. J. Murphy, P. Virtue, M. Elad, J. M. Pauly, S. S. Vasanawala, and M. Lustig, “Espirit—an eigenvalue approach to autocalibrating parallel mri: where sense meets grappa,” _Magnetic resonance in medicine_ , vol. 71, no. 3, pp. 990–1001, 2014.
  * [40] X. Glorot and Y. Bengio, “Understanding the difficulty of training deep feedforward neural networks,” in _Proceedings of the thirteenth international conference on artificial intelligence and statistics_. JMLR Workshop and Conference Proceedings, 2010, pp. 249–256.
  * [41] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” in _3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings_ , Y. Bengio and Y. LeCun, Eds., 2015.
  * [42] T. Zhang, J. M. Pauly, and I. R. Levesque, “Accelerating parameter mapping with a locally low rank constraint,” _Magnetic resonance in medicine_ , vol. 73, no. 2, pp. 655–661, 2015.
  * [43] O. Maier, J. Schoormans, M. Schloegl, G. J. Strijkers, A. Lesch, T. Benkert, T. Block, B. F. Coolen, K. Bredies, and R. Stollberger, “Rapid t1 quantification from high resolution 3d data with model-based reconstruction,” _Magnetic resonance in medicine_ , vol. 81, no. 3, pp. 2072–2089, 2019.
  * [44] F. Zhuang, Z. Qi, K. Duan, D. Xi, Y. Zhu, H. Zhu, H. Xiong, and Q. He, “A comprehensive survey on transfer learning,” _Proceedings of the IEEE_ , vol. 109, no. 1, pp. 43–76, 2020.