# Span of regularization for solution of inverse problems, with application to
magnetic resonance relaxometry

Chuan Bi  National Institute on Aging, Baltimore, MD 21224, U.S.A  Yvonne M.
Ou  Department of Mathematical Sciences, University of Delaware, Newark, DE
19716, U.S.A  Mustapha Bouhrara  National Institute on Aging, Baltimore, MD
21224, U.S.A  Richard G. Spencer  National Institute on Aging, Baltimore, MD
21224, U.S.A

> We present a fundamentally new regularization method for the solution of the
> Fredholm integral equation of the first kind, in which we incorporate
> solutions corresponding to a range of Tikhonov regularizers into the end
> result. This method identifies solutions within a much larger function
> space, spanned by this set of regularized solutions, than is available to
> conventional regularization methods. Each of these solutions is regularized
> to a different extent. In effect, we combine the stability of solutions with
> greater degrees of regularization with the resolution of those that are less
> regularized. In contrast, current methods involve selection of a single, or
> in some cases several, regularization parameters that define an optimal
> degree of regularization. Because the identified solution is within the span
> of a set of differently-regularized solutions, we call this method span of
> regularization, or SpanReg. We demonstrate the performance of SpanReg
> through a non-negative least squares analysis employing a Gaussian basis,
> and demonstrate the improved recovery of bimodal Gaussian distribution
> functions as compared to conventional methods. We also demonstrate that this
> method exhibits decreased dependence of the end result on the optimality of
> regularization parameter selection. We further illustrate the method with an
> application to myelin water fraction mapping in the human brain from
> experimental magnetic resonance imaging relaxometry data. We expect SpanReg
> to be widely applicable as an effective new method for regularization of
> inverse problems.

## Significance Statement

Estimating model parameters from data is notoriously difficult, with results
that can exhibit extreme sensitivity to noise. There is a large literature on
limiting this sensitivity through methods known as regularization; these
require an optimal amount of regularization to be defined so that the derived
solution is not over-smoothed. All non-optimal solutions are discarded.
However, there is substantial information content in these discarded
solutions. Accordingly, we present a fundamentally new method, SpanReg, that
incorporates solutions for a range of regularizations into the end result and
shows substantial improvements over existing methods. We apply this to
magnetic resonance imaging myelin mapping in the human brain. We expect
SpanReg to be widely applicable as an effective method for regularization.

##  1 Introduction

###  1.1 Regularization of Inverse Problems Arising from the Fredholm Equation
of the First Kind

A large class of inverse problems arise from a discretized Fredholm integral
equation of the first kind. With the introduction of upper and lower bounds
based on prior knowledge, this may be written as

$y(t)=\int_{\tau_{L}}^{\tau_{U}}A(t,\tau)f(\tau)d\tau$ (1.1)

where y(t) is the measurement or observation, $A(t,\tau)$ is called the kernel
and $f(\tau)$ is the distribution function (DF) to be determined. A
discretization of Eq. (1.1) leads to

$y(t_{i})=\sum_{j=1}^{n}A(t_{i},\tau_{j})f(\tau_{j})\Delta\tau,\quad
i=1,2,\cdots,m$ (1.2)

and

$\mathbf{y}=\mathbf{A}\mathbf{f}$ (1.3)

where $\mathbf{y}\in\mathbb{R}^{m}$, $\mathbf{A}\in\mathbb{R}^{m\times n}$
with $A_{i,j}:=A(t_{i},\tau_{j})\Delta\tau$, and
$\mathbf{f}\in\mathbb{R}^{n}$. Here we assume a uniform discretization
$\triangle\tau:=\frac{\tau_{L}-\tau_{U}}{n}$ and
$\tau_{i}=\tau_{L}+(i-1)\triangle\tau$ for convenience. Depending on the
application, soft or hard constraints on $f(\tau)$ may include, for example,
degree of smoothness, values for upper and lower bounds, and amount of total
variation.

One important application of this formulation is to magnetic resonance
relaxometry (MRR), which estimates the distribution of relaxation times
$T_{2}$ within a sample. We emphasize that this serves as a paradigm for a
much larger class of physical phenomena leading to essentially identical
mathematical considerations. Physically, $T_{2}$, also called the transverse
relaxation time, represents the time constant for decay of transverse
magnetization. The range of possible $T_{2}$ values depends on the sample or
tissue under study, but typical values may range from $\sim$10 ms to 2000 ms
in biomedical studies, with smaller values corresponding to more solid or
rigid tissues. Shorter values require specialized techniques to detect,
depending on available hardware and experimental protocols. The distribution
function $f(T_{2})$ of such values can be of great importance in
characterizing the physical properties and the material composition of a
sample. In MRR, the kernel function is of the special form
$A(t,T_{2})=\exp(-\frac{t}{T_{2}})$ and the Fredholm integral equation in
(1.1) is a Laplace transform. Explicitly, the signal model is

$y(t)=\int e^{-\frac{t}{T_{2}}}f(T_{2})dT_{2}.$ (1.4)

This defines the integrated signal from an ensemble of decaying exponentials,
with the contribution for a given value of the decay constant $T_{2}$
determined by $f(T_{2})$. Thus, as the time $t$ increases, the overall signal
decays into the noise with a mixture of time constants, the distribution of
which is to be determined. Practically, the observed data
$\mathbf{y}_{\text{ob}}\in\mathbb{R}^{m}$ is obtained at measurement times
$\\{t_{i}\\}_{i=1}^{m}$, and the goal is to determine the $T_{2}$ distribution
$\mathbf{f}\in\mathbb{R}^{n}$ as defined by the discretized problem:

$\mathbf{y}_{\text{ob}}=\mathbf{A}\mathbf{f}+\mbox{\boldmath$\omega$},\quad\mathbf{f}\geq\mathbf{0}$
(1.5)

where $\mathbf{A}\in\mathbb{R}^{m\times n}$ is the kernel matrix with entries
$A_{i,j}=e^{-\frac{t_{i}}{T_{2,j}}}\Delta T_{2}$ and $\omega$ is additive
random noise. We define the SNR of the signal by

$\text{SNR}=\frac{\max\left|\mathbf{y}_{\text{ob}}\right|}{\text{RMS}(\mbox{\boldmath$\omega$})}$
(1.6)

where RMS$(\mbox{\boldmath$\omega$})$ is the root mean square noise amplitude.
The positivity constraint $\mathbf{f}\geq\mathbf{0}$ arises from the physical
requirement that $\mathbf{f}$ represents the volume fraction of materials in
the sample. The matrix $\mathbf{A}$ inherits the smoothing property of the
integral operator in (1.1), and can exhibit a condition number that is so
large that any direct inversion of the noisy data $\mathbf{y}_{ob}$ without
regularization, such as through non-negative least squares (NNLS) analysis
[34], can be extremely unstable.

From the perspective of inverse problems, determination of the discretized
$T_{2}$ distribution $\mathbf{f}$ can be seen as a version of the inverse
Laplace transform (ILT), a classic ill-posed problem. An important method for
addressing the numerical instability inherent in this process is through
Tikhonov regularization. In the context of NNLS, where $\mathbf{f}$ is
everywhere non-negative, the native problem

$\mathbf{f}_{LS}=\mathop{\rm
argmin}_{\mathbf{f}\geq\mathbf{0}}\left\lVert\mathbf{A}\mathbf{f}-\mathbf{y}_{\text{ob}}\right\rVert_{2}^{2}$
(1.7)

is replaced by the closely related problem

$\mathbf{f}_{\lambda}=\mathop{\rm
argmin}_{\hat{\mathbf{f}}\geq\mathbf{0}}\left\\{\left\lVert\mathbf{A}\hat{\mathbf{f}}-\mathbf{y}_{\text{ob}}\right\rVert_{2}^{2}+\lambda^{2}\left\lVert\hat{\mathbf{f}}\right\rVert_{2}^{2}\right\\}.$
(1.8)

The second term serves to penalize large values of the norm of the recovered
DF, limiting sensitivity to noise. The regularization parameter $\lambda$ acts
to titrate the relative importance of the two terms in Eq. (1.8). The
recovered DF $\hat{\mathbf{f}}$ may be highly dependent on this parameter;
selection of an optimal $\lambda$ remains the topic of active research. This
formulation in effect replaces the original ill-posed problem, Eq. (1.7), by
Eq. (1.8), a more well-behaved but different problem [21]. Additionally, prior
application-dependent assumptions regarding $\hat{\mathbf{f}}$, such as
sparsity or smoothness, can be introduced into the minimization by
introduction of suitable alternative or additional regularization terms [32,
34].

###  1.2 Methods of Parameter Selection

The problem of selecting $\lambda$ in Tikhonov regularization has been studied
for decades, with no universal approach having been identified. Classical
methods such as the Morozov discrepancy principle (DP) [28], the L-curve [20,
19, 18, 9], and generalized cross-validation (GCV) [15, 35, 36] define
criteria by which to select $\lambda$. For example, the DP seeks a value of
$\lambda$ such that the size of the first term in Eq. (1.8), called the
residual, is matched to the noise level in the observed data. The L-curve
method identifies a value of $\lambda$ that defines the corner of the L-shaped
curve defined by the log-log plot of the solution norm against the norm of the
residual. GCV selects lambda using an expression based on leave-one-out cross-
validation. All of these methods seek to identify one optimal $\lambda$ and
its corresponding regularized solution, and discard all other solutions.
Variations on this include the elastic net, which incorporates two different
$\lambda$’s serving what are essentially distinct roles, but still identifies
single optimal values for each [38]. A survey of methods can be found in [2].

The central theme of the present work is that discarding all results except
the reconstruction corresponding to a single selected $\lambda$ leads to loss
of information that could contribute to the accurate recovery of the DF,
$\mathbf{f}$. We have observed that the effects of regularization, for
example, variations of widths, amplitudes, and shapes of recovered
distributions $\mathbf{f}_{\lambda_{j}}$, depend on the underlying
distribution $\mathbf{f}$. This motivates the notion that the solutions
corresponding to values of $\lambda$ other than the one selected may contain
additional information, so that improved results may be obtained by
incorporating these solutions into the final determination of $\mathbf{f}$
Accordingly, we describe a new method, termed span of regularization, or
SpanReg, for which the recovered distribution function is a linear combination
of regularized solutions across a range of $\lambda$’s. This expands the space
of functions from which the desired DF is drawn.

###  1.3 Gaussian Mixture Representation and Application to Determining
$T_{2}$ Distribution Functions

We use a Gaussian mixture representation as a dictionary to describe the
unknown DF, $\mathbf{f}$, as required for our analysis. One approach would be
to write the DF as the linear combination of a finite set of Gaussian
functions $g(\mu_{i},\sigma_{i})$, with $\mu_{i}$ and $\sigma_{i}$
representing the unknown mean and standard deviation (SD) of a given element
of that set. The discretization of the $g_{i}$’s along the abscissa follows
from the discretization of $\mathbf{f}$ along the abscissa, that is, the
choice of the set of abscissa $T_{2}$ values [6]. With a prior assumption of
the number $M$ of Gaussian components required for an adequate description,
the determination of $\mathbf{f}$ can be recast as the non-linear least
squares problem:

$\mathop{\rm
argmin}_{\left\\{\sigma_{i},\mu_{i}\right\\}}{\left\lVert\mathbf{y}_{\text{ob}}-\mathbf{A}\sum_{i=1}^{M}\mathbf{g}(\sigma_{i},\mu_{i})\right\rVert_{2}^{2}}$

Alternatively, by establishing a dictionary of Gaussian functions of specified
$\mu_{i}$ and $\sigma_{i}$ [12, 17] and incorporating the non-negativity of
$\mathbf{f}_{\text{true}}$, we instead have the problem:

$\mathop{\rm
argmin}_{\mathbf{c}\geq\mathbf{0},\sum_{i}c_{i}=1}{\left\lVert\mathbf{y}_{\text{ob}}-\left(\mathbf{A}\mathbf{G}\right)\mathbf{c}\right\rVert_{2}^{2}}$

where the columns of the matrix $\mathbf{G}$,
$\left\\{\mathbf{g}_{i}\right\\}$, represent an element of the dictionary and
$\mathbf{c}$ is the vector of coefficients respectively assigned to these
elements. The first approach above has the advantage of significantly reducing
the effective dimensionality of the problem [30], but requires the solution of
the highly nonlinear problem of determining the $2M$ variables $\sigma_{i}$
and $\mu_{i}$. On the other hand, the second approach requires only sufficient
knowledge of the system under study to permit a reasonable selection of
dictionary functions; nonlinearity is introduced through the positivity
constraint. This problem becomes increasingly ill-posed with increasing $M$,
and will in general be severely ill-posed due both to the required value of
$M$ and the non-orthogonality of the Gaussian dictionary functions.

##  2 The Span of Regularization Method

###  2.1 Motivation

To motivate SpanReg, we first discuss the artificial case in which
$\mathbf{f}\geq\mathbf{0}$ is known, from which a noise-corrupted signal
$\mathbf{y}_{\text{ob}}=\mathbf{A}\mathbf{f}+\mbox{\boldmath$\omega$}$ can be
created; we now determine $\mathbf{f}_{\lambda_{j}}$ for specific values of
$\lambda_{j}$ within a set of regularization parameters $\Lambda$ by solving
Eq. (1.8).

As indicated above, the various available methods for selection of $\lambda$
in effect evaluate results for several $\lambda$’s and retain only the one
corresponding to an optimal value, discarding the others. Thus, we define our
approximation of $\mathbf{f}$ as a linear combination of regularized
solutions:

$\mathbf{f}\approx\sum_{j=1}^{N}\alpha_{j}\mathbf{f}_{\lambda_{j}}$ (2.1)

where $\left\\{\alpha_{j}\right\\}_{1}^{N}$ is the solution to the least
squares problem

$\alpha=\left(\alpha_{1},\alpha_{2},\cdots,\alpha_{N}\right)=\mathop{\rm
argmin}_{\alpha\geq\mathbf{0}}\left\lVert\mathbf{f}-\sum_{j=1}^{N}\alpha_{j}\mathbf{f}_{\lambda_{j}}\right\rVert_{2}^{2}.$
(2.2)

When $\mathbf{f}$ is known, the set of $\alpha_{j}$’s defined in Eq. (2.2),
with each $\alpha_{j}$ corresponding to an element of a linear combination of
regularized solutions, can always be selected to provide an improvement in
signal reconstruction as compared to selection of any single value of lambda.
This motivates the SpanReg approach for reconstructing an unknown
distribution. However, an algorithm for selecting the best set
$\left\\{\alpha_{j}\right\\}_{{j=}1}^{N}$ in the approximation:

$\mathbf{f}\approx\sum_{j=1}^{N}\alpha_{j}\mathbf{f}_{\lambda_{j}}=:\mathbf{f}_{\alpha},\quad\text{where
}\alpha_{j}\geq 0.$ (2.3)

for an unknown DF, $\mathbf{f}$, must now be developed.

An additional relationship arises from the representation of $\mathbf{f}$ as a
linear combination of Gaussian dictionary functions
$\left\\{\mathbf{g}_{i}\right\\}_{{i=}1}^{M}$, which we impose with a non-
negativity constraint on the expansion coefficients:

$\mathbf{f}\approx\sum_{i=1}^{M}c_{i}\mathbf{g}_{i}=:\mathbf{f}_{\mathbf{c}},\quad\text{where
}c_{i}\geq 0.$ (2.4)

We consider the case in which the integral over $\mathbf{f}$ equals $1$, so
that it represents a probability distribution function (PDF). Since $\int
g\left(\sigma_{i},\mu_{i}\right)=1$, we have the requirement that $\sum
c_{i}=1$.

###  2.2 Analytic Framework

####  2.2.1 Theory

To link the two expressions (2.3) and (2.4), we define a symbolic inversion
operator $\mathbf{A}^{-1}_{\lambda_{j}}$ as follows: given a vector of noisy
observations $\mathbf{y}_{\text{ob}}$,
$\mathbf{A}^{-1}_{\lambda_{j}}\mathbf{y}_{\text{ob}}:=\mathbf{f}_{\lambda_{j}}$,
where $\mathbf{f}_{\lambda_{j}}$ is the solution of (1.8) with
$\lambda=\lambda_{j}$, $j=1,\cdots,N$. This operator is not the usual matrix
pseudoinverse because of the non-negativity constraints.

We apply this to the noisy observations that would correspond to the signal
generated by the $\mathbf{g}_{i}$, obtaining a set of associated DF’s
$\\{\mathbf{g}_{i,\lambda_{j}}\\}$, $i=1,\cdots,M$, $j=1,\cdots,N$ according
to Eq. (1.8):

$\mathbf{g}_{i,\lambda_{j}}=\mathbf{A}^{-1}_{\lambda_{j}}\left(\mathbf{A}\mathbf{g}_{i}+\mbox{\boldmath$\omega$}\right).$
(2.5)

where $\omega$ represents a noise realization exhibiting the same noise level,
defined by RMS($\omega$), as the observed data $\mathbf{y}_{\text{ob}}$. Note
that while determining noise amplitude within an experimental signal can be
highly problematic, in our examples of decaying exponentials this information
is directly available from data collected from the signal tail.

We now seek to approximate both $\mathbf{f}_{\alpha}$ and
$\mathbf{f}_{\mathbf{c}}$ by the $\\{\mathbf{g}_{i,\lambda_{j}}\\}$. We first
write

$\mathbf{g}_{i}\approx\sum_{j=1}^{N}\beta_{ij}\mathbf{g}_{i,\lambda_{j}},\quad\beta_{ij}\geq
0.$ (2.6)

The expansion coefficients $\beta_{ij}$ would be Kronecker $\delta$’s in the
absence of noise and regularization. As $\omega$ is random, the corresponding
$\\{\mathbf{g}_{i,\lambda_{j}}\\}$ and $\\{\beta_{ij}\\}$ are random
variables. Let $\langle\cdot\rangle$ denote the ensemble average of a random
variable over $n_{\text{run}}$ noise realizations. In particular,

$\langle\mathbf{g}_{i,\lambda_{j}}\rangle=\frac{1}{n_{\text{run}}}\sum_{k=1}^{n_{\text{run}}}\mathbf{g}_{i,\lambda_{j}}^{(k)},\quad\langle\beta_{ij}\rangle=\frac{1}{n_{\text{run}}}\sum_{k=1}^{n_{\text{run}}}\beta_{ij}^{(k)}$

where the superscript $(k)$ indicates the $k$-th noise realization. We now
apply a representation analogous to Eq. (2.4) to define regularized
approximations:

$\mathbf{f}_{\lambda_{j}}\approx\sum_{i=1}^{M}x_{ij}\langle\mathbf{g}_{i,\lambda_{j}}\rangle$
(2.7)

where the $\\{x_{ij}\\}$ can be obtained through least squares analysis.

Thus, given the dictionary of Gaussian distributions
$\left\\{\mathbf{g}_{i}\right\\}_{i=1}^{M}$, we can find the coefficients
$\left\\{\langle\beta_{ij}\rangle\right\\}$ and $\left\\{x_{ij}\right\\}$ from
the approximations to the regularized solutions
$\left\\{\mathbf{f}_{\lambda_{j}}\right\\}$ and
$\left\\{\langle\mathbf{g}_{i,\lambda_{j}}\rangle\right\\}$, respectively. Now
equating (2.3) and (2.4) and using (2.6) and (2.7), we arrive at an expression
containing only $\\{\langle\beta_{ij}\rangle\\}$, $\\{x_{ij}\\}$ and
$\\{\langle\mathbf{g}_{i,\lambda_{j}}\rangle\\}$:

$\displaystyle\sum_{j=1}^{N}\alpha_{j}\mathbf{f}_{\lambda_{j}}$
$\displaystyle\approx\sum_{i=1}^{M}c_{i}\mathbf{g}_{i}$ (2.8)
$\displaystyle\Rightarrow\mathbf{f}_{\alpha}:=\sum_{j=1}^{N}\alpha_{j}\sum_{i=1}^{M}x_{ij}\langle\mathbf{g}_{i,\lambda_{j}}\rangle$
$\displaystyle\approx\sum_{i=1}^{M}c_{i}\sum_{j=1}^{N}\langle\beta_{ij}\rangle\langle\mathbf{g}_{i,\lambda_{j}}\rangle=:\mathbf{f}_{\mathbf{c}}$

The coefficients $\bm{\alpha}$ and $\mathbf{c}$ are obtained by solving the
least squares problem:

$\left\\{\begin{aligned} (\bm{\alpha}^{*},\mathbf{c}^{*})=\mathop{\rm
argmin}\left\lVert\mathbf{f}_{\alpha}-\mathbf{f}_{\mathbf{c}}\right\rVert_{2}\\\
\text{subject to
}\bm{\alpha^{*}}\geq\mathbf{0},\;\mathbf{c}^{*}\geq\mathbf{0},\text{ and
}\sum_{i}c^{*}_{i}=1,\end{aligned}\right.$ (2.9)

Our expression for the final recovered $\mathbf{f}^{*}$ is then:

$\mathbf{f}^{*}_{\alpha}=\sum_{j=1}^{N}\alpha_{j}{{}^{*}}\mathbf{f}_{\lambda_{j}}$
(2.10)

Eq. 2.10 is the main result of our analysis, defining the desired recovered DF
in terms of a linear combination of differently-regularized solutions based on
the observed data. Alternatively, the corresponding result for
$\mathbf{f}^{*}_{c}$ may also be used with essentially equivalent results.

####  2.2.2 Numerical Implementation

We define the following notation:

$\displaystyle\mathbf{L}_{\alpha}$
$\displaystyle=\begin{pmatrix}\begin{pmatrix}\langle\mathbf{g}_{1,\lambda_{1}}\rangle&\langle\mathbf{g}_{2,\lambda_{1}}\rangle&\cdots&\langle\mathbf{g}_{M,\lambda_{1}}\rangle\end{pmatrix}&,\cdots,&\begin{pmatrix}\langle\mathbf{g}_{1,\lambda_{N}}\rangle&\langle\mathbf{g}_{2,\lambda_{N}}\rangle&\cdots&\langle\mathbf{g}_{M,\lambda_{N}}\rangle\end{pmatrix}\end{pmatrix}\in\mathbb{R}^{n\times
MN}$ $\displaystyle\mathbf{L}_{c}$
$\displaystyle=\begin{pmatrix}\begin{pmatrix}\langle\mathbf{g}_{1,\lambda_{1}}\rangle&\langle\mathbf{g}_{1,\lambda_{2}}\rangle&\cdots&\langle\mathbf{g}_{1,\lambda_{N}}\rangle\end{pmatrix}&,\cdots,&\begin{pmatrix}\langle\mathbf{g}_{M,\lambda_{1}}\rangle&\langle\mathbf{g}_{M,\lambda_{2}}\rangle&\cdots&\langle\mathbf{g}_{M,\lambda_{N}}\rangle\end{pmatrix}\end{pmatrix}\in\mathbb{R}^{n\times
MN}$ $\displaystyle\mathbf{x}_{\text{vec}}$
$\displaystyle=\begin{pmatrix}\begin{pmatrix}x_{11}&x_{21}&\cdots&x_{M1}\end{pmatrix}&,\cdots,&\begin{pmatrix}x_{1N}&x_{2N}&\cdots&x_{MN}\end{pmatrix}\end{pmatrix}\in\mathbb{R}^{MN}$
$\displaystyle\bm{\beta}_{\text{vec}}$
$\displaystyle=\begin{pmatrix}\begin{pmatrix}\langle\beta_{1,\lambda_{1}}\rangle&\langle\beta_{1,\lambda_{2}}\rangle&\cdots&\langle\beta_{1,\lambda_{N}}\rangle\end{pmatrix}&,\cdots,&\begin{pmatrix}\langle\beta_{M,\lambda_{1}}\rangle&\langle\beta_{M,\lambda_{2}}\rangle&\cdots&\langle\beta_{M,\lambda_{N}}\rangle\end{pmatrix}\end{pmatrix}\in\mathbb{R}^{MN}$

where for each pair of fixed $i$ and $j$, the regularized approximation
$\langle\mathbf{g}_{i,\lambda_{j}}\rangle$ is a column vector of dimension
$\mathbb{R}^{n}$ and the elements $x_{ij}$ and
$\langle\beta_{i\lambda_{j}}\rangle$ are scalars. Thus, the lengths of the
vectors $\mathbf{x}_{\text{vec}}$ and $\bm{\beta}_{\text{vec}}$ are equal to
the number of columns in $\mathbf{L}_{\alpha}$ and $\mathbf{L}_{c}$,
respectively. Moreover, we write for the unknowns:

$\bm{\alpha}_{\text{vec}}=\begin{pmatrix}\begin{pmatrix}\alpha_{1}&\alpha_{1}&\cdots&\alpha_{1}&\alpha_{1}\end{pmatrix}_{M}&,\cdots,&\begin{pmatrix}\alpha_{N}&\alpha_{N}&\cdots&\alpha_{N}&\alpha_{N}\end{pmatrix}_{M}\end{pmatrix}^{T}\in\mathbb{R}^{MN}$
$\mathbf{c}_{\text{vec}}=\begin{pmatrix}\begin{pmatrix}c_{1}&c_{1}&\cdots&c_{1}&c_{1}\end{pmatrix}_{N}&,\cdots,&\begin{pmatrix}c_{M}&c_{M}&\cdots&c_{M}&c_{M}\end{pmatrix}_{N}\end{pmatrix}^{T}\in\mathbb{R}^{MN}$

Then Eqs. (2.3) and (2.4) can be written:

$\mathbf{f}_{\alpha}=\mathbf{L}_{\alpha}\cdot\text{diag}\left(\mathbf{x}_{\text{vec}}\right)\cdot\bm{\alpha}_{\text{vec}}$
$\mathbf{f}_{c}=\mathbf{L}_{c}\cdot\text{diag}\left(\bm{\beta}_{\text{vec}}\right)\cdot\mathbf{c}_{\text{vec}}$

where $\cdot$ denotes the usual matrix-vector multiplication and diag
indicates the diagonal matrix formed from the vector argument. Eq. (2.9) can
then be re-formulated as:

$\left\\{\begin{aligned} (\mathbf{c}^{*},\bm{\alpha}^{*})&=\mathop{\rm
argmin}\left\lVert\mathbf{L}_{\alpha}\cdot\text{diag}\left(\mathbf{x}_{\text{vec}}\right)\cdot\bm{\alpha}_{\text{vec}}-\mathbf{L}_{c}\cdot\text{diag}\left(\bm{\beta}_{\text{vec}}\right)\cdot\mathbf{c}_{\text{vec}}\right\rVert_{2}\\\
&\text{subject to
}\mathbf{c}\geq\mathbf{0},\bm{\alpha}\geq\mathbf{0},\sum_{i}c_{i}=1\end{aligned}\right.$
(2.11)

By writing the solution in the stacked form:

$\mathbf{s}=\left(\alpha_{1},\alpha_{2},\cdots,\alpha_{N},c_{1},c_{2},\cdots,c_{M}\right)^{T}\in\mathbb{R}^{N+M,}$
(2.12)

(2.11) can be expressed as a conventional LS problem for
$\mathbf{s}^{*}\in\mathbb{R}^{N+M}$

$\left\\{\begin{aligned} \mathbf{s}^{*}&=\mathop{\rm
argmin}_{\mathbf{s}\geq\mathbf{0}}\left\lVert\mathbf{B}\mathbf{s}\right\rVert_{2}\\\
&\text{such that }\sum_{j=N+1}^{N+M}s_{j}=1\end{aligned}\right.$ (2.13)

with

$\mathbf{B}=\mathbf{L}_{\alpha}\cdot\text{diag}\left(\mathbf{x}_{\text{vec}}\right)\cdot\mathbf{TT}_{\alpha}-\mathbf{L}_{c}\cdot\text{diag}\left(\bm{\beta}_{\text{vec}}\right)\cdot\mathbf{TT}_{c}\in\mathbb{R}^{n\times(N+M)}$
(2.14) $\mathbf{TT}_{\alpha}=\mathbf{I}_{N}\otimes\begin{pmatrix}1\\\ 1\\\
\vdots\\\ 1\end{pmatrix}_{M}\begin{pmatrix}\mathbf{I}_{N}&\mathbf{0}_{N\times
M}\end{pmatrix}\in\mathbb{R}^{MN\times(N+M)}$ (2.15)

and

$\mathbf{TT}_{c}=\mathbf{I}_{M}\otimes\begin{pmatrix}1\\\ 1\\\ \vdots\\\
1\end{pmatrix}_{N}\begin{pmatrix}\mathbf{0}_{M\times
N}&\mathbf{I}_{M}\end{pmatrix}\in\mathbb{R}^{MN\times(N+M)}$ (2.16)

where $\otimes$ is the Kronecker tensor product, $\mathbf{I}_{M}$,
$\mathbf{I}_{N}$ are the identity matrices of rank $M$ and $N$, respectively,
and $\mathbf{0}_{N\times M}$ is the zero matrix of size $N\times M$.

The computational procedure can be divided into two parts, which we refer to
as the offline part, meaning independent of the actual data set, and the data-
dependent online part. In the offline computation,
$\left\\{\mathbf{g}_{i}\right\\}$,
$\left\\{\langle\mathbf{g}_{i,\lambda_{j}}\rangle\right\\}$ and
$\left\\{\langle\beta_{i,\lambda_{j}}\rangle\right\\}$ are determined only
once. For the online part, for each noisy measurement
$\mathbf{y}_{\text{ob}}$, the corresponding
$\left\\{\mathbf{f}_{\lambda_{j}}\right\\}$, $\left\\{x_{ij}\right\\}$ and
$(\mathbf{c}^{*},\bm{\alpha}^{*})$ can be obtained with the desired solution
given by (2.10). Note that the $\mathbf{g}_{i,\lambda_{j}}$, and hence the
$\beta_{i,\lambda_{j}}$ in (2.6), are noise-dependent. For computations, we
use the ensemble average over $n_{\text{run}}$ of realizations of Eq. (2.6).

The pseudocode for the offline and online computations reads as follows:

Data: Pre-determined $\left\\{\lambda_{j}\right\\}$,
$\left\\{\mathbf{g}_{i}\right\\}$, $n_{\text{run}}$ and noise SNR

Result: $\left\\{\langle\mathbf{g}_{i,\lambda_{j}}\rangle\right\\}$ and
$\left\\{\langle\beta_{i,\lambda_{j}}\rangle\right\\}$

initialization;

for  _$k=1,2,\cdots,n_{\text{run}}$ ,_ do

generate noise realization $\omega_{k}$ for  _$i=1,2,\cdots,M$_ do

$\mathbf{y}_{\text{clean}}^{(i)}=A\mathbf{g}_{i}$;

$\mathbf{y}_{\text{noisy}}^{(i)}(\omega_{k})=\mathbf{y}_{\text{clean}}^{(i)}+\omega_{k}$;

for  _$j=1,2,\cdots,N$_ do

$\mathbf{g}_{i,\lambda_{j}}(\omega_{k})=\mathop{\rm
argmin}_{\mathbf{g}\geq\mathbf{0}}\left\\{\left\lVert\mathbf{A}\mathbf{g}_{i}-\mathbf{y}_{\text{noisy}}^{(i)}(\omega_{k})\right\rVert_{2}^{2}+\lambda_{j}^{2}\left\lVert\mathbf{g}_{i}\right\rVert_{2}^{2}\right\\}$;

end for

end for

$\beta_{ij}(\omega_{k})=\mathop{\rm argmin}_{b_{ij}\geq
0}\left\lVert\mathbf{g}_{i}-\sum_{j=1}^{N}b_{ij}\mathbf{g}_{i,\lambda_{j}}(\omega_{k})\right\rVert_{2}$

end for

$\langle\mathbf{g}_{i,\lambda_{j}}\rangle=\frac{1}{n_{\text{run}}}\sum_{k=1}^{n_{\text{run}}}\mathbf{g}_{i,\lambda_{j}}(\omega_{k})$
and
$\langle\beta_{i,\lambda_{j}}\rangle=\frac{1}{n_{\text{run}}}\sum_{k=1}^{n_{\text{run}}}\beta_{ij}(\omega_{k})$

Algorithm 1 Offline Computation

Data: $\mathbf{y}_{\text{ob}}$, $\left\\{\lambda_{j}\right\\}$,
$\left\\{\langle\mathbf{g}_{i,\lambda_{j}}\rangle\right\\}$,
$\left\\{\langle\beta_{i,\lambda_{j}}\rangle\right\\}$

Result: $\mathbf{f}^{*}$

initialization;

for  _$j=1,2,\cdots,N$_ do

$\mathbf{f}_{\lambda_{j}}=\mathop{\rm
argmin}_{\hat{\mathbf{f}}\geq\mathbf{0}}\left\\{\left\lVert\mathbf{A}\hat{\mathbf{f}}-\mathbf{y}_{\text{ob}}\right\rVert_{2}^{2}+\lambda_{j}^{2}\left\lVert\hat{\mathbf{f}}\right\rVert_{2}^{2}\right\\}$;

$\mathbf{x}_{\cdot,j}=\mathop{\rm argmin}_{x_{ij}\geq
0}\left\lVert\mathbf{f}_{\lambda_{j}}-\sum_{i=1}^{M}x_{ij}\langle\mathbf{g}_{i,\lambda_{j}}\rangle\right\rVert_{2}$;

end for

$(\mathbf{c}^{*},\bm{\alpha}^{*})=\mathop{\rm
argmin}_{\mathbf{c},\bm{\alpha}\geq
0,\sum_{i}c_{i}=1}\left\lVert\sum_{j=1}^{N}\alpha_{j}\sum_{i=1}^{M}x_{ij}\langle\mathbf{g}_{i,\lambda_{j}}\rangle-\sum_{i=1}^{M}c_{i}\sum_{j=1}^{N}\langle\beta_{i,\lambda_{j}}\rangle\langle\mathbf{g}_{i,\lambda_{j}}\rangle\right\rVert_{2}$;

and $\mathbf{f}^{*}_{\alpha}=\sum_{j}\alpha^{*}_{j}\mathbf{f}_{\lambda_{j}}$;

Algorithm 2 Online Computation

##  3 Applications of SpanReg to One-dimensional Magnetic Resonance
Relaxometry

We illustrate the application of SpanReg to the inverse problem of MRR for
transverse relaxation in one dimension, whose discretized version is described
in Section 1.1. The objective is to approximate the distribution of transverse
relaxation times $T_{2}$ within a sample from the signal $\mathbf{y}_{ob}$. In
the context of our formalism, this distribution function takes the role of
$\mathbf{f}$. $T_{2}$ values will be expressed in milliseconds (ms).

###  3.1 Parameter Settings for the Implementation of SpanReg

Several user-defined parameters must be selected before proceeding with the
off-line computation. We first fix the values of $N$, $M$ and
$n_{\text{run}}$, which are the number of regularization parameters, the
number of Gaussian functions in the dictionary and the number of noise
realizations, respectively. The values of $\lambda_{j}$, $j=1,\cdots,N$ and
the means and standard deviations (SD) of the Gaussian functions
$\left\\{\mathbf{g}_{i}\right\\}_{i=1}^{M}$ in the dictionary must also be
selected.

The regularization parameters are determined based on L-curve analysis, with
particular attention to parameters that provide a wide range of distinct
solutions to Eq. (1.8). In the following simulations, we chose $N$ = 16 with
$\lambda_{1},\cdots,\lambda_{N}$ logarithmically spaced over the interval
$[10^{-6},10^{1}]$. The $T_{2}$ axis range is from 1 - 200 ms, discretized at
1 ms intervals, i.e. $\triangle T_{2}=1$. We chose a dictionary consisting of
three families of Gaussians, each of which has its mean values equally spaced
along the $T_{2}$ axis, and has a specified SD. The three families
respectively consist of 160 members with SD$=2$ ms, 40 members with SD$=3$ ms,
and 20 members with SD$=4$ ms; $M=160+40+20=220$.

###  3.2 Comparison of SpanReg and Classical Parameter Selection Methods

We evaluated the performance of SpanReg for the reconstruction of a range of
DFs consisting of two Gaussian functions with means
$\left(\mu_{1},\mu_{2}\right)$ and standard deviations
$\left(\sigma_{1},\sigma_{2}\right):$

$f_{\text{sim}}(T_{2};\\{\mu_{1},\mu_{2},\sigma_{1},\sigma_{2}\\})=\frac{1}{\sqrt{2\pi\sigma_{1}^{2}}}e^{-\frac{(T_{2}-\mu_{1})^{2}}{2\sigma_{1}^{2}}}+\frac{1}{\sqrt{2\pi\sigma_{2}^{2}}}e^{-\frac{(T_{2}-\mu_{2})^{2}}{2\sigma_{2}^{2}}}.$
(3.1)

We examine DF’s with $\sigma_{1}$ taking values in the set
$\\{2+0.75k,\,k=0,1,2,3,4\\}$ (ms) and $\sigma_{2}=3\sigma_{1}$. The
separation between the Gaussians is defined by their ratio of peak separation
(RPS), defined as RPS=$\frac{\mu_{2}}{\mu_{1}}$. We evaluate DF’s with
$RPS=1+0.75k$, $k=0,1,2,3,4$ for each $(\sigma_{1},\sigma_{2})$ pair. Note
that none of the Gaussian functions in these distributions is included in the
Gaussian dictionary $\left\\{\mathbf{g}_{i}\right\\}$. We compare the recovery
of each of the 25 resulting DF’s (Fig. 1) obtained from SpanReg and from the
DP; here and elsewhere, by DP we indicate NNLS using the DP to set the
Tikhonov regularization parameter. This is a particularly suitable comparison
since for both methods, the noise level in the data must be known or
estimated. As noted, the noise level can be estimated with good accuracy in
the case of the superposition of exponentially decaying signals such as in an
MRR experiment. The discretization along the $T_{2}$ axis is defined as above,
with $n=200$ evenly spaced values within the range $[1,200]$ ms. The
measurement time vector is defined as $m=150$ sampling times evenly
distributed in $t\in[0.3,400]$ ms. The additive noise level is set to result
in SNR of $500$. The range of $\lambda$’s used for Tikhonov regularization was
from $10^{-6}$ to $10$ with $N=16$ logarithmically spaced values.

For the implementation of DP, the selected $\lambda_{DP}$ is determined by the
following criterion [21]:

$\text{Choose }\lambda=\lambda_{DP}\text{ such that
}\left\lVert\mathbf{A}\mathbf{f}_{\lambda_{DP}}-\mathbf{y}_{\text{ob}}\right\rVert_{2}=\nu_{DP}\left\lVert\mbox{\boldmath$\omega$}\right\rVert_{2},$
(3.2)

where $\nu_{DP}\geq 1$ and
$\left\lVert\mbox{\boldmath$\omega$}\right\rVert_{2}\approx\sqrt{m}\sigma(\mbox{\boldmath$\omega$})$.
We chose the safety factor $\nu_{DP}$ to be $\nu_{DP}=1.05$.

Fig. 1 shows the comparison of SpanReg with the DP for recovery of the
indicated underlying DF for a single noise realization. As seen, SpanReg
exhibits greater ability to resolve the two components of the DF and
accurately model their amplitudes.

Figure 1: Comparison of the recovered distributions from SpanReg (red) and the
discrepancy principle (blue) with the true underlying distributions (black).
Each panel represents a $T_{2}$ distribution consisting of the sum of two
Gaussian functions as in Eq. (3.1). The value of $\sigma_{1}$, the standard
deviation of the left-most component, increases from 2 ms (top row) to 5 ms
(bottom row) with $\sigma_{2}=3\sigma_{1}$ throughout. Across each row, the
ratio of peak separation (RPS) increases from 1 (left-most column) to 4
(right-most column), with the left-most distribution centered at $\mu=35$ ms.

Figure 2 presents the comparison of the two reconstruction methods in terms of
heat maps of mean relative error, defined respectively for SpanReg and the DP
by

$\varepsilon_{\text{SpanReg}}=\frac{\left\lVert\mathbf{f}^{*}_{\text{SpanReg}}-\mathbf{f}_{\text{sim}}\right\rVert_{2}}{\left\lVert\mathbf{f}_{\text{sim}}\right\rVert_{2}},\,\varepsilon_{\text{DP}}=\frac{\left\lVert\mathbf{f}^{*}_{\text{DP}}-\mathbf{f}_{\text{sim}}\right\rVert_{2}}{\left\lVert\mathbf{f}_{\text{sim}}\right\rVert_{2}}\mbox{
and
}\varepsilon_{\text{diff}}=\varepsilon_{\text{SpanReg}}-\varepsilon_{\text{DP}}$
(3.3)

where $\mathbf{f}^{*}$ is the reconstructed distribution and
$\mathbf{f}_{\text{sim}}$ is defined in Eq. (3.1). Thus, in terms of relative
errors, a negative $\varepsilon_{\text{diff}}$ indicates the superiority of
SpanReg reconstruction.

(a)  (b)

(c)

Figure 2: Heat maps showing relative errors for SpanReg (upper-left),
discrepancy principle (upper-right), and their differences
$\varepsilon_{\text{diff}}$ (bottom), as defined in the text.

We see that by this summary metric of relative error, SpanReg outperforms DP
across all DF’s studied. Indeed, the errors for SpanReg are fairly constant
across distribution, while the errors for DP increase as RPS increases over
the illustrated range. Thus, as the component centers deviate from each other,
including when closely-spaced, SpanReg is more able to accurately resolve
them. In addition, both methods tend to perform better as the widths of the
Gaussian distributions $\sigma$ of the target DF increase; this highlights the
smoothing characteristics of the $L_{2}$ norm.

Additional insights are provided by comparing the heat map for the relative
error metric described in Eq. (3.3), of use in summarizing the quality of
reconstructions across a wide range of target DFs, with Fig. 1. From the
latter we see that for closely-spaced Gaussian components, represented by the
second column of Fig. 1, the reconstructions are not only quantitatively, but
also qualitatively, different. The DP reconstruction incorrectly provides a
single-component reconstruction in a number of cases, while SpanReg is clearly
able to resolve the two underlying components. Therefore, it is important to
examine the actual recovered $T_{2}$ distributions in addition to comparing
the single summary metric defined by relative error.

This is further highlighted in Figs. 3 and 4, comparing SpanReg and DP on two
target distributions with respectively a greater and a lesser RPS. We show the
results of reconstructions over $10$ noise realizations, as well as a
comparisons of the corresponding signals generated by these reconstructions.

Figure 3: Reconstruction of the sum of two Gaussians with
$(\mu_{1},\sigma_{1}$)=(30 ms, 3 ms) and $(\mu_{2},\sigma_{2})$=(120 ms, 5ms),
shown in black in the top and middle rows; SNR=500. Each column corresponds to
a single noise realization. Top: DP (blue), Middle: SpanReg (red), Bottom:
signals generated from the DP solution and the SpanReg solution. Figure 4:
Reconstruction of the sum of two Gaussians with $(\mu_{1},\sigma_{1}$)=(30 ms,
3 ms) and $(\mu_{2},\sigma_{2})$=(50 ms, 5 ms), shown in black in the top and
middle rows; SNR=500. Each column corresponds to a single noise realization.
Top: DP (blue), Middle: SpanReg (red), Bottom: signals generated from the DP
solution and the SpanReg solution.

In Fig. 3, SpanReg, as compared to the DP, provides a much more accurate
reconstruction of the two components in the DF, although both methods
accurately recover two resolved components. In the much more ill-posed problem
illustrated in Fig. 4, SpanReg clearly outperforms the DP in terms of
stability with respect to noise and ability to resolve two closely-spaced
peaks. DP fails to do this for all of the 10 noise realizations, while SpanReg
succeeds in 7 out of 10 cases. Decreasing the DP safety factor improves the
resolution performance of the DP, but introduces spurious peaks and greater
noise sensitivity.

###  3.3 Stability With Respect to Choice of Regularization Parameters

Selection of the optimal regularization parameter $\lambda$ is of critical
importance in Tikhonov regularization, with deviations from that optimal value
potentially resulting in substantial differences in the reconstructed DF. Fig.
6 from the Supplementary Information shows the solutions obtained from
Tikhonov regularization, and from SpanReg over different ranges of
$\lambda$’s, including in all cases $\lambda_{\text{opt}}$. In this simulation
with a known underlying distribution, $\lambda_{\text{opt}}$ is determined by
the $L_{2}$ error metric describing the difference between the underlying DF
and the reconstructed DF:

$\lambda_{\text{opt}}=\mathop{\rm
argmin}_{\lambda}\left\\{\left\lVert\mathbf{f}_{\lambda}-\mathbf{f}_{\text{sim}}\right\rVert_{2}\right\\}.$
(3.4)

Once $\lambda_{opt}$ is determined, stability is quantified by comparing the
DF recovered using the different values of $\lambda$ from the set
$\mathcal{S}=\left\\{\frac{\lambda_{\text{opt}}}{2^{4}},\cdots,\frac{\lambda_{\text{opt}}}{2^{1}},\frac{\lambda_{\text{opt}}}{2^{0}},\frac{\lambda_{\text{opt}}}{2^{-1}},\cdots,\frac{\lambda_{\text{opt}}}{2^{-5}}\right\\}$.
To apply this analysis to the DP, we simply compare solutions obtained with
the different values of $\lambda$. We also show the differences between the
results for these $\lambda$’s and the result for
$\mathbf{f}_{\lambda_{\text{opt}}}$; larger values in these difference plots
indicate less stability. For SpanReg, we perform the analysis using different
triads of adjacent values of $\lambda$ from the $\mathcal{S}$. We form the
reconstruction using the three $\lambda$’s centered around
$\lambda_{\text{opt}}$, denoted
$\mathbf{f}_{\text{SpanReg},\lambda_{\text{opt}}}$:

$\mathbf{f}_{\text{SpanReg},\lambda_{\text{opt}}}:=\alpha_{1}\mathbf{f}_{\frac{\lambda_{\text{opt}}}{2}}+\alpha_{2}\mathbf{f}_{\lambda_{\text{opt}}}+\alpha_{3}\mathbf{f}_{\frac{\lambda_{\text{opt}}}{2^{-1}}},$
(3.5)

where the values $\left\\{\alpha_{1},\alpha_{2},\alpha_{3}\right\\}$ are
obtained from SpanReg. Then, to evaluate stability with respect to choice of
$\lambda$, we plot results for SpanReg using three shifted values of
$\lambda$’s. For example, for an integer $-3\leq j\leq 4$,

$\mathbf{f}_{\text{SpanReg},\frac{\lambda_{\text{opt}}}{2^{j}}}:=\alpha_{1}\mathbf{f}_{\text{SpanReg},\frac{\lambda_{\text{opt}}}{2^{j+1}}}+\alpha_{2}\mathbf{f}_{\text{SpanReg},\frac{\lambda_{\text{opt}}}{2^{j}}}+\alpha_{3}\mathbf{f}_{\text{SpanReg},\frac{\lambda_{\text{opt}}}{2^{j-1}}}.$

The comparison of the results for shifted triples of $\lambda$s defines the
stability of SpanReg with respect to choice of regularization parameter. We
also plot the difference between the DF recovered by adjacent shifted triplets
to the the recovery obtained using the triplet centered on
$\mathbf{f}_{\lambda_{\text{opt}}}$, i.e. we define
$\Delta\mathbf{f}_{\frac{\lambda_{\text{opt}}}{2^{j}}}=\mathbf{f}_{\lambda_{\text{opt}}}-\mathbf{f}_{\frac{\lambda_{\text{opt}}}{2^{j}}}$,
and similarly,
$\Delta\mathbf{f}_{\text{SpanReg},\frac{\lambda_{\text{opt}}}{2^{j}}}=\mathbf{f}_{\text{SpanReg},\lambda_{\text{opt}}}-\mathbf{f}_{\text{SpanReg},\frac{\lambda_{\text{opt}}}{2^{j}}}$.
These plots quantify the changes in the recovered DF resulting from a
suboptimal selection of $\lambda$ in SpanReg. A larger magnitude in these
difference plots corresponds to lower stability with respect to choice of
$\lambda$ and vice versa.

An example of this is shown in the Supplementary Information, Fig. 6, where
results for the illustrated two-component Gaussian DF are shown for for
optimally regularized DP solutions, and solutions obtained with sub-optimal
regularization (left panels). Corresponding solutions obtained with SpanReg
using different sub-sequences of regularized solutions are also shown (right
panels). Departures from the optimal regularization solutions obtained for
Tikhonov regularization (lower left panel) are substantially larger than those
using SpanReg (lower right panel). The corresponding $L_{2}$ norms of the
differences are
$\left(\left\lVert\Delta\mathbf{f}_{\lambda_{1}}\right\rVert,\left\lVert\Delta\mathbf{f}_{\lambda_{2}}\right\rVert,\left\lVert\Delta\mathbf{f}_{\lambda_{4}}\right\rVert,\left\lVert\Delta\mathbf{f}_{\lambda_{5}}\right\rVert\right)=\left(0.23,0.11,0.08,0.13\right)$
for Tikhonov regularization and
$\left(\left\lVert\Delta\mathbf{f}_{\text{SpanReg},\lambda_{1}}\right\rVert,\left\lVert\Delta\mathbf{f}_{\text{SpanReg},\lambda_{2}}\right\rVert,\left\lVert\Delta\mathbf{f}_{\text{SpanReg},\lambda_{4}}\right\rVert,\left\lVert\Delta\mathbf{f}_{\text{SpanReg},\lambda_{5}}\right\rVert\right)=\left(0.14,0.06,0.05,0.08\right)$
for SpanReg, respectively, indicating a roughly 60$\%$ improvement in
stability using SpanReg as compared to the DP.

###  3.4 Myelin Water Fraction Mapping of the Human Brain

Having shown the effectiveness of SpanReg on a wide range of simulated data,
we now demonstrate its application to brain imaging. In particular, we compare
the results of SpanReg, NNLS regularization using the DP, and non-regularized
NNLS for conventional myelin water fraction (MWF) mapping in the human brain.
In brief, this method involves acquiring a decaying bi-exponential signal for
each image pixel, from which a DF of spin-lattice relaxation times $T_{2}$ is
recovered. Given the high SNR of these images, we restrict our attention to
the Gaussian noise approximation. Note again that the simulations above were
targeted to parameter ranges appropriate for this problem of in vivo myelin
mapping.

3D gradient and spin-echo (GRASE) images were obtained from the brain of a
healthy 49-year-old female using $32$ echoes at $TE_{i}=i\times TE$, where
$TE=11.3$ ms for $i=1,2,\cdots,32$. The notation $TE$, standing for echo time,
is conventional in magnetic resonance, and corresponds to the measurement
times in Eq. (1.4). The acquisition sequence had additional parameters of
$TR=1000$ ms, echo planar imaging acceleration factor of 3, field of view
$278\text{ mm}\times 200\text{ mm}\times 30\text{ mm}$, acquisition matrix of
size $185\times 133\times 10$, acquisition voxel size $=1.5\text{ mm}\times
1.5\text{ mm}\times 3\text{ mm}$, reconstructed to $=1\text{ mm}\times 1\text{
mm}\times 3\text{ mm}$ using zero filling in k-space. Scan time was
approximately 10 minutes. A 3T Philips MRI system (Achieva, Best, the
Netherlands), equipped with an internal quadrature body coil for transmission
and an eight-channel phased-array head coil for reception, was used for
acquisition.

We employ both unfiltered and NESMA-filtered datasets. NESMA [7] is a highly
effective nonlocal denoising image filter that has been shown to increase the
quality of myelin water fraction (MWF) mapping. The SNR of the unfiltered
dataset is $\sim 10-300$, as defined by Eq. (1.6), with variation across
pixels due primarily to the variation in image amplitude.

While the ground truth for brain MWF in vivo is unknown, we can still
implement numerical experiments using this data to compare different
reconstruction approaches. First, we defined the MWF as the integral of the
$T_{2}$ DF between the initial abscissa values of $T_{2}=6$ ms and $T_{2}=40$
ms [6, 7] for each pixel and created the corresponding map. To compare two
different MWF maps $\mathbf{A}$ and $\mathbf{B}$, we define the scaled
absolute difference (SAD) by:

$SAD(\mathbf{A},\mathbf{B})=\frac{S_{\mathbf{A}-\mathbf{B}}}{S_{\mathbf{A}}}$
(3.6)

where $S_{\mathbf{M}}=\sum_{i,j}\left|M_{i,j}\right|$ is the sum of the
absolute values of the elements of a matrix $M_{i,j}$. In this way, if
$SAD(\mathbf{A},\mathbf{B}_{1})$, as a comparison between a reference map
$\mathbf{A}$ and an approximation $\mathbf{B}_{1}$, is smaller than
$SAD(\mathbf{A},\mathbf{B}_{2})$, we say that the reconstructed MWF map
$\mathbf{B}_{1}$ exhibits superior accuracy to $\mathbf{B}_{2}$.

Because absolute magnetic resonance signal amplitudes are arbitrarily scaled
by various sample characteristics, experimental settings and instrument
properties, it is necessary to normalize the data so that $\int
f(T_{2})dT_{2}=1$ before applying SpanReg. This represents normalization of
the PDF represented by the DF $f(T_{2})$. Then, from Eq. (1.4) written in
terms of $TE$, the signal must be normalized to its initial value, that is
$y(TE=0)=1$. This is problematic for two reasons. First, for experimental
reasons, the spin-echo experiment in particular has measurement times at
$i\times TE$ with $i=1$ rather than $i=0$. This limitation does not apply to
other similarly-modeled experiments, such as diffusion-sensitizing pulse
sequences. In addition, however, in any experiment, all data points are
corrupted by noise so that normalization of the acquired signal by its initial
value is not equivalent to normalizing the underlying signal, which is what is
desired; the underlying and observed signals differ by noise at each data
point. Instead, we proceed by setting the normalization factor to
$\left\lVert\mathbf{f}_{LS}\right\rVert_{1}$ according to the following:

$\begin{split}y_{LS}\left(TE=0\right)&\approx\left(e^{-0/T_{2,1}},e^{-0/T_{2,2}},\cdots,e^{-0/T_{2,n}}\right)\left(f_{LS,1},f_{LS,2},\cdots,f_{LS,n}\right)^{T}\Delta
T_{2}\\\ &=\sum_{j=1}^{n}f_{LS,j}\Delta
T_{2}=\left\lVert\mathbf{f}_{LS}\right\rVert_{1}\Delta T_{2}\end{split}$ (3.7)

where the NNLS estimate of $\mathbf{f}_{LS}\geq\mathbf{0}$ is determined from
non-normalized data according to Eq. (1.7).

We note that in effect, Eq. (1.7) is equivalent to Eq. (1.5), which is the
original problem. In the current context, we are not representing the solution
to Eq. (1.7) to be an accurate estimate of the entire DF $\mathbf{f}$, but
rather as the best available estimate of $f_{TE=0}$. This approach introduces
no additional bias though use of regularization. In our imaging analyses, the
observed data $\mathbf{y}_{\text{ob}}$ is divided by
$\left\lVert\mathbf{f}_{LS}\right\rVert_{1}$ for each pixel. This analysis was
not required for the simulation experiments above, where we illustrated
SpanReg in the more usual case of signals with an initial abscissa value of
zero. In Sections 3.4.1 and 3.4.2, we implement this approach to reconstruct
SpanReg-based MWF maps.

We present two distinct analyses for the brain dataset. First, in section
3.4.1, we define a high-quality reference MWF map, based on pixel-wise
reconstructions from SpanReg on NESMA-filtered data. This serves as a standard
of comparison for three reconstruction methods, SpanReg and NNLS with DP
regulariztion and without regulariztion, for two different levels of SNR.
Then, in the following section 3.4.2, we work directly with NESMA-filtered
(high SNR) and non-filtered (low SNR) images using all three reconstruction
methods to demonstrate their performance with respect to SNR; a comparison is
made between results obtained at high and low SNR for each method separately.

####  3.4.1 Comparison of Reconstruction Methods with Respect to a Reference
Image

We obtained the reference MWF map shown in the uppermost panel of Fig. 5 as
follows. First, we obtained a stack of successively $T_{2}$-weighted images
from the data acquisition sequence as described above; for each pixel, there
is therefore a corresponding decay curve. However, before processing these
curves as described in detail above, we first apply the NESMA filter to the
stack of images, resulting in pixel-wise decay curves of greatly improved SNR.
These are then inverted using SpanReg to obtain a $T_{2}$ distribution for
each pixel. The reference MWF map is formed from these distribution as
described above, by taking the integral of the DF’s up to $T_{2}$=40 ms for
each pixel separately.

The MWF maps created from high and low SNR data for comparison, shown in the
same figure, were generated as follows. The high-quality $T_{2}$ distribution
functions used to define the reference map were used to generate noiseless
decay curves, to which Gaussian noise was added to achieve SNR$=800$ (high
SNR) and SNR$=200$ (low SNR) in the decays. The decay curve at each pixel to
corresponds to $\mathbf{y}_{\text{ob}}$ in Eq. (1.5). We then reconstructed
the $T_{2}$ DF’s using SpanReg, and NNLS with DP regularization and without
regularization. MWF maps were constructed from these and compared to the
reference map using the SAD metric, Eq. (3.6).

Visual inspection shows that the SpanReg analysis yields a result more closely
resembling the reference image than either of the other two methods.

(a)  (b)

Figure 5: Uppermost image: reference MWF map; MWF maps generated from high SNR
(top row; SNR=800) and low SNR (bottom row; SNR=200) data using SpanReg, NNLS
with the DP, and non-regularized NNLS.

Quantitatively, the scaled absolute difference (SAD) between the reference map
and the derived MWF maps from the high SNR dataset are: 0.32 (compared to
SpanReg), 0.38 (DP), and 0.59 (NNLS). The corresponding values for the low SNR
dataset are: 0.39 (SpanReg), 0.58 (NNLS with DP regularization), and 0.68
(non-regularized NNLS). Thus, by this metric, SpanReg outperforms DP for the
high SNR and low SNR datasets by $16\%$ resp. $33\%$, and NNLS by $46\%$ resp.
$43\%$.

####  3.4.2 Comparison of Performance for High Versus Low SNR Data

We now examine the sensitivity to noise of the reconstruction methods. Even
with spatially invariant noise, SNR varies across an image because of signal
strength variations. Ideally the entire SpanReg algorithm would be applied
separately for each pixel according to its SNR. This is not practical for
typical images; in our case, the images consist of $\sim 25,000$ pixels.
Therefore, we divided the full range of image SNR into 18 bins over the range
of 10 to 800, encompassing both unfiltered and NESMA-filtered pixel-wise decay
signals. The same Gaussian basis functions were used in all cases, but
different SpanReg parameter values and functions
$\left\\{\beta_{i,\lambda_{j}}\right\\}$,
$\left\\{\mathbf{g}_{i,\lambda_{j}}\right\\}$ were obtained for each of these
18 SNR bins based on its SNR. The corresponding
$\left\\{\beta_{i,\lambda_{j}}\right\\}$ and
$\left\\{\mathbf{g}_{i,\lambda_{j}}\right\\}$ were used to reconstruct the
data from each pixel with SNR within that bin. MWF for each pixel was then
derived from the resulting distributions as described above, by integrating
the DF’s up to 40 ms. The resulting MWF maps are displayed in the
Supplementary Information, Fig. 7.

We first note that the results for SpanReg are visually more similar at the
two noise levels than are the results for the other reconstruction methods,
indicating greater stability with respect to noise. In addition, the MWF
distributions recovered by DP and LS show prominent edge effects surrounding
the ventricles, as well as myelin voids towards the brain periphery. In
contrast, the MWF maps obtained with SpanReg show a much more qualitatively
normative pattern. Quantitatively, the SAD between the high and low SNR MWF
maps, $SAD\left(\mathbf{M}_{\text{high SNR}},\mathbf{M}_{\text{low
SNR}}\right)$, for the three reconstruction methods are 0.35 for SpanReg, 0.43
for DP, and 0.63 for NNLS. This is consistent with the results shown in
Section 3.4.1.

##  4 Discussion

In this paper, we present a new approach, SpanReg, to obtaining a regularized
solution of the discretized Fredholm equation of the first kind. This
technique should be widely applicable for inverse problems requiring
regularization to obtain meaninful solutions. The idea of the method is to
incorporate information content over a range of regularized solutions rather
than to attempt to identify a single optimal regularization parameter,
$\lambda$. For inverse problems of a type arising extensively in MR
relaxometry, we provided simulations demonstrating greater accuracy of SpanReg
as compared to conventional Tikhonov regularization using the DP for selection
of $\lambda$. We also showed that SpanReg is more stable with respect to the
choice of regularization parameters, in the sense that the departure from an
optimal solution is smaller when the solution arrived at using the optimal
$\lambda$ is not incorporated. The method was also demonstrated on
experimental biomedical MRR data to evaluate the MWF in the human brain; this
is a very challenging problem of great current interest.

The framework put forward in this work is to provide a more general form of
the recovered DF: $\mathbf{f}^{*}=\sum_{j}\alpha_{j}\mathbf{f}_{\lambda_{j}}$;
conventional Tikhonov regularization incorporates just a single term, optimal
in some sense, of this sum. In addition to relaxing the restriction of
isolating the single best regularized solution, our method also provides a
flexible form for the construction of the desired solution. To obtain the
coefficients $\left\\{\alpha_{j}\right\\}$, we formulated the solution
$\sum_{j}\alpha_{j}\mathbf{f}_{\lambda_{j}}$ as a linear combination of
multiple regularized solutions via an over-complete dictionary of Gaussian
functions $\left\\{\mathbf{g}_{i,\lambda_{j}}\right\\}$. The algorithm can be
separated into offline and online computations, where the offline part is
specific to a given SNR and basis set, but is independent of a specific
instance of experimental data. The offline calculation requires the
determination of the noise-corrupted basis functions and parameters, e.g.
fixing the number $M$ of Gaussian dictionary functions $\mathbf{g}_{i}$, the
number $N$ of regularization parameters used by the algorithm, and performing
the actual computation required to obtain the noise-corrupted DF,
$\left\\{\mathbf{g}_{i,\lambda_{j}}\right\\}$.The online calculation is
dependent upon a specific data set and is much less computation-intensive.

This formulation exhibits a condition number $cond(\mathbf{B})=O(10^{5})$ in
Eq. (2.13) for the simulation analysis in Section 3.2, compared to the
condition number of $cond(\mathbf{A})$ in the NNLS problem, which is
$O(10^{20})$.

In this study, we demonstrated that the SpanReg method achieves higher
accuracy than the conventional DP parameter selection method, which also
requires knowledge of noise level in the observation $\mathbf{y}_{\text{ob}}$.
For the prototypical bimodal DF’s studied, SpanReg exhibits a substantially
greater ability to resolve the two components. With more disparate spaced
components, SpanReg also achieved improved reconstruction of positions and
amplitudes of the components. Further, SpanReg is not constrained to the
selection of a single optimal value of $\lambda$; in fact, in conventional
practice, this optimal value will depend upon the specific selection criterion
implemented. Thus, application of different conventional methods is virtually
guaranteed to recover different DF’s. In contrast, our more general
formulation forms the reconstruction through a linear combination of
differently-weighted solutions with appropriate weights, and so is less
dependent upon the exact degree of regularization of any one of these.

Previous studies have also presented important reconstruction methods based on
Tikhonov regularization with multiple regularization parameters. In distinct
contrast to our work,these have employing a modified regularization penalty
from which is derived a single optimal solution. These methods include
multiple regularization [4, 29], multi-parameter regularization [10, 25, 37,
14, 11], and multi-penalty regularization [13, 1]. This literature generally
studies the same optimization problem as we have worked with, with the
somewhat more general form:

$\mathbf{f}^{*}=\mathop{\rm
argmin}_{\mathbf{f}}\left\\{\left\lVert\mathbf{y}_{\text{ob}}-\mathbf{A}\mathbf{f}\right\rVert_{2}^{2}+\sum_{j=1}^{N}\lambda_{j}\mathcal{L}_{j}\left(\mathbf{f}\right)\right\\}$
(4.1)

where $\mathcal{L}_{j}$ are regularization operators and $\lambda_{j}$ are
regularization parameters associated with the corresponding operator. One of
the most well-known of these is elastic-net regularization (EN) [38], where
$\mathcal{L}_{1}\left(\mathbf{f}\right)=\left\lVert\mathbf{f}\right\rVert_{1}$,
and
$\mathcal{L}_{2}\left(\mathbf{f}\right)=\left\lVert\mathbf{f}\right\rVert_{2}^{2}$.
The EN is able to promote sparsity while also being able to handle smoother
distributions. Other types of penalty terms have also been proposed involving
smoothness constraints; these include penalties of the form
$\left\lVert\nabla\mathbf{f}\right\rVert_{1}$ [4, 25, 37, 29, 14, 11],
$\left\lVert\nabla\mathbf{f}\right\rVert_{2}^{2}$ [25, 13, 14], and other
penalties of the form $\left\lVert\mathbf{W}\mathbf{f}\right\rVert$ [4, 10,
29, 11], where $\mathbf{W}$ may be, for example a projection or transformation
matrix. Unlike SpanReg, all of these methods seek optimal regularization
parameters for a specific form of regularization rather than incorporating
differently-regularized solutions.

An interesting alternative method more closely related to ours deviates from
the formalism of Eq. (4.1) by constructing a linear combination of regularized
solutions [22], as in our treatment. The components forming the recovered DF
are obtained from different regularization methods, with, as in conventional
treatments and distinct from SpanReg, a single optimal regularizer being
identified for each. Nevertheless, this work introduces the idea of combining
DF’s, each of which in effect reflects a potentially desirable quality, into a
final derived DF.

A variety of parameter selection techniques have been presented for the
conventional types of multi-parameter regularization described above,
generally based on the classical approaches for single-parameter
regularization optimization. These include the multi-parameter discrepancy
principle [25], L-hypersurface [3, 4], and GCV-multi [11]. Other methods [10]
such as simple grid search [32] and the balancing principle [1] have been
implemented. However, as emphasized above, none of these studies presents the
notion of including sub-optimal regularized solutions as we have proposed and
demonstrated with SpanReg. Thus, the current work is fundamentally different
from what has appeared in the literature in; we form a linear combination of
sub-optimal solutions, in contrast to previous work which has used
combinations of regularization terms leading to a single “optimal” solution,
or [22] a linear combination of optimal solutions. SpanReg allows us to
combine the stability of more highly-regularized solutions with the resolution
of less-regularized solutions; the potential for the efficacy of this
construction was provided by our observations of the different responses to
regularization of very similar ill-posed problems.

SpanReg should be widely applicable to a wide range of linear inverse problems
based on the Fredholm equation of the first kind. As one class of examples, we
demonstrated its applicability to MRR investigations of brain, where the DF is
based on the $T_{2}$ distribution of tissue. However, the same formulation
also applies to diffusion MRI [5], where application is simpler due to the
availability of the $b=0$ data point, where $b$ is the conventional symbol for
the combination of gradient strengths and timing that defines diffusion
sensitization. In this case, the formalism of Eq. (3.7) is not required.
Similar comments apply to the closely-related problem of $T_{1}$ MRR [23].
Further, MRR applications extend far beyond those in biomedicine, and include
the food sciences [24], engineering [16], and the petrochemical industry [33].
There are also many applications of SpanReg outside of MR studies, including
fluorescence analysis [31].

Although we have presented a novel and effective method for combining multiple
degrees of regularization into an ill-conditioned problem, certain limitations
remain. We have not fully explored selection of an optimal basis, although the
Gaussian basis set used has been shown to perform well for reconstructing a
range of challenging DF’s. In addition, there are a number of user-selected
parameters for SpanReg, and we have not yet explored these systematically.
These include the number $N$ of regularized solutions to incorporate into the
DF as well as their associated regularization parameter values $\lambda_{j}$;
the size $M$ and details of the Gaussian dictionary used to represent the DF
for each of these solutions; and the number $n_{\text{run}}$ of noise
realizations used for averaging to render the results robust with respect to
noise realization. In this paper, our considerations have been based on the
goal of maintaining reasonable conditioning of the the LS problem in Eq.
(2.6); as demonstrated, this empirical approach has worked well. We also
emphasize that the requirement for knowledge of SNR is a limitation of SpanReg
in many contexts, although for the large category of problems based on
decaying signals, including MRR and related experiments in MR, SNR estimates
are generally available. This consideration also applies to certain
conventional parameter selection methods such as DP. In spite of these
limitation, we have shown that SpanReg can be applied to datasets with
varying, though known, levels of SNR, such as in the brain MWF analysis. This
does require a lengthier offline computation, with separate calculations
required across a range of SNR values. The degree to which a given range of
SNR values must be discretized also remains an open question, though the
success of the binning procedure described in the brain MRI analysis support
the notion that this discretization need not be unduly fine. This is a
critical finding, indicating that the noise estimate required for SpanReg can
be approximate. The degree to which this holds remains to be quantified.

In addition to theoretical developments and demonstration on simulated data,
we have provided an analysis of in vivo brain MRI data. We compare the
performance of SpanReg to two other methods for recovering DF’s from decay
curves for application to brain MWF mapping. This presentation is not intended
to provide an extensive comparison of SpanReg with other state-of-the-art
methods for assessing MWF; development of these methods is a topic of major
current interest, with rapid introduction of new approaches [26, 8, 7]. A full
comparison with these is beyond the scope of the present manuscript, which is
to introduce and provide an initial application of SpanReg.

In addition to addressing the current limitations of SpanReg as outlined
above, potential extensions include exploration of alternative forms of
Tikhonov and other forms of regularization penalties. In addition, SpanReg can
be extended in a straightforward fashion to 2D and higher dimensional MRR
analysis [27], for which the governing inverse problem is a higher-dimensional
inverse Laplace transform.

In conclusion, we have proposed a new approach to determining regularized
solutions to the Fredholm integral equation of the first kind by incorporating
the information content of non-optimal solutions, and have demonstrated its
efficacy in simulations and through application to MRR of the human brain.
SpanReg should be widely applicable throughout the field of inverse problems,
presenting an alternative to Tikhonov and related forms of regularization

## Acknowledgments

This work was supported by the Intramural Research Program of the National
Institute on Aging of the National Institutes of Health. We thank Pak-Wing Fok
for particularly helpful discussions.

## References

  * [1] Abhishake and S. Sivananthan.  Multi-penalty regularization in learning theory.  Journal of Complexity, 36:141–165, 2016.
  * [2] Frank Bauer and Mark A. Lukas.  Comparing parameter choice methods for regularization of ill-posed problems.  Mathematics and Computers in Simulation, 81(9):1795–1841, 2011\.
  * [3] Murat Belge, Misha E. Kilmer, and Eric L. Miller.  Simultaneous multiple regularization parameter selection by means of the L-hypersurface with applications to linear inverse problems posed in the wavelet transform domain.  In Ali Mohammad-Djafari, editor, Bayesian Inference for Inverse Problems, volume 3459, pages 328 – 336. International Society for Optics and Photonics, SPIE, 1998.
  * [4] Murat Belge, Misha E Kilmer, and Eric L Miller.  Efficient determination of multiple regularization parameters in a generalized L-curve framework.  Inverse Prblems, 18(4):1161–1183, 2002.
  * [5] Dan Benjamini, Michal E. Komlosh, and Peter J. Basser.  Imaging local diffusive dynamics using diffusion exchange spectroscopy MRI.  Physical Review Letters, 118(15):158003, 2017.
  * [6] Jean-Marie Bonny, Amidou Traore, Mustapha Bouhrara, Richard G. Spencer, and Guilhem Pages.  Parsimonious discretization for characterizing multi-exponential decay in magnetic resonance.  NMR in Biomedicine, 33(12):e4366, 2020.
  * [7] Mustapha Bouhrara, David A. Reiter, Michael C. Maring, Jean-Marie Bonny, and Richard G. Spencer.  Use of the nesma filter to improve myelin water fraction mapping with brain mri.  Journal of Neuroimaging, 28(6):640–649, 2018.
  * [8] Mustapha Bouhrara and Richard G. Spencer.  Rapid simultaneous high-resolution mapping of myelin water fraction and relaxation times in human brain using bmc-mcdespot.  NeuroImage, 147:800–811, 2017.
  * [9] D. Calvetti, S. Morigi, L. Reichel, and F. Sgallari.  Tikhonov regularization and the L-curve for large discrete ill-posed problems.  Journal of Computational and Applied Mathematics, 123(1):423–446, 2000.
  * [10] Zhongying Chen, Yao Lu, Yuesheng Xu, and Hongqi Yang.  Multi-parameter Tikhonov regularization for linear ill-posed operator equations.  Journal of Computational Mathematics, 26(1):37–55, 2008.
  * [11] Julianne Chung and Malena I Español.  Learning regularization parameters for general-form Tikhonov.  Inverse Prblems, 33(7):074004, 2017.
  * [12] Christophe Couvreur and Yoram Bresler.  Dictionary-based decomposition of linear mixtures of Gaussian processes.  In 1996 IEEE International Conference on Acoustics, Speech, and Signal Processing Conference Proceedings, volume 5, pages 2519–2522. IEEE, 1996\.
  * [13] Massimo Fornasier, Valeriya Naumova, and Sergei V. Pereverzyev.  Parameter choice strategies for multipenalty regularization.  SIAM Journal on Numerical Analysis, 52(4):1770–1794, 2014.
  * [14] Silvia Gazzola and Lothar Reichel.  A new framework for multi-parameter regularization.  BIT Numerical Mathematics, 56(3):919–949, 2016.
  * [15] Gene H Golub, Michael Heath, and Grace Wahba.  Generalized cross-validation as a method for choosing a good ridge parameter.  Technometrics, 21(2):215–223, 1979.
  * [16] Jessie Greener, Hartwig Peemoeller, Changho Choi, Rick Holly, Eric J. Reardon, Carolyn M. Hansson, and Mik M. Pintar.  Monitoring of hydration of white cement paste with proton NMR spin–spin relaxation.  Journal of the American Ceramic Society, 83(3):623–627, 2000.
  * [17] So Hyan Han, E Ackerstaff, Radka Stoyanova, S Carlin, Wei Huang, JA Koutcher, JK Kim, G Cho, G Jang, and Hyungjoon Cho.  Gaussian mixture model-based classification of dynamic contrast enhanced MRI data for identifying diverse tumor microenvironments: preliminary results.  NMR in Biomedicine, 26(5):519–532, 2013.
  * [18] Per Hansen and Dianne O’leary.  The Use of the L-curve in the Regularization of Discrete Ill-Posed Problems.  SIAM J. Sci. Comput., 14:1487–1503, November 1993.
  * [19] Per Christian Hansen.  Rank-Deficient and Discrete Ill-Posed Problems: Numerical Aspects of Linear Inversion.  Society for Industrial and Applied Mathematics, Philadelphia, January 1987\.
  * [20] Per Christian. Hansen.  Analysis of Discrete Ill-Posed Problems by Means of the L-Curve.  SIAM Review, 34(4):561–580, December 1992.
  * [21] Per Christian Hansen.  Discrete inverse problems: insight and algorithms.  SIAM, 2010.
  * [22] Michiel E. Hochstenbach and Lothar Reichel.  Combining approximate solutions for linear discrete ill-posed problems.  Journal of Computational and Applied Mathematics, 236(8):2179–2185, 2012.
  * [23] Martin D Hürlimann, Lauren Burcaw, and Yi-Qiao Song.  Quantitative characterization of food products by two-dimensional D–T2 and T1–T2 distribution functions in a static gradient.  Journal of colloid and interface science, 297(1):303–311, 2006\.
  * [24] Emrah Kirtil and Mecit H Oztop.  1H nuclear magnetic resonance relaxometry and magnetic resonance imaging and applications in food science and processing.  Food Engineering Reviews, 8(1):1–22, 2016.
  * [25] Shuai Lu and Sergei V Pereverzev.  Multi-parameter regularization and its numerical realization.  Numerische Mathematik, 118(1):1–31, 2011.
  * [26] Alex L MacKay and Cornelia Laule.  Magnetic resonance of myelin water: an in vivo marker for myelin.  Brain plasticity, 2(1):71–91, 2016.
  * [27] J. Mitchell, T.C. Chandrasekera, and L.F. Gladden.  Numerical estimation of relaxation and diffusion distributions in two dimensions.  Progress in Nuclear Magnetic Resonance Spectroscopy, 62:34–50, 2012\.
  * [28] Vladimir Alekseevich Morozov.  Methods for solving incorrectly posed problems.  Springer Science & Business Media, 2012.
  * [29] Varun P. Gopi, P. Palanisamy, Khan A. Wahid, Paul Babyn, and David Cooper.  Multiple regularization based MRI reconstruction.  Signal Processing, 103:103–113, 2014.
  * [30] Ashish Raj, Sneha Pandya, Xiaobo Shen, Eve LoCastro, Thanh D Nguyen, and Susan A Gauthier.  Multi-compartment T2 relaxometry using a spatially constrained multi-Gaussian model.  PLoS One, 9(6):e98391, 2014.
  * [31] Matthias Reisser, Johannes Hettich, Timo Kuhn, Achim P. Popp, Andreas Große-Berkenbusch, and J. Christof M. Gebhardt.  Inferring quantity and qualities of superimposed reaction rates from single molecule survival time distributions.  Scientific Reports, 10(1):1758, Feb 2020.
  * [32] Christiana Sabett, Ariel Hafftka, Kyle Sexton, and Richard G Spencer.  L1, Lp, L2, and elastic net penalties for regularization of Gaussian component distributions in magnetic resonance relaxometry.  Concepts in Magnetic Resonance Part A, 46(2):e21427, 2017.
  * [33] Yi-Qiao Song and Ravinath Kausik.  NMR application in unconventional shale reservoirs – A new porous media research frontier.  Progress in Nuclear Magnetic Resonance Spectroscopy, 112-113:17–33, 2019.
  * [34] Richard G. Spencer and Chuan Bi.  A tutorial introduction to inverse problems in magnetic resonance.  NMR in Biomedicine, 33(12):e4315, 2020.
  * [35] Grace Wahba.  Practical approximate solutions to linear operator equations when the data are noisy.  SIAM Journal on Numerical Analysis, 14(4):651–667, 1977.
  * [36] Grace Wahba.  Spline models for observational data, volume 59.  SIAM, 1990.
  * [37] Zewen Wang.  Multi-parameter Tikhonov regularization and model function approach to the damped Morozov principle for choosing regularization parameters.  Journal of Computational and Applied Mathematics, 236(7):1815–1832, 2012.
  * [38] Hui Zou and Trevor Hastie.  Regularization and variable selection via the elastic net.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 67(2):301–320, 2005.

##  5 Supplementary information

(a)  (b)

(c)  (d)

Figure 6: Stability with respect to choices of $\lambda$’s, using a single
value of $\lambda$ for Tikhonov regularization and triplets of $\lambda$’s for
SpanReg. In all cases, SNR=500. The true distribution is
$\frac{2}{3}g(30,80)+\frac{1}{3}g(3,5)$. The optimal solution is
$f_{\lambda_{3}}$. Upper left: Solutions obtained with optimal and sub-optimal
regularizations. Upper right: Solutions obtained with SpanReg with
subsequences of regularized solutions. Bottom left: Departures from optimal
regularized solution in the upper left panel. Bottom right: Departures from
corresponding SpanReg solution in the upper right panel. Figure 7: MWF maps
reconstructed from original imaging data, to which two different levels of
noise were added.