# Inversion of the Laplace Transform of Point Masses

Michael Mckenna1, Hrushikesh N. Mhaskar2 , Richard G. Spencer1 Institute of
Mathematical Sciences, Claremont Graduate University, Claremont, CA 91711
(United States).The research of HNM is supported in part by NSF grant DMS
2012355 and ONR grants N00014-23-1-2394, N00014-23-1-2790. The research of MM
and RGS was funded in part by the National Institute on Aging, Intramural
Research Program. email<EMAIL_ADDRESS>
(1 National Institute on Aging, National Institutes of Health, Laboratory of
Clinical Investigation, Magnetic Resonance Imaging and Spectroscopy Section,
Baltimore MD 21224
2Institute of Mathematical Sciences, Claremont Graduate University, Claremont,
CA 91711 U.S.A.
)

###### Abstract

Motivated by applications in magnetic resonance relaxometry, we consider the
following problem: Given samples of a function
$t\mapsto\sum_{k=1}^{K}A_{k}\exp(-t\lambda_{k})$, where $K\geq 2$ is an
integer, $A_{k}\in\mathbb{R}$, $\lambda_{k}>0$ for $k=1,\cdots,K$, determine
$K$, $A_{k}$’s and $\lambda_{k}$’s. Our approach is to transform this function
into another function of the same form where $\lambda_{k}$’s are replaced by
$i\lambda_{k}$. For this purpose, we study the least square approximation
using polynomials weighted by the Gaussian weight, and use the fact that
Hermite functions are eigenfunctions of the Fourier transform. We provide a
detailed analysis of the effect of noise in the data.

##  1 Introduction

Multiexponential decay signals, consisting of the superposition of two or more
exponentials of the form:

$F_{0}(t)=\sum_{k=1}^{K}A_{k}\exp(-t\lambda_{k})+\epsilon_{0}(t),\qquad 0\leq
t\leq T,$ (1.1)

where $T>0$, $K$ is a positive integer, $A_{k},\lambda_{k}>0$ for
$k=1,\cdots,K$, and $\epsilon_{0}(t)$’s are i.i.d. random variables with mean
$0$, arise in innumerable settings in the natural sciences and mathematics
[34], including optics [18], astrophysics [56], food science [37], geophysics
[59], and biology [17, 47]. The main problem is to determine the parameters
$a_{k}$ and $\lambda_{k}$. The literature is less focused on determination of
the number of components $K$, although that may also be of substantial
interest and can be addressed by a number of statistical methods. The value of
$K$ is often assumed.

Our main motivating example is magnetic resonance relaxometry, in which the
proton nuclei of water are first excited with radio frequency pulses and then
exhibit an exponentially decaying electromagnetic signal. These applications
most often have $K=2$, although applications with $K\geq 3$ appear
occasionally. In general, biological tissue will contain multiple states of
water defined by, among other factors, molecular mobility, often modeled as
the co-existence of discrete compartments [29]. In principle, these
compartments can then be separately mapped based on the fact that they exhibit
different MR transverse decay time constants, or $T_{2}$ values, as follows. A
particularly important example is myelin mapping in the brain. Myelin forms an
insulating sheath surrounding neuronal axons in the brain and spinal cord and
is critical for the transmission of electrical signals [35]. Water within the
myelin sheath itself exhibits limited mobility and therefore, according to MR
theory, exhibits a more rapid electromagnetic signal decay as compared to the
more mobile water. The more rapidly decying signal component is designated as
the myelin water fraction (MWF). Thus, the MR signal is composed of sub-
signals from both the MWF and the more mobile water, with these two
contributions separately quantified by making use of their different
relaxation time constants. MWF mapping in the human brain based on this
principle was pioneered by Mackay, [41], who focused on the related model of a
discretized continuous range of tissue decay time constants. MWF mapping is
applied widely in the investigation of cognitive decline [16, 28], Alzheimer’s
disease [23, 36], neurodegenerative diseases such as multiple sclerosis [39],
and in studies of normative aging that are characterized by or exhibit myelin
loss [10, 11, 12, 13, 16, 24, 39, 51, 64].

In addition to widespread application, these techniques are still undergoing
important refinement, due to the complexity of separating the signal generated
by the MWF from the signal generated by the more mobile water The difficulty
in this and related multiexponential analyses may be understood in the
following way. It is well-known that the solution of the Fredholm equation of
the first kind (FEFK)

$g(t)=\int_{a}^{b}K(t,\tau)f(\tau)d\tau$ (1.2)

is an ill-posed inverse problem, with determination of the unknown function,
$f$, multiplying the Fredholm kernel K, being highly sensitive to noise. When
the kernel function is the exponential, $K(t,\tau)=exp(-t/\tau)$, parameter
extraction is a special case of the notoriously ill-posed inverse Laplace
transform (ILT)[19]. If the support of the unknown distribution function is
the sum of point masses, we recover the multiexponential decay problem (1.1),
which inherits its ill-posedness from the FEFK. In the discretized case, these
difficulties are highlighted by the condition number of the matrix defining
the discretization of the FEFK, which can become enormous with fine
discretization, even for only two monoexponential components [8, 60]. This
effect of increasingly finer discretization may be understood by the fact that
it brings the discretized problem closer to the ill-posed native integral
equation of the FEFK.

This behavior is in marked contrast to that of the Fourier transform; as a
unitary transform, with appropriate normalization, the FT enjoys the optimal
condition number of unity; that is, $\|\mathfrak{F}^{*}\|\|\mathfrak{F}\|=1$,
where $\|\cdot\|$ is the $L^{2}$ operator norm [60]. Given this, the
transformation of the ILT to an equivalent Fourier problem is an appealing
approach, as was first (to our knowledge) described by Gardner. This allowed
for the visualization of decay constants as peaks on a frequency spectrum
[27]. This work was later built upon by Provencher through use of a of tapered
window function to reduce sidelobes in the frequency domain defining the rates
$\lambda_{k}$ [54].

Here we present a further development in this area, based on the properties of
the Hermite functions. Our approach is in some sense reminiscent of previous
work on the expansion of multiexponential signals by orthogonal polynomials
[33], including an earlier approach for inversion of the Laplace transform
through the use of the Jacobi and Laguerre polynomials. [49]. In our case, we
find that the fact that the Hermite functions are eigenfunctions of the
Fourier transform permits a particularly natural transition between the
original problem domain and the Fourier domain. Furthermore, unlike many
previous mathematical treatments of multi-component multiexponential decay
signals, we incorporate experimental noise into our analysis.

The organization of this paper is as follows. We describe some related work in
Section 2. Section 3 presents definitions and notation. Section 4 introduces
our main idea, which is expansion of the multiexponential function in terms of
Hermite polynomials and transformation to Fourier space. Section 5 presents a
general theory for computation of the coefficients in the expansion of Hermite
polynomials from discrete values of $F_{0}$. We describe the process of least
square approximation and related quantities in Section 6. Based on these
considerations, Section 7 describes our main theorems concerning the errors in
the various approximations, including errors due to noise. The theorems in
Section 7 are proved in Section 9. The procedure for obtaining $\lambda_{k}$’s
from the Fourier transform is based on our earlier work, and is summarized in
the Appendix. A numerical implementation is given in Section 8.

##  2 Related work

Many methods have been proposed, developed, and implemented over the past 200
years for exponential analysis [7, 18, 31, 33, 34, 50, 52, 53, 56]. Methods
range from the historical Prony method, still used in modified form, to neural
network approaches [32, 48, 55, 66]. Extensive review material may be found in
[33, 34]. Given theoretical bounds on the ability to define these parameters
as a function of signal-to-noise ratio (SNR), measurement intervals and total
duration, and number of underlying exponentials [3, 4, 6, 5], development of
further methodology represents attempts to extract the maximum information
from a given signal, or to reformulate the problem to circumvent these
limitations [7, 30].

We emphasize that in this work, we are analyzing the expression (1.1), which
has been applied to both myelin quantification and other tissue
characterization studies in MRI [14, 15, 20, 40, 41, 57, 58]. A related body
of published work in the magnetic resonance literature instead aims to recover
a spectrum of $T_{2}$ values; the underlying physics is similar for the two
approaches, although the mathematical analysis differs between them.

##  3 Notation

The orthonormalized Hermite functions are defined (cf. [61, Formula (5.5.3)])
by

$\psi_{k}(x)=(\pi^{1/2}2^{k}k!)^{-1/2}(-1)^{k}\exp(-x^{2}/2)\left(\frac{d}{dx}\right)^{k}\exp(-x^{2}),\qquad
k\in{\mathbb{Z}}_{+}.$ (3.1)

For $n>0$, we denote

$\Pi_{n}=\mathsf{span}\\{\psi_{k}:k<n\\}.$

Each element of $\Pi_{n}$ is a function of the form $x\mapsto
P(x)\exp(-x^{2}/2)$ for some polynomial $P$ of degree $<n$. Such a function
will be called a _weighted polynomial of degree $<n$_. We note that the
notation does not require $n$ to be an integer; we interpret
$\Pi_{n}=\Pi_{\lfloor n\rfloor}$ if $n$ is not an integer. This simplifies
some notation later.

The Fourier transform of a function $f:{\mathbb{R}}\to{\mathbb{R}}$ is defined
formally by

$\mathfrak{F}(f)(\omega)=\frac{1}{\sqrt{2\pi}}\int_{\mathbb{R}}f(t)\exp(-i\omega
t)dt.$ (3.2)

It is known (cf. [43, Theorem 4.3.2]; the notation is different) that

$\mathfrak{F}(\psi_{k})(\omega)=(-i)^{k}\psi_{k}(\omega),\qquad
k\in{\mathbb{Z}}_{+},\ \omega\in{\mathbb{R}}.$ (3.3)

Next, we introduce some concepts and notation related to various measures
which we will use in this paper.

Rather than writing sums of the form $\sum_{j=1}^{m}w_{j}g(y_{j})$, it is
convenient to use the Stieltjes’ integral notation. Let $\delta_{x}$ be the
Dirac delta at $x$. As a functional (or measure), this is expressed as $\int
gd\delta_{x}=g(x)$. We define a (possibly complex) measure
$\nu=\sum_{j=1}^{m}w_{j}\delta_{y_{j}}$ (described as a measure that
associates the mass $w_{j}$ with each $y_{j}$) and write $\int g(y)d\nu(y)$ or
even shorter, $\int gd\nu$, instead of $\sum_{j=1}^{m}w_{j}g(y_{j})$. The
_support_ of $\nu$ (denoted by $\mathsf{supp}(\nu)$) is
$\\{y_{j}\\}_{j=1}^{m}$, and $\nu$ is said to be _supported on_ any set that
contains the support. The quantity $\sum_{j=1}^{m}|w_{j}|$ is the _total
variation_ of $\nu$, denoted by $\|\nu\|_{TV}$. There is a measure associated
with $\nu$ called the total variation measure, denoted by $|\nu|$. For
example, for the above measure $\nu$,

$\int gd|\nu|=\sum_{j=1}^{m}|w_{j}|g(y_{j}).$

The Laplace transform $\mathcal{L}(\nu)$ and the Fourier transform
$\mathfrak{F}(\nu)$ for $\nu=\sum_{j=1}^{m}w_{j}\delta_{y_{j}}$ are defined by

$\mathcal{L}(\nu)(t)=\int\exp(-tu)d\nu(u)=\sum_{j=1}^{m}w_{j}\exp(-ty_{j}),\quad\mathfrak{F}(\nu)(\omega)=\sum_{j=1}^{m}w_{j}\exp(-i\omega
y_{j}).$ (3.4)

If $\nu$ is any possibly signed measure, and $A\subseteq{\mathbb{R}}$ is
$\nu$-measurable, we will write formally

$\|f\|_{\nu;p,A}=\begin{cases}\displaystyle\left\\{\int_{A}|f|^{p}d|\nu|\right\\}^{1/p},&\mbox{
if $1\leq p<\infty$,}\\\\[4.30554pt] \sup_{t\in\mathsf{supp}(\nu)\cap
A}|f(t)|,&\mbox{ if $p=\infty$.}\end{cases}$ (3.5)

The space $L^{p}(\nu)$ comprises measurable functions $f$ for which
$\|f\|_{\nu;p}<\infty$. If $A={\mathbb{R}}$, we will omit the mention of $A$
from the notation. When $\nu$ is the Lebesgue measure on ${\mathbb{R}}$, we
will omit its mention and write $\|f\|_{p}$ in place of $\|f\|_{\nu;p}$,
$L^{p}$ instead of $L^{p}(\nu)$, etc. As usual, we identify functions which
are equal $|\nu|$-almost everywhere.

For $f\in L^{p}$, we define the degree of approximation of $f$ by

$E_{p,n}(f)=\min_{P\in\Pi_{n}}\|f-P\|_{p}.$ (3.6)

Constant convention:

In the sequel, $c,c_{1},\cdots$ will denote generic positive constants, whose
values may be different at different occurrences, even within a single
formula. The notation $A\lesssim B$ means $A\leq cB$, $A\gtrsim B$ means
$B\lesssim A$, $A\sim B$ means $A\lesssim B\lesssim A$. For example, in a
formula of the form

$A_{n}\lesssim n^{c}\exp(-c_{1}n)\lesssim\exp(-cn),$

there are implied constants in $\lesssim$, and the constants $c,c_{1}$ in the
middle term are different from the constant $c$ in the last term.

Finally, we recall a fundamental result from [43, Theorem 6.2.4, Theorem
6.1.6].

###### Proposition 3.1

For any $\delta>0$, $n>0$, $1\leq p\leq\infty$, and $P\in\Pi_{n}$,

$\|P\|_{p,{\mathbb{R}}\setminus[-\sqrt{2n}(1+\delta),\sqrt{2n}(1+\delta)]}\leq\exp(-cn)\|P\|_{p,[-\sqrt{2n}(1+\delta),\sqrt{2n}(1+\delta)]}.$
(3.7)

In particular,

$\|P\|_{p}\sim\|P\|_{[-2\sqrt{n},2\sqrt{n}]}.$ (3.8)

##  4 Fundamental Idea

To explain our main idea, we first define a noiseless signal by (cf. (1.1))

$F(t)=\sum_{k=1}^{K}a_{k}\exp(-t\lambda_{k}),\qquad t\in{\mathbb{R}}.$ (4.1)

With the measure

$\mu^{*}(t)=\sum_{k=1}^{K}a_{k}\delta_{\lambda_{k}},$ (4.2)

where $\delta_{\lambda}$ denotes the Dirac delta supported at $\lambda$, we
observe that $f$ is the Laplace transform of $\mu^{*}$. The fundamental idea
is to convert this Laplace transform of $\mu^{*}$ into the Fourier transform
of another measure $\tilde{\mu}$ supported on $\\{\lambda_{k}\\}$ as in (4.5)
below. Towards this end, we introduce another function

$f(t)=F(t)\exp(-t^{2}/2)=\sum_{k=1}^{K}a_{k}\exp(\lambda_{k}^{2}/2)\exp(-(t+\lambda_{k})^{2}/2).$
(4.3)

and observe the fact that
$\mathfrak{F}(\exp(-\cdot^{2}/2))(\omega)=\exp(-\omega^{2}/2)$. We obtain:

$\mathfrak{F}(f)(-\omega)=\exp(-\omega^{2}/2)\sum_{k=1}^{K}a_{k}\exp(\lambda_{k}^{2}/2)\exp(-i\omega\lambda_{k}).$
(4.4)

With

$\tilde{\mu}(u)=\sum_{k=1}^{K}a_{k}\exp(\lambda_{k}^{2}/2)\delta_{\lambda_{k}},$
(4.5)

we observe

$\mathfrak{F}(\tilde{\mu})(\omega)=\exp(\omega^{2}/2)\mathfrak{F}(f)(-\omega).$
(4.6)

We will describe in the appendix a method to estimate the $\lambda_{k}$’s from
samples of $\mathfrak{F}(\tilde{\mu})(\omega)$ using our theory of localized
trigonometric polynomial kernels [21, 26, 45, 46].

In order to compute $\mathfrak{F}(\tilde{\mu})(\omega)$, our basic idea is to
expand $f$ in terms of the Hermite functions, and use (3.3).

There are three technical barriers here: one is that the signal $F_{0}$ as
defined in (1.1) is defined only on a finite interval. The second is that we
know only finitely many values of $F_{0}$, not a formula for $F_{0}$. Finally,
there is noise.

We will address the first barrier in this section, the second in Section 5,
and study the effect of noise in Section 7.

To remedy the barrier of knowing the signal only a finite interval, we will
shift the signal to the left, so that the resulting interval on which the
observations lie can contain all the quadrature points which are needed to
evaluate the integrals from the values of the functions. To carry out this
program, we introduce a tunable parameter $L>0$, write $R=T-L$, and define

$\displaystyle F_{\mbox{shifted}}(t)$
$\displaystyle=F_{0}(t+L)=\begin{cases}\sum_{k=1}^{K}A_{k}\exp(-L\lambda_{k})\exp(-t\lambda_{k})+\epsilon(t),&\mbox{
if $-L\leq t\leq R$,}\\\ 0&\mbox{ otherwise,}\end{cases}$ (4.7)
$\displaystyle=\begin{cases}F(t)+\epsilon(t),&\mbox{ if $-L\leq t\leq R$,}\\\
0&\mbox{ otherwise,}\end{cases}$

where

$a_{k}=A_{k}\exp(-L\lambda_{k}),\qquad\epsilon(t)=\epsilon_{0}(t+L).$ (4.8)

We observe that with $f$ defined in (4.3) and

$E(t)=\begin{cases}\epsilon_{0}(t+L)\exp(-t^{2}/2),&\mbox{if $t\in[-L,R]$,}\\\
0&\mbox{otherwise},\end{cases}$ (4.9)

we have

$F_{\mbox{shifted}}(t)\exp(-t^{2}/2)=\begin{cases}f(t)+E(t),&\mbox{ if $-L\leq
t\leq R$,}\\\ 0&\mbox{ otherwise}.\end{cases}$ (4.10)

To summarize, starting from the signal (1.1), we define $a_{k}$’s by (4.8),
$F$ by (4.1), $f$ by (4.3), and modify the noise $\epsilon_{0}$ in (1.1) as
per (4.9). In the rest of the paper, we will treat $f$ as the signal, and $E$
as the noise.

##  5 Discretization formulas

Another issue that must be addressed is the fact that in the context of signal
analysis, we know the values of $f$ only at finitely many points on $[-L,R]$;
these are the points at which the signal is sampled. Therefore, the expansion
of $f$ must be calculated in terms of Hermite functions using only these
values. This motivates our use of the Gauss quadrature formula to compute the
coefficients $\int_{\mathbb{R}}f(t)\psi_{k}(t)dt$. However, it is a delicate
matter to compute the quadrature nodes and weights [62, 63]. Moreover, this
requires the knowledge of $f$ at exactly the Gauss quadrature nodes, which is
not the case in practice. The existence of quadrature formulas based on
arbitrary nodes are known to exist [44], but it will be seen that for our
purposes we do not need quadrature formulas exact for integrating products of
Hermite functions of some order. Rather, it is sufficient to approximate $f$
using least squares approximation, provided the points at which $f$ is
evaluated are sufficiently dense in $[-L,R]$. This is encapsulated in the
following definition.

###### Definition 5.1

A sequence $\\{\nu_{m}\\}_{m=1}^{\infty}$ of possibly signed measures is
called a strongly Marcinkiewicz Zygmund sequence if each of the following two
conditions is satisfied.

$\int_{\mathbb{R}}|P(t)|^{2}d|\nu_{m}|(t)\sim\int_{\mathbb{R}}|P(t)|^{2}dt,\qquad
P\in\Pi_{m}.$ (5.1) $\|\nu_{m}\|_{TV}\sim m^{1/2}.$ (5.2)

By an abuse of notation, we may say that $\nu_{m}$ is a strongly Marcinkiewicz
Zygmund measure of order $m$, and write $\nu_{m}\in SMZ(m)$.

###### Example 5.1

Obviously, the Lebesgue measure on ${\mathbb{R}}$ satisfies (5.1) for every
$m\geq 1$, but not the condition (5.2). However, in view of Proposition 3.1,
the Lebesgue measure restricted to $[-2\sqrt{m},2\sqrt{m}]$ belongs to
$SMZ(m)$. $\blacksquare$

###### Example 5.2

Let $x_{m,m}<\cdots<x_{1,m}$ be the zeros of $\psi_{m}$, and $\lambda_{k,m}$
be the corresponding weights in the Gauss quadrature formula:

$\sum_{k=1}^{m}\lambda_{k,m}\exp(x_{k,m}^{2}/2)P(x_{k,m})=\int_{\mathbb{R}}P(t)\exp(-t^{2}/2)dt,\qquad
P\in\Pi_{2m}.$ (5.3)

Writing $w_{k,m}=\lambda_{k,m}\exp(x_{k,m}^{2})$, it follows that

$\sum_{k=1}^{m}w_{k,m}|P(x_{k,m})|^{2}=\int|P(t)|^{2}dt,\qquad P\in\Pi_{m}.$

Let $\nu_{m}$ be the measure that associates the mass $w_{k,m}$ with each
$x_{k,m}$. Clearly, $\nu_{m}$ satisfies (5.1). In view of [43, Theorem 8.2.7]
(used with $p=2$, $b=0$), $\nu_{m}$ satisfies (5.2) as well, so that
$\nu_{m}\in SMZ(m)$. $\blacksquare$

###### Example 5.3

There exist $c,c_{1}$ with the following property [2]: Let
$\\{x_{k,m}\\}_{k=1}^{m}$ be any real numbers, $x_{m,m}<\cdots<x_{1,m}$, such
that

$[x_{m,m},x_{1,m}]\subseteq[-\sqrt{2m}(1+cm^{-2/3}),\sqrt{2m}(1+cm^{-2/3})],$
(5.4)

and

$|x_{k,m}-x_{k+1,m}|\leq c_{1}m^{-1/2}.$ (5.5)

Let $\nu_{m}$ be the measure that associates the mass $x_{k,m}-x_{k+1,m}$ with
$x_{k,m}$. Then for $1\leq p<\infty$,

$\left|\|P\|_{\nu_{m};p}^{p}-\|P\|_{p}^{p}\right|\leq(1/4)\|P\|_{p}^{p},\qquad
P\in\Pi_{m}.$ (5.6)

Clearly, $\nu_{m}\in SMZ(m)$. $\blacksquare$

##  6 Least squares approximation

Given integers $m\geq n\geq 1$, a measure $\nu\in SMZ(m)$, and $f\in
C_{0}({\mathbb{R}})$, we can define the least square approximation to $f$ from
$\Pi_{n}$ by

$S_{n}(\nu;f)=\mathop{\hbox{{arg
min}}}_{P\in\Pi_{n}}\|f-P\|_{\nu;2}^{2}=\mathop{\hbox{{arg
min}}}_{P\in\Pi_{n}}\int|f(x)-P(x)|^{2}d|\nu|(x).$ (6.1)

_In the sequel, we will assume that $\nu$ is a positive measure, so that
$\nu=|\nu|$._

###### Remark 6.1

We note that $S_{n}(\nu;f)$ is defined entirely in terms of the restriction of
$f$ to $\mathsf{supp}(\nu)$. Thus, with the notation as in Section 3, and $f$
as defined in (4.3), if $\mathsf{supp}(\nu)\subset[-L,R]$, then
$S_{n}(\nu;f)=S_{n}(\nu;f|_{[-L,R]})$. $\blacksquare$

We will need an explicit expression for $S_{n}(\nu;f)$ in the form $\int
f(y)K_{n}(\nu;x,y)d\nu(y)$ for a reproducing kernel $K_{n}(\nu;x,y)$ (cf.
(6.8)).

To obtain an expression for this kernel, we define the Gram matrix
$\mathbf{G}$ by

$\mathbf{G}_{\ell,k}=\int\psi_{\ell}(x)\psi_{k}(x)d\nu(x),\qquad\ell,k=0,\cdots,m-1.$
(6.2)

A quadratic form for $\mathbf{G}$ is

$\sum_{\ell,k}d_{\ell}\overline{d_{k}}\mathbf{G}_{\ell,k}=\int\left|\sum_{j}d_{j}\psi_{j}\right|^{2}d\nu.$
(6.3)

Since

$\int\left|\sum_{j}d_{j}\psi_{j}(x)\right|^{2}dx=\sum_{j}|d_{j}|^{2},$

the estimate (5.1) and the Raleigh-Riesz principle shows that $\mathbf{G}$ is
positive definite, and all of its singular values are bounded from above and
below by constants independent of $m$. In particular, $\mathbf{G}$ is
invertible and positive definite, with all of the singular values of
$\mathbf{G}^{-1}$ also bounded from above and below by constants independent
of $m$. Thus, there exists a lower triangular matrix $\mathbf{L}$ so that the
Cholesky decomposition

$\mathbf{G}^{-1}=\mathbf{L}\mathbf{L}^{*}$ (6.4)

holds. Moreover, for all $\\{d_{\ell}\\}\subset{\mathbb{C}}$,

$\sum_{\ell,k}d_{\ell}\overline{d_{k}}\mathbf{G}_{\ell,k}\sim\sum_{\ell}|d_{\ell}|^{2}.$
(6.5)

Let

$\tilde{p}_{j}=\sum_{\ell=0}^{j}\mathbf{L}_{\ell,j}\psi_{\ell},\qquad
j=0,\cdots,m-1.$ (6.6)

Using the definition of the Gram matrix and the Cholesky decomposition in
(6.4), it is easy to verify that

$\int\tilde{p}_{j}(x)\tilde{p}_{k}(x)d\nu(x)=\delta_{j,k},\qquad
j,k=0,\cdots,m-1.$ (6.7)

Next, we define the reproducing kernel

$K_{n}(\nu;x,y)=\sum_{j=0}^{n-1}\tilde{p}_{j}(x)\tilde{p}_{j}(y)=\sum_{\ell,\ell^{\prime}=0}^{n-1}\mathbf{G}^{-1}_{\ell,\ell^{\prime}}\psi_{\ell}(x)\psi_{\ell^{\prime}}(y),\qquad
x,y\in{\mathbb{R}}.$ (6.8)

The Fourier transform of $K_{n}(\nu;\cdot,y)$ is given by

$\widehat{K}_{n}(\nu;x,y)=\mathfrak{F}(K_{n}(\nu;\cdot,y))(x)=\sum_{j=0}^{n-1}\mathfrak{F}(\tilde{p}_{j})(x)\tilde{p}_{j}(y)$
(6.9)

It is not difficult to verify that for $x\in{\mathbb{R}}$,

$S_{n}(\nu;f)(x)=\int
f(y)K_{n}(\nu;x,y)d\nu(y),\qquad\mathfrak{F}(S_{n}(\nu;f))(\omega)=\int
f(y)\widehat{K}_{n}(\nu;\omega,y)d\nu(y).$ (6.10)

Computationally, the procedure to find $S_{n}(\nu;f)$ is as follows. We define

$\hat{f}(\nu;\ell)=\int f(y)\psi_{\ell}(y)d\nu(y),\qquad\ell=0,\cdots,n-1,$
(6.11)

and solve the linear system of equations

$\sum_{k=0}^{n-1}\mathbf{G}_{\ell,k}d_{k}(f)=\hat{f}(\nu;\ell),\qquad\ell=0,\cdots,n-1.$
(6.12)

for $d_{k}(f)$’s. Then $S_{n}(\nu;f)$ is given by
$S_{n}(\nu;f)=\sum_{k=0}^{n-1}d_{k}(f)\psi_{k}$.

###### Remark 6.2

In the case when $\nu$ is the Lebesgue measure on ${\mathbb{R}}$,
$S_{n}(\nu;f)$ is the partial sum of the Hermite expansion. The corresponding
reproducing kernel is given by

$K_{n}(x,y)=\sum_{k=0}^{n-1}\psi_{k}(x)\psi_{k}(y),\qquad x,y\in{\mathbb{R}}.$
(6.13)

$\blacksquare$

##  7 Main theorems

As remarked earlier, if $\nu\in SMZ(m)$ is supported on $[-L,R]$, then

$S_{n}(\nu;F_{\mbox{shifted}}\exp(-\cdot^{2}/2))=S_{n}(\nu;f)+S_{n}(\nu;E),$
(7.1)

can be calculated using observed values of $F_{0}$.

Our first theorem is a general theorem, giving estimates on the error in
approximating $\mathfrak{F}(f)$ by the Fourier transform of the least square
approximation to $f$; i.e., $\mathfrak{F}(S_{n}(\nu;f))$. The first estimate
(7.2) has a desired weaker dependence on $n$, but the estimate is in
comparison with the best polynomial approximation of $f$ from $\Pi_{n/2}$
rather than from $\Pi_{n}$. This is useful when this degree of approximation
decreases polynomially with $n$. The second estimate (7.3) compares the error
in terms of the best approximation from $\Pi_{n}$ itself (the same degree as
$S_{n}(\nu;f))$. This exhibits a more rapid rate of convergence as a function
of $n$ when this degree of approximation decreases faster than any polynomial
in $1/n$, as is the case for (4.3) (cf. (9.27)), our main concern in this
paper.

###### Theorem 7.1

Let $f\in C_{0}\cap L^{1}$, $m\geq n\geq 1$ be integers and, $\nu\in SMZ(m)$
be a positive measure. Then

$\|\mathfrak{F}(f)-\mathfrak{F}(S_{n}(\nu;f))\|_{\infty}\lesssim\|f-S_{n}(\nu;f)\|_{1}\lesssim
E_{1,n/2}(f)+nE_{\infty,n/2}(f).$ (7.2)

Alternately,

$\|\mathfrak{F}(f)-\mathfrak{F}(S_{n}(\nu;f))\|_{\infty}\lesssim\|f-S_{n}(\nu;f)\|_{1}\lesssim
n^{1/2}E_{1,n}(f)+n^{3/2}E_{\infty,n}(f).$ (7.3)

When $f$ is a discrete Laplace transform as in (4.3), we need to compare the
desired Fourier transform in (4.6) with the Fourier transform of
$S_{n}(\nu;f)$. Here, $f$ being an entire function of exponential type, the
estimate (7.3) exhibits a more rapid rate of convergence with respect to $n$.

Although the following theorem is just a corollary of Theorem 7.1, we list it
as a separate theorem in view of its relation with the central problem of this
paper.

###### Theorem 7.2

Let $f$ be as defined in (4.3), $m\geq n\geq 1$ be integers, $\nu\in SMZ(m)$,
and $\mathsf{supp}(\nu)\subseteq[-L,R]$. Let $\Lambda=\max_{1\leq\ell\leq
K}\lambda_{\ell}$. Then

$\|\mathfrak{F}(f)-\mathfrak{F}(S_{n}(\nu;f))\|_{\infty}\lesssim\|f-S_{n}(\nu;f)\|_{1}\lesssim\frac{n^{3/4}}{\sqrt{n!}}(\Lambda/\sqrt{2})^{n}\exp(\Lambda^{2}/4)\sum_{\ell=1}^{K}|a_{\ell}|.$
(7.4)

Next, we turn our attention to the noise part, $S_{n}(\nu;E)$.

###### Theorem 7.3

Let $m\geq n\geq 1$ be integers, $\nu\in SMZ(m)$,
$\mathsf{supp}(\nu)\subseteq[-L,R]$. With the constant $C$ as in Lemma 9.1(b),
let

$C\sqrt{m}\geq 2n\geq c\log m$ (7.5)

for a suitable choice of $c$. We assume that the expected value of each $E(y)$
is $0$, and

$|E(y)|\leq S,\qquad y\in{\mathbb{R}}.$ (7.6)

Then for any $\delta\in(0,1)$, we have with probability $>1-\delta$

$\left\|\int
E(y)\hat{K}_{n}(\nu;w,y)d\nu(y)\right\|_{\infty}\lesssim\left(\frac{\log
m}{m}\right)^{1/4}\sqrt{\log\left(\frac{\log m}{\delta}\right)}.$ (7.7)

###### Remark 7.1

We can relax the assumption (7.6) to assume that $E$ is sub-Gaussian. However,
this just means that (7.6) holds with a high probability, and complicates the
notation without adding further insight. $\blacksquare$

The following corollary is a simple consequence of Theorem 7.3 and the
Höffding’s inequality Lemma 9.6.

###### Corollary 7.1

With the conditions as in Theorem 7.3, if we repeat the calculations with $M$
instances of $E(y)$: $E_{1}(y),\cdots,E_{M}(y)$, then with probability
$>1-\delta$

$\left\|\frac{1}{M}\sum_{j=1}^{M}\int
E_{j}(y)\hat{K}_{n}(w,y)d\nu(y)\right\|_{\infty}\lesssim
M^{-1/2}\left(\frac{\log m}{m}\right)^{1/4}\sqrt{\log\left(\frac{M\log
m}{\delta}\right)}.$ (7.8)

##  8 Numerical Results

Based on the above analysis, we define a procedure for determining the
$\lambda_{j}$’s in (1.1) which we term “Laplace to Fourier (L2F).” Most of our
examples will illustrate results for a biexponential signal with parameters
$A_{1}$=0.5, $A_{2}$=0.5, $T_{21}$=10ms, and $T_{22}$=50ms, with corresponding
$\lambda_{1}=1/T_{21}$ and $\lambda_{2}=1/T_{22}$; these parameters may of
course be varied to explore behavior in different regions of parameter space.

The first step in the L2F procedure is to expand the modified biexponential
signal in the time domain, (4.3), in Hermite polynomials. These polynomials,
up to degree 32, were produced using the Chebfun numerical computing package
[25]. The coefficients in the expansion were found through integration of
(6.11), using the Gauss quadrature expression defined in (6.12). Quadrature
nodes and weights were also produced using Chebfun packages. Having
established an approximation for (4.3), we proceeded to find an approximation
to its Fourier transform, (4.6).

As indicated in (6.10) the approximation of this FT follows from (3.3),
expressing the property of the Hermite functions as eigenfunction of the FT,
with eigenvalues $(i)^{degree}$. We now filter this signal using a custom-
designed finite-support filter (A.3) [26, 45, 46] to reduce sidelobes in the
frequency domain, with the goal of identifying dominant peaks at the positions
of the desired $\lambda$ values. However, we found that for realistic values
of $\lambda$ such as those that arise in MRI physics, the large range of
$\omega$ values in (4.6) required for resolution precluded this demonstration
(Figs. 1, 2). As an alternative, the signal shift described above enabled us
to identify the $\lambda$ corresponding to the component with the greater of
the two $\lambda$’s; the other component was substantially more attenuated by
the multiplicative exponential introduced in (4.3). Our approach then was to
evaluate successive shifts to identify this dominant peak without the need for
numerical optimization (e.g. nonlinear least squares analysis), thereby
reducing the optimization problem from its original form with four unknowns,
two of which enter nonlinearly, to three unknowns, one of which enters
nonlinearly.

Following estimation of the larger of the two time constants, which by
convention we denote by $T_{22}$, we now complete our L2F method by estimating
the remaining three parameters using NLLS. The reduction in dimensionality
from four to three through the above procedure is expected to result in
improved estimation of the remaining parameters. Thus, $T_{22}$ is inserted as
a fixed parameter into the NLLS objective function and the remaining
parameters, $A_{1}$,$A_{2}$, and $T_{21}$, are estimated using any one of the
numerous available methods. We used the MATLAB lsqnonlin function, based on
the trust region reflective algorithm with initial guesses for the $a_{k}$
taken as random values between 0 and 1 based on their interpretation as
component fractions that sum to one. The initial guess for $T_{21}$ was set to
a random value between 0ms and 300ms based on known properties of myelin-
constrained water in MRI. Feasible set constraints were set between 0 and 1
for coefficients and 1 and 300ms for $T_{21}$. For the native four-parameter
NLLS estimation problem, the same constraints and initial guesses used for
$T_{21}$ were implemented for $T_{22}$; again, for L2F, $T_{22}$ was obtained
from peak-finding without any other least-squares procedure.

The approximation of (4.3) using an expansion with Hermite polynomials is
shown in Figure 1. As seen, the Gauss-Hermite quadrature approach performs
extremely well for approximation of the time-domain signal. This provides the
basis for using of these coefficients to approximate the Fourier transform of
equation (4.3).

Figure 1: Approximation of a modified biexponential (4.3) decay curve with
parameters: $A_{1}=0.5$, $A_{2}=0.5$, $T_{21}=10$ms, $T_{22}=50$ms. The error
for this approximation, shown in the right panel, remains small throughout the
signal duration.

The Fourier transform of the modified signal is presented in (4.6). The
approximation of this signal, (6.10), is nearly indistinguishable from the
original in as seen in the difference plot in Fig. 2, right-hand panel. We
find that the approximation maintains fidelity for frequencies in the domain
of roughly [-6,6], although the displayed range in the difference plot (right-
hand panel) is [-3,3] in order to preserve detail in the region about zero.
Similar results were obtained for the imaginary part of the FT.

The restricted range over which the approximation to the FT remains accurate
presents severe constraints on our ability to resolve closely spaced signals.
As an alternative, we make use of the signal shift as described in (4.7) to
attenuate the amplitude of the peak corresponding to the smaller value of
$\lambda$, which is $\lambda_{2}$ according to our convention in which
$T_{22}>T_{21}$, so that the value of $\lambda_{1}$ can be readily identified
as the dominant remaining peak in the FT. A visualization of this procedure,
which is detailed in Section A, is shown in Figure 3.

Figure 2: The real part of the Fourier transform of the modified biexponential
signal (4.6)and its approximation through expansion with Hermite polynomials
(6.10). The same signal parameters detailed in Figure 1 are used. The fidelity
of this approximation to the FT of the signal supports our attempt to extract
a $\lambda$ value. However, the error in our approximation increases rapidly
outside of a relatively small frequency window (right-hand panel). Figure 3:
Illustration of $T_{22}$ estimation.The same signal parameters detailed in
Figure 1 are used. The vertical red line shows the theoretical value of
$\lambda_{1}$ in terms of $x$ (see Section A) while the red circle shows the
location of the point of maximum value (the estimated value of $\lambda$). For
these parameters, a $T_{22}$ value of approximately 49.03 was found, close to
the correct underlying value of 50.

Further numerical results were generated for time decay constants which are
more closely spaced, and hence more problematic to extract (see Table 1 and
Table 2). To model experimental data, we performed our analysis with different
levels of Gaussian noise added to the original biexponential signal. Thus,
results in these tables were obtained by adding noise to a biexponential
signal, with SNR defined as the signal peak (initial) value divided by the
standard deviation of the added Gaussian noise. Parameters were then extracted
from this signal using the NLLS and L2F procedures described above. This
process was performed for noiseless signals as well as over 500 noise
realizations for SNR values of $10^{6}$, $10^{5}$, and $10^{4}$. The mean,
standard deviation (StDev), and root mean squared error (RMSE) with respect to
the correct underlying value for each estimated parameter across noise
realizations are shown in the Tables.

Table 1 shows results for decay constants of $T_{21}=40$ms and $T_{22}=60$ms,
with $A_{1}=A_{2}=0.5$. For the noiseless signal, there is no variability in
the prediction of $T_{22}$ using the L2F method. Variability of L2F’s
remaining parameters-as well as all of the NLLS predictions-is attributable to
the randomness of initial guesses and the presence of local minima in the
optimization procedure used to obtain those estimates. In the noiseless case
and the nearly-noiseless case of SNR=$10^{6}$, our results are superior to
those obtained from NLLS. However, with increasing noise, error in the L2F
method increases rapidly due to the extent of shifting required to extract a
single decay constant; recall that, due to of the rapidly decaying exponential
in (4.3), this shift attenuates the amplitude of the peak corresponding to the
larger of the two $\lambda$’s (i.e. the smaller of the two $T_{2}$’s),
permitting the identification of the smaller $\lambda$ (i.e. the larger
$T_{2}$). Given that shifting also attenuates the peak related to the smaller
$\lambda$ (though not to the same extent as the larger $\lambda$), we see that
introducing too much of a shift causes noise in the signal to become dominant.

Table 1: Parameter estimation for varying levels of noise, where
$T_{21}=40$ms, $T_{22}=60$ms, and $A_{1}$=$A_{2}$=0.5. The mean, standard
deviation (StDev), and root mean squared error (RMSE) for each parameter are
calculated from 500 independent noise realizations. Results are shown for four
values of SNR. In the leftmost column, we see that for a noisless signal the
L2F method has smaller StDev and RMSE results compared to NLLS for all four
parameters-indicating the superiority of this method for this ideal case.
Performance of L2F is seen to decrease even for the very high SNR value of
SNR=$10^{6}$.For SNR $\leq$ $10^{5}$, L2F is seen to perform worse than NLLS
for all parameters, as shown by the corresponding values of StDev and RMSE.
Table 2: Parameter estimation for varying levels of noise, where
$T_{21}=50$ms, $T_{22}=60$ms, and $A_{1}$=$A_{2}$=0.5. Abbreviations are the
same as those in Table 1. In the case of a noisless signal, for all
parameters, L2F outperforms NLLS in terms of StDev and RMSE values. As noise
increases, L2F performance worsens, due to the increased StDev of the initial
$T_{22}$ estimate that is the basis for the method.

It is notable that NLLS performs essentially equivalently across all of these
levels of noise, all of which represent signals with very high SNR. However,
while the noiseless behavior of L2F, and its behavior for SNR = $10^{6}$, were
superior to NLLS, L2F exhibited more rapid degradation of overall performance
as SNR decreased. We attribute this decline in performance to the increasing
variability in the estimation of $T_{22}$. We see that in the noiseless case,
there is no variability in the $T_{22}$ prediction. However, upon the
introduction of even a small amount of noise, L2F’s estimation of $T_{22}$
exhibits considerable variability, as seen from the given StDev values for the
estimates. Recall that the remaining parameters are estimated using the fixed
value of $T_{22}$ recovered from L2F, so that variability in the latter leads
variability in these remaining parameters:$A_{1}$, $A_{2}$, and $T_{21}$.
Viewing 1, as SNR decreases, standard deviation and RMSE values are fairly
constant for the NLLS method when comparing the results of parameter
estimation at high and low SNRs. Performing the same comparison with the L2F
method, we observe that while L2F outperforms NLLS for the highest values of
SNR, its peformance degrades and is not as robust as NLLS for noisier signals.

Table 2 presents the same analysis for a case in which the decay constants of
the two signal components are yet more closely spaced:$T_{21}=50$ms and
$T_{22}=60$ms. Again, the L2F approach outperforms NLLS for noiseless signals
and for near-infinite SNR due to its stable prediction of the $T_{22}$
parameter. This indicated the theoretical efficacy of L2F even in the regime
where NLLS is unable to effectively distinguish between closely-spaced decay
constants. As noise increases, the increasing StDev of L2F’s $T_{22}$
prediction reduces the performance of this method. $T_{21}$ appears to suffer
the most and the coefficients $A_{1}$ and $A_{2}$ are also affected.
Continuing our analysis of the SNR=$10^{4}$ case, we see that while the StDevs
of these coefficients are smaller than those predicted by NLLS, bias in the
values of these predictions is the driver of their poor performance with the
L2F method.

Finally, we recognize that the multiexponential model we are working with in
this paper represents a potential oversimplification of the signal model in an
actual MRI experiment. In particular, the quantification of the estimated
parameters, including in the related method of T2 relaxation spectrum analysis
mentioned above in the Introduction, depends on several complexities in the
signal generation process and has been explored in previous work[1, 37, 38,
42, 65] . However, while such experiments, and in particular applications to
neuroscience, form the basis for our motivation in this study, we are here
concerned with the classic, and widely-applicable, multiexponential signal as
written in (1.1) and as interpreted throughout our manuscript.

##  9 Proofs

In order to prove Theorem 7.1, we need some preparation. First, we recall some
facts from the theory of weighted polynomial approximation. Parts (a) and (b)
are given (in a greater generality) in [43, Theorem 3.2.5]. Part (c) follows
from Proposition 3.1 and (9.1). Part (d) is easy to deduce using [43, Thoerem
3.4.2(b), Lemma 3.4.4(b)].

###### Lemma 9.1

Let $n\geq 1$ be an integer.
(a) We have

$K_{n}(x,x)=\sum_{j=0}^{n-1}\psi_{j}(x)^{2}\lesssim n^{1/2},\qquad
x\in{\mathbb{R}}.$ (9.1)

(b) There exists $C>0$ such that for $|x|\leq Cn^{1/2}(1+c/n)$,

$K_{n}(x,x)\gtrsim\sqrt{n}.$ (9.2)

(c) If $|y|>2n^{1/2}$, then

$|\hat{K}_{n}(w,y)|\lesssim\exp(-cn),\qquad w\in{\mathbb{R}}.$ (9.3)

(d) For any $P\in\Pi_{n}$, $1\leq p\leq\infty$, we have the Bernstein
inequality:

$\|P^{\prime}\|_{p}\lesssim\sqrt{n}\|P\|_{p}.$ (9.4)

###### Lemma 9.2

Let $m\geq 1$ be an integer, and $\nu\in SZM(m)$. With $C$ as in Lemma 9.1(b),
we have for $|x|\leq C\sqrt{m}$, and for some constant $C_{1}>0$,

$\nu([x-C_{1}/\sqrt{m},x+C_{1}/\sqrt{m}])\lesssim m^{-1/2}.$ (9.5)

Proof. Let $|x|\leq C\sqrt{m}$. Using the Bernstein inequality (9.4), the
Schwarz inequality, (9.1), and (9.2), we see that

$\displaystyle|K_{m}(x,y)-K_{m}(x,x)|$
$\displaystyle\leq|x-y|\|K_{m}^{\prime}(x,\circ)\|_{\infty}\lesssim\sqrt{m}|x-y|\|K_{m}(x,\circ)\|_{\infty}$
$\displaystyle\lesssim\sqrt{m}|x-y|\left\\{\sum_{k=0}^{m-1}\psi_{k}(x)^{2}\right\\}^{1/2}\left\\{\max_{y\in{\mathbb{R}}}\sum_{k=0}^{m-1}\psi_{k}(y)^{2}\right\\}^{1/2}$
$\displaystyle\lesssim m|x-y|K_{m}(x,x)^{1/2}\lesssim\sqrt{m}K_{m}(x,x)|x-y|.$

Consequently, for a judiciously chosen $C_{1}>0$, we have for $|x-y|\leq
C_{1}/\sqrt{m}$:

$\sqrt{m}\lesssim(1/2)K_{m}(x,x)\leq
K_{m}(x,y)\leq(3/2)K_{m}(x,x)\lesssim\sqrt{m}.$ (9.6)

Hence, using (5.1) and (9.1) again,

$\displaystyle m\nu([x-C_{1}/\sqrt{m},x+C_{1}/\sqrt{m}])$
$\displaystyle\leq\int_{x-C_{1}/\sqrt{m}}^{x+C_{1}/\sqrt{m}}K_{m}(x,y)^{2}d\nu(y)\leq\int_{\mathbb{R}}K_{m}(x,y)^{2}d\nu(y)$
(9.7)
$\displaystyle\lesssim\int_{\mathbb{R}}K_{m}(x,y)^{2}dy=K_{m}(x,x)\lesssim\sqrt{m}.$

This proves (9.5). $\blacksquare$

The following lemma establishes bounds on the integral norms of the kernels
defined in (6.8) and (6.9).

###### Lemma 9.3

Let $m\geq n\geq 1$ be integers, such that $C^{2}m\geq 4n$, where $C$ as in
Lemma 9.1(b), and $\nu\in SMZ(n)$. Then for every $x\in{\mathbb{R}}$,

$\sum_{j=0}^{n-1}|\tilde{p}_{j}(x)|^{2}\sim\sum_{j=0}^{n-1}\psi_{j}(x)^{2}\lesssim
n^{1/2},\quad$ (9.8)
$\int|\hat{K}_{n}(\nu;x,y)|^{2}d\nu(y)=\sum_{j=0}^{n-1}|\mathfrak{F}(\tilde{p}_{j})(x)|^{2}\sim\sum_{j=0}^{n-1}\psi_{j}(x)^{2}\lesssim
n^{1/2}.$ (9.9)

Moreover,

$\int|K_{n}(\nu;x,y)|d\nu(y)\lesssim n^{1/2},$ (9.10)
$\int|K_{n}(x,y)|d\nu(y)\lesssim n^{1/2}.$ (9.11)

Since the Lebesgue measure satisfies (5.1) for every $n$, we have

$\max\left(\int|K_{n}(x,y)|dy,\int|K_{n}(\nu;x,y)|dy\right)\lesssim
n^{1/2},\qquad n=1,2,\cdots.$ (9.12)

Moreover

$\max\left(\int|\hat{K}_{n}(\nu;x,y)|dy,\int|\hat{K}_{n}(\nu;x,y)|d\nu(y)\right)\lesssim
n^{1/2}$ (9.13)

Proof. In view of (6.8),

$\sum_{j=0}^{n-1}|\tilde{p}_{j}(x)|^{2}=\sum_{\ell,\ell^{\prime}=0}^{n-1}\mathbf{G}^{-1}_{\ell,\ell^{\prime}}\psi_{\ell}(x)\psi_{\ell^{\prime}}(x),$
(9.14)

which is a quadratic form for $\mathbf{G}^{-1}$. So, the estimate (9.8)
follows by using (6.5) with $\psi_{\ell}(x)$ in place of $d_{\ell}$. The
second estimate in (9.8) is the same as (9.1). The estimates in (9.9) are
obtained similarly, noting that

$\hat{K}(\nu;x,y)=\sum_{\ell,\ell^{\prime}=0}^{n-1}\mathbf{G}^{-1}_{\ell,\ell^{\prime}}(-i)^{\ell}\psi_{\ell}(x)\psi_{\ell^{\prime}}(x)$

To prove (9.10), we note that since $2\sqrt{n}\leq C\sqrt{m}$, Lemma 9.2 shows
that

$|\nu|([-2\sqrt{n},2\sqrt{n}])\lesssim\sqrt{nm}\sqrt{1/m}\lesssim\sqrt{n}.$

Hence, using Proposition 3.1, with $K_{n}(\nu;x,y)$ in place of $P$ and $p=1$,
followed by Schwarz inequality and (9.8),we deduce that

$\displaystyle\int_{\mathbb{R}}|K_{n}(\nu;x,y)|d|\nu|(y)$
$\displaystyle\lesssim\int_{-2\sqrt{n}}^{2\sqrt{n}}|K_{n}(\nu;x,y)|d|\nu|(y)$
$\displaystyle\lesssim\left\\{\int_{-2\sqrt{n}}^{2\sqrt{n}}|K_{n}(\nu;x,y)|^{2}d|\nu|(y)\right\\}^{1/2}\left\\{|\nu|([-2\sqrt{n},2\sqrt{n}])\right\\}^{1/2}\lesssim
n^{1/2}.$

The estimates (9.11) (and hence, (9.12)) follow similary. The estimate (9.13)
follows similarly as well, using (9.9)) and (6.7). $\blacksquare$

We note that we are working with a least squares approximation, but require
estimates in the uniform norm, while the Fourier transform of a function is
bounded uniformly by the $L^{1}$ norm of the function. To facilitate the
transitions among these norms, we recall the following lemma [43, Theorem
4.2.4] relating the various norms for weighted polynomials.

###### Lemma 9.4

For any integer $n\geq 1$ and $P\in\Pi_{n}$, we have

$\|P\|_{1}\lesssim n^{1/2}\|P\|_{\infty}\lesssim n\|P\|_{1}.$ (9.15)

Next, we obtain some constructions which approximate a function $f$ in
different norms at the same time; i.e., constructions which do not depend upon
the norm in which the degree of approximation is measured. It is not difficult
to construct a near-best approximation to $f$ from $\Pi_{n}$. Let
$n^{\prime}=\lfloor n/2\rfloor$. For $f\in L^{1}+L^{\infty}$, we write

$\hat{f}(k)=\int f(y)\psi_{k}(y)dy,$
$s_{\ell}(f)=\sum_{k=0}^{\ell-1}\hat{f}(k)\psi_{k},$

and

$V_{n}(f)=\frac{1}{n-n^{\prime}}\sum_{\ell=n^{\prime}}^{n-1}s_{\ell}(f).$
(9.16)

Part (a) of the following lemma is proved in [43, Lemma 4.1.5]. Part (b)
follows easily from (9.12) and the Riesz-Thorin interpolation theorem.

###### Lemma 9.5

If $1\leq p\leq\infty$, $f\in L^{p}$, then
(a)

$E_{p,n}\leq\|f-V_{n}(f)\|_{p}\lesssim E_{p,n/2}(f).$ (9.17)

and
(b)

$E_{p,n}\leq\|f-s_{n}(f)\|_{p}\lesssim n^{1/2}E_{p,n}(f).$ (9.18)

With this preparation, we are now ready to prove Theorem 7.1.

Proof of Theorem 7.1. Lemma 9.4 followed by Lemma 9.3 leads to

$\|S_{n}(\nu;f)\|_{1}\lesssim n^{1/2}\|S_{n}(\nu;f)\|_{\infty}\lesssim
n\|f\|_{\infty}.$ (9.19)

Next, we observe that for any $P\in\Pi_{n}$,

$S_{n}(\nu;P)(x)=\int P(y)K_{n}(\nu;x,y)d\nu(y)=P(x).$ (9.20)

Hence, using (9.19), we have for _any_ $P\in\Pi_{n}$,

$\|f-S_{n}(\nu;f)\|_{1}=\|f-P-
S_{n}(\nu;f-P)\|_{1}\leq\|f-P\|_{1}+\|S_{n}(\nu;f-P)\|_{1}\lesssim\|f-P\|_{1}+n\|f-P\|_{\infty}.$
(9.21)

We use this estimate with $V_{n}(f)$ in place of $P$ and use Lemma 9.5(a) to
arrive at (7.2). If we use $s_{n}(f)$ in place of $P$ in (9.21) instead, then
Lemma 9.5(b) leads to (7.3). $\blacksquare$

Next, we turn to the proof of Theorem 7.2, which again requires some further
preparation.

The generating function for Hermite polynomials is given by (cf. [61, Formula
(5.5.7)]))

$\exp(2xw-w^{2})=\pi^{1/4}\sum_{k=0}^{\infty}\frac{2^{k/2}}{\sqrt{k!}}h_{k}(x)w^{k}.$
(9.22)

Using $w=-\lambda/2$, this yields

$\exp(-x\lambda)=\pi^{1/4}\exp(\lambda^{2}/4)\sum_{k=0}^{\infty}\frac{(-1)^{k}}{2^{k/2}\sqrt{k!}}h_{k}(x)\lambda^{k},$

so that

$g_{\lambda}(x)=\exp(-x\lambda-x^{2}/2)=\pi^{1/4}\exp(\lambda^{2}/4)\sum_{k=0}^{\infty}\frac{(-1)^{k}}{2^{k/2}\sqrt{k!}}\psi_{k}(x)\lambda^{k}.$
(9.23)

Next, we recall that

$\|\psi_{k}\|_{1}\lesssim k^{1/4},\qquad\|\psi_{k}\|_{\infty}\lesssim
k^{-1/4}.$ (9.24)

We observe further using the Schwarz inequality that for any
$x\in{\mathbb{R}}$,

$\sum_{k=0}^{\infty}\frac{x^{k}}{\sqrt{k!}}=\sum_{k=0}^{\infty}\frac{(\sqrt{2}x)^{k}}{\sqrt{k!}}2^{-k/2}\leq\left\\{\sum_{k=0}^{\infty}\frac{(\sqrt{2}x)^{2k}}{k!}\right\\}^{1/2}\left\\{\sum_{k=0}^{\infty}2^{-k}\right\\}^{1/2}=\exp(x^{2}).$
(9.25)

Finally, we note that

$\frac{(n+k-1)!}{(n-1)!}=k!\left(\genfrac{}{}{0.0pt}{}{n+k-1}{n-1}\right)\geq
k!.$

Using (9.24), (9.25), and the above estimate, we deduce that

$\displaystyle\exp(\lambda^{2}/4)E_{1,n}(g_{\lambda})$
$\displaystyle\lesssim\sum_{k=n}^{\infty}\frac{k^{1/4}}{2^{k/2}\sqrt{k!}}\lambda^{k}\leq\sum_{k=n}^{\infty}\frac{1}{k^{1/4}\sqrt{(k-1)!}}(\lambda/\sqrt{2})^{k}\leq
n^{-1/4}\sum_{k=n}^{\infty}\frac{1}{\sqrt{(k-1)!}}(\lambda/\sqrt{2})^{k}$
(9.26)
$\displaystyle\leq\frac{(\lambda/\sqrt{2})^{n}}{n^{1/4}\sqrt{(n-1)!}}\sum_{k=0}^{\infty}\frac{1}{\sqrt{k!}}(\lambda/\sqrt{2})^{k}\lesssim\frac{(\lambda/\sqrt{2})^{n}}{n^{1/4}\sqrt{(n-1)!}}\exp(\lambda^{2}/2).$

We estimate $E_{\infty,n}(g_{\lambda})$ similarly and conclude that

$E_{1,n}(g_{\lambda})\lesssim\frac{n^{1/4}}{\sqrt{n!}}(\lambda/\sqrt{2})^{n}\exp(\lambda^{2}/4),\qquad
E_{\infty,n}(g_{\lambda})\lesssim\frac{n^{-1/4}}{\sqrt{n!}}(\lambda/\sqrt{2})^{n}\exp(\lambda^{2}/4).$
(9.27)

With this preparation, Theorem 7.2 is now easy to prove.

Proof of Theorem 7.2. We apply Theorem 7.1 with $f$ defined in (4.3), taking
$\nu$ to be supported on $[-L,R]$ (i.e., $\mathsf{supp}(\nu)\subseteq[-L,R]$).
We note that in view of Remark 6.1, we may assume that $f$ is defined on the
entire real axis rather than just on $[-L,R]$. We note that

$f(t)=\sum_{k=1}^{K}a_{k}g_{\lambda_{k}}(t).$

Hence, Theorem 7.2 is a simple consequence of Theorem 7.1 and the estimates
(9.27). $\blacksquare$

The next objective of this section is to estimate $\int
E(y)\hat{K}_{n}(w,y)d\nu(y)$. We will use the following lemma, Höffding’s
inequality ([9, Theorem 2.8]), several times. Note that this inequality
quantifies the degree to which the sum of the indicated random variables can
exceed a given value; it is an example of a concentration inequality.

###### Lemma 9.6

(Höffding’s inequality) Let $X_{1},\cdots,X_{M}$ be i.i.d. random variables
with mean $0$, and each $X_{k}\in[a_{k},b_{k}]$. Then for any $t>0$,

$\mathsf{Prob}\left\\{\left|\sum_{k=1}^{M}X_{k}\right|\geq t\right\\}\leq
2\exp\left(-\frac{t^{2}}{\sum_{k=1}^{M}(b_{k}-a_{k})^{2}}\right).$ (9.28)

We also need a covering lemma (see, e.g., [22]).

###### Lemma 9.7

There exists a set ${\mathcal{C}}$ with at most $\sim n$ elements such that
for any $P\in\Pi_{n}$,

$\|P\|_{\infty}\sim\max_{y\in{\mathcal{C}}}|P(y)|.$ (9.29)

With these preliminaries, we can now prove Theorem 7.3.

Proof of Theorem 7.3.

We observe first that in view of (5.2) and (3.8),

$\left\|\int_{|y|\geq
C\sqrt{m}}E(y)\hat{K}_{n}(\circ,y)d\nu(y)\right\|_{\infty}\lesssim\left\|\int_{|y|\geq
C\sqrt{m}}|E(y)||\hat{K}_{n}(\circ,y)|d\nu(y)\right\|_{\infty}\lesssim
Sm^{1/2}\exp(-cn).$ (9.30)

Now, let $w\in{\mathbb{R}}$ be fixed for the time being (until (9.36) is
proved). With $C_{1}$ as in Lemma 9.2, we divide the interval
$[-C\sqrt{m},C\sqrt{m}]$ into subintervals $I_{k}$, each of length at most
$2C_{1}/\sqrt{m}$, so that $|\nu|(I_{k})\lesssim m^{-1/2}$. We denote the
characteristic (indicator) function of $I_{k}$ by $\chi_{k}$. We then consider
the random variables

$X_{k}=\int_{I_{k}}E(y)\hat{K}_{n}(w,y)d\nu(y)=\int
E(y)\chi_{k}(y)\hat{K}_{n}(w,y)d\nu(y),$ (9.31)

so that

$\int_{|y|\leq Cm^{1/2}}E(y)\hat{K}_{n}(w,y)d\nu(y)=\sum_{k}X_{k}.$ (9.32)

Since the expected value of each $E(y)$ is $0$, that of each $X_{k}$ is $0$ as
well. Further, for each $k$,

$|X_{k}|=\left|\int E(y)\chi_{k}(y)\hat{K}_{n}(w,y)d\nu(y)\right|\leq
S\int|\chi_{k}(y)\hat{K}_{n}(w,y)|d\nu(y).$ (9.33)

Using the Schwarz inequality and Lemma 9.2, we see that

$\left(\int|\chi_{k}(y)\hat{K}_{n}(w,y)|d\nu(y)\right)^{2}\leq\nu(I_{k})\left\\{\int_{I_{k}}|\hat{K}_{n}(w,y)|^{2}d\nu(y)\right\\}\lesssim
m^{-1/2}\left\\{\int_{I_{k}}|\hat{K}_{n}(w,y)|^{2}d\nu(y)\right\\};$ (9.34)

i.e., with

$r_{k}=Sm^{-1/4}\left\\{\int_{I_{k}}|\hat{K}_{n}(w,y)|^{2}d\nu(y)\right\\}^{1/2},$
$|X_{k}|\lesssim r_{k}\qquad\mbox{for all }k.$

Using (9.9), we conclude that

$\sum_{k}r_{k}^{2}\lesssim S^{2}(n/m)^{1/2}.$ (9.35)

So, Höffding’s inequality implies that

$\mathsf{Prob}\left\\{\left|\int_{|y|\leq
Cm^{1/2}}E(y)\hat{K}_{n}(w,y)d\nu(y)\right|\geq
t\right\\}=\mathsf{Prob}\left\\{\left|\sum_{k}X_{k}\right|\geq t\right\\}\leq
2\exp\left(-\frac{t^{2}}{S^{2}(n/m)^{1/2}}\right).$ (9.36)

We note that this estimate holds for each $w\in{\mathbb{R}}$, and
$\sum_{k}X_{k}\in\Pi_{n}$. In view of the covering lemma, Lemma 9.7, this
leads to

$\mathsf{Prob}\left\\{\left\|\int_{|y|\leq
Cm^{1/2}}E(y)\hat{K}_{n}(\circ,y)d\nu(y)\right\|_{\infty}\geq
t\right\\}\lesssim n\exp\left(-\frac{ct^{2}}{S^{2}(n/m)^{1/2}}\right).$ (9.37)

Choosing

$t=cS(n/m)^{1/4}\sqrt{\log(n/\delta)},$

we deduce that with probability $>1-\delta$,

$\left\|\int_{|y|\leq
Cm^{1/2}}E(y)\hat{K}_{n}(\circ,y)d\nu(y)\right\|_{\infty}\lesssim
S(n/m)^{1/4}\sqrt{\log(n/\delta)}.$

Together with (9.30), we obtain (7.7) with a judicious choice of $n$.
$\blacksquare$

##  10 Conclusions

The analysis of multiexponential decaying signals is a longstanding problem in
mathematics and physics. We have introduced a novel method for parameter
estimation in this setting based on the fact that the Hermite polynomials are
eigenfunctions of the Fourier transform. This enables us to avoid direct use
of the notoriously ill-posed inverse Laplace transform for multiexponential
analysis. We present numerical results for the special case of biexponential
decay, and demonstrate the success of the method for a wide range of decay
constants, including in regimes where the more traditional method of NLLS
fails. We also note the sensitivity of the method to noise, as a sobering
reminder of its genesis in the inverse Laplace transform.

##  Appendix A Fourier-based estimation of the $\lambda$’s

In this section, we summarize how to obtain the values of $\lambda_{k}$ from
the Fourier transform $\mathfrak{F}(\tilde{\mu})(\omega)$ as in (4.6).

We consider a measure of the form

$\mu^{\prime}=\sum_{k=1}^{K}b_{k}\delta_{\lambda_{k}}.$

When we sample $\mathfrak{F}(\mu^{\prime})(\omega)$ at $\omega=\delta\ell$,
$|\ell|<N$ for some $\delta$ and write $\omega_{k}=\lambda_{k}\delta$, we
obtain the Fourier coefficients of a periodic measure $\mu$:

$\hat{\mu}(j)=\sum_{k=1}^{K}b_{k}\exp(-ij\omega_{k}),\qquad
j=-N+1,\cdots,N-1.$ (A.1)

Of course, it is assumed that the points $\omega_{k}\in[0,2\pi]$. In this
section $\|\cdot\|$ will denote $|(\cdot)\mbox{ mod }2\pi|$. The material in
this section is based on [26, 45, 46].

Clearly, the periodic measure $\mu$ can be recovered formally by its Fourier
series, and one is tempted to use the partial sum to take into account that
only finitely many Fourier coefficients are known. However, a judiciously
filtered Fourier sum avoids sidelobes and yields greater accuracy.
Accordingly, we introduce such a sum (see (A.3)) using a kernel. Let
$h:{\mathbb{R}}\to[0,1]$ be an infinitely differentiable, even function, such
that $h(t)=1/2$ for $|t|\leq 1/2$ and $h(t)=0$ for $|t|\geq 1$. We introduce
the kernel

$\Phi_{N}(t)=\sum_{\ell\in{\mathbb{Z}}}h\left(\frac{|\ell|}{N}\right)\exp(i\ell
t)=\sum_{|\ell|<N}h\left(\frac{|\ell|}{N}\right)\exp(i\ell t),\qquad
t\in{\mathbb{R}},\ N>0.$ (A.2)

It is not difficult to see that

$\sigma_{N}(\mu)(x)=\sum_{|\ell|<N}h\left(\frac{|\ell|}{N}\right)\hat{\mu}(\ell)\exp(i\ell
x)=\sum_{k=1}^{K}a_{k}\Phi_{N}(x-\omega_{k}).$ (A.3)

Now, the kernel $\Phi_{N}$ has the localization property that

$|\Phi_{N}(t)|\lesssim\frac{\Phi_{N}(0)}{\max(1,(N\|t\|)^{S})},\qquad N\geq
1,\ t\in{\mathbb{R}}.$ (A.4)

If $\eta=\min_{j\not=k}\|(\omega_{j}-\omega_{k}\|$, then it is clear that if
$N\eta\geq 4$ and $\|x-\omega_{k}\|\geq\eta/4$ for all $k$, then

$|\sigma_{N}(\mu)(x)|\lesssim\frac{\Phi_{N}(0)}{(N\eta)^{S}}\sum_{k=1}^{K}|a_{k}|.$
(A.5)

On the other hand, for any $x$ which does not satisfy the above condition,
there is a unique $j$ such that $\|x-\omega_{j}\|<\eta/4$. We have

$|\sigma_{N}(\omega_{j})|\geq\Phi_{N}(0)|a_{j}|-c\frac{\Phi_{N}(0)}{(N\eta)^{S}}\sum_{k=1}^{K}|a_{k}|$

for a suitable constant $c>0$. Using the Bernstein inequality, a concentration
inequality related to Hoeffding’s inequality that incorporates the variance of
the summed random variables, for the trigonometric polynomials, it follows
that for sufficiently large values of $N$, the prominent peaks of
$|\sigma_{N}(x)|$ appear within $c_{1}/N$ of the values of the $\omega_{j}$’s.
We refer the reader to the above cited literature for further details.

## References

  * [1] A. Akhondi-Asl, O. Afacan, R. V. Mulkern, and S. K. Warfield.  T2-relaxometry for myelin water fraction extraction using Wald distribution and extended phase graph.  In Medical Image Computing and Computer-Assisted Intervention–MICCAI 2014: 17th International Conference, Boston, MA, USA, September 14-18, 2014, Proceedings, Part III 17, pages 145–152. Springer.
  * [2] D. Batenkov, L. Demanet, and H. N. Mhaskar.  Stable soft extrapolation of entire functions.  Inverse Problems, 35(1):015011, 2019.
  * [3] M. Bertero, P. Boccacci, and E. R. Pike.  On the recovery and resolution of exponential relaxation rates from experimental data: a singular-value analysis of the Laplace transform inversion in the presence of noise.  Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences, 383(1784):15–29, 1982.
  * [4] M. Bertero, P. Boccacci, and E. R. Pike.  On the recovery and resolution of exponential relaxation rates from experimental data ii. The optimum choice of experimental sampling points for Laplace transform inversion.  Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences, 393(1804):51–65, 1984.
  * [5] M. Bertero, P. Brianzi, P. Parker, and E. R. Pike.  Resolution in diffraction-limited imaging, a singular value analysis.  Optica Acta: International Journal of Optics, 31(2):181–201, 1984\.
  * [6] M. Bertero, P. Brianzi, and E. R. Pike.  On the recovery and resolution of exponential relaxational rates from experimental data: Laplace transform inversions in weighted spaces.  Inverse Problems, 1(1):1, 1985.
  * [7] C. Bi, K. Fishbein, M. Bouhrara, and R. G. Spencer.  Stabilization of parameter estimates from multiexponential decay through extension into higher dimensions.  Scientific Reports, 12(1), 2022.
  * [8] J. Bonny, A. Traore, M. Bouhrara, R. G. Spencer, and G. Pages.  Parsimonious discretization for characterizing multi‐exponential decay in magnetic resonance.  NMR in Biomedicine, 33(12), 2020.
  * [9] S. Boucheron, G. Lugosi, and P. Massart.  Concentration inequalities: A nonasymptotic theory of independence.  Oxford university press, 2013.
  * [10] M. Bouhrara, J. S. Alisch, N. Khattar, R. W. Kim, A. C. Rejimon, L. E. Cortina, W. Qian, L. Ferrucci, S. M. Resnick, and R. G. Spencer.  Association of cerebral blood flow with myelin content in cognitively unimpaired adults.  BMJ Neurology Open, 2(1), 2020.
  * [11] M. Bouhrara, L. E. Cortina, N. Khattar, A. C. Rejimon, S. Ajamu, D. S. Cezayirli, and R. G. Spencer.  Maturation and degeneration of the human brainstem across the adult lifespan.  Aging, 13(11):14862, 2021.
  * [12] M. Bouhrara, L. E. Cortina, A. C. Rejimon, N. Khattar, C. Bergeron, J. Bergeron, D. Melvin, L. Zukley, and R. G. Spencer.  Quantitative age-dependent differences in human brainstem myelination assessed using high-resolution magnetic resonance mapping.  Neuroimage, 206:116307, 2020.
  * [13] M. Bouhrara, R. W. Kim, N. Khattar, W. Qian, C. M. Bergeron, D. Melvin, L. M. Zukley, L. Ferrucci, S. M. Resnick, and R. G. Spencer.  Age‐related estimates of aggregate g‐ratio of white matter structures assessed using quantitative magnetic resonance neuroimaging.  Human Brain Mapping, 42(8):2362–2373, 2021.
  * [14] M. Bouhrara, D. A. Reiter, H. Celik, J.-M. Bonny, V. Lukas, K. W. Fishbein, and R. G. Spencer.  Incorporation of Rician noise in the analysis of biexponential transverse relaxation in cartilage using a multiple gradient echo sequence at 3 and 7 tesla.  Magnetic resonance in medicine, 73(1):352–366, 2015.
  * [15] M. Bouhrara, D. A. Reiter, and R. G. Spencer.  Bayesian analysis of transverse signal decay with application to human brain.  Magnetic resonance in medicine, 74(3):785–802, 2015.
  * [16] M. Bouhrara, A. C. Rejimon, L. E. Cortina, N. Khattar, C. M. Bergeron, L. Ferrucci, S. M. Resnick, and R. G. Spencer.  Adult brain aging investigated using BMC-mcDESPOT–based myelin water fraction imaging.  Neurobiology of Aging, 85:131–139, 2020.
  * [17] K. R. Brownstein and C. E. Tarr.  Importance of classical diffusion in nmr studies of water in biological cells.  Physical Review A, 19(6):2446–2453, 1979.
  * [18] J. Brübach, J. Janicka, and A. Dreizler.  An algorithm for the characterisation of multi-exponential decay curves.  Optics and Lasers in Engineering, 47(1):75–79, 2009.
  * [19] H. Celik, M. Bouhrara, D. A. Reiter, K. W. Fishbein, and R. G. Spencer.  Stabilization of the inverse Laplace transform of multiexponential decay through introduction of a second dimension.  Journal of Magnetic Resonance, 236:134–139, 2013.
  * [20] E. M. Charles-Edwards and M. d. Nandita.  Diffusion-weighted magnetic resonance imaging and its application to cancer.  Cancer imaging, 6(1):135, 2006.
  * [21] C. K. Chui and H. N. Mhaskar.  Signal decomposition and analysis via extraction of frequencies.  Applied and Computational Harmonic Analysis, 40(1):97–136, 2016\.
  * [22] A. Cloninger and H. N. Mhaskar.  Cautious active clustering.  Applied and Computational Harmonic Analysis, 54:44–74, 2021.
  * [23] D. C. Dean, S. A. Hurley, S. R. Kecskemeti, J. P. O’Grady, C. Canda, N. J. Davenport-Sis, C. M. Carlsson, H. Zetterberg, K. Blennow, S. Asthana, M. A. Sager, S. C. Johnson, A. L. Alexander, and B. B. Bendlin.  Association of amyloid pathology with myelin alteration in preclinical Alzheimer disease.  JAMA Neurology, 74(1):41, 2017.
  * [24] G. S. Drenthen, W. H. Backes, A. P. Aldenkamp, G. J. Op’t Veld, and J. F. Jansen.  A new analysis approach for T2 relaxometry myelin water quantification: Orthogonal matching pursuit.  Magnetic Resonance in Medicine, 81(5):3292–3303, 2019.
  * [25] T. A. Driscoll, N. Hale, and L. N. Trefethen.  Chebfun Guide.  Pafnuty Publications, 2014.
  * [26] F. Filbir, H. N. Mhaskar, and J. Prestin.  On the problem of parameter estimation in exponential sums.  Constructive Approximation, 35(3):323–343, 2012.
  * [27] D. G. Gardner, J. C. Gardner, G. Laush, and W. W. Meinke.  Method for the analysis of multicomponent exponential decay curves.  The Journal of Chemical Physics, 31(4):978–986, 1959.
  * [28] Z. Gong, M. Bilgel, M. Kiely, C. Triebswetter, L. Ferrucci, S. M. Resnick, R. G. Spencer, and M. Bouhrara.  Lower myelin content is associated with more rapid cognitive decline among cognitively unimpaired individuals.  Alzheimer’s & Dementia, 19(7):3098–3107, 2023.
  * [29] S. J. Graham, P. L. Stanchev, and M. J. Bronskill.  Criteria for analysis of multicomponent tissue relaxation data.  Magnetic Resonance in Medicine, 35(3):370–378, 1996.
  * [30] G. Hampton, R. Neff, T. Hecht, P.-W. Fox, and R. Spencer.  Conditon number phase transitions for model selection in magnetic resonance relaxometry.  American Physical Society: April Meeting, 2023.
  * [31] R. Harman and F. Fairman.  Exponential approximation via closed-form Gauss-Newton method.  IEEE Transactions on Circuit Theory, 20(4):361–369, 1973.
  * [32] L. Héliot and A. Leray.  Simple phasor-based deep neural network for fluorescence lifetime imaging microscopy.  Scientific Reports, 11(1):23858, 2021.
  * [33] K. Holmström and J. Petersson.  A review of the parameter estimation problem of fitting positive exponential sums to empirical data.  Applied Mathematics and Computation, 126(1):31–61, 2002.
  * [34] A. A. Istratov and O. F. Vyvenko.  Exponential analysis in physical phenomena.  Review of Scientific Instruments, 70(2):1233–1257, 1999.
  * [35] O. Jahn, S. Tenzer, and H. B. Werner.  Myelin proteomics: Molecular anatomy of an insulating sheath.  Molecular Neurobiology, 40(1):55–72, 2009.
  * [36] E. Kavroulakis, P. G. Simos, G. Kalaitzakis, T. G. Maris, D. Karageorgou, I. Zaganas, S. Panagiotakis, M. Basta, A. Vgontzas, and E. Papadaki.  Myelin content changes in probable Alzheimer’s disease and mild cognitive impairment: Associations with age and severity of neuropsychiatric impairment.  Journal of Magnetic Resonance Imaging, 47(5):1359–1372, 2018.
  * [37] E. Kirtil and M. H. Oztop.  Food Engineering Reviews, 8(1):1–22, 2016.
  * [38] D. Kumar, S. Siemonsen, C. Heesen, J. Fiehler, and J. Sedlacik.  Noise robust spatially regularized myelin water fraction mapping with the intrinsic B1‐error correction based on the linearized version of the extended phase graph model.  Journal of Magnetic Resonance Imaging, 43(4):800–817, 2016.
  * [39] C. Laule, I. M. Vavasour, G. R. W. Moore, J. Oger, D. K. B. Li, D. W. Paty, and A. L. Mackay.  Water content and myelin water fraction in multiple sclerosis.  Journal of Neurology, 251(3):284–293, 2004.
  * [40] A. MacKay, C. Laule, I. Vavasour, T. Bjarnason, S. Kolind, and B. Mädler.  Insights into brain microstructure from the T2 distribution.  Magnetic resonance imaging, 24(4):515–525, 2006.
  * [41] A. Mackay, K. Whittall, J. Adler, D. Li, D. Paty, and D. Graeb.  In vivo visualization of myelin water in brain by magnetic resonance.  Magnetic Resonance in Medicine, 31(6):673–677, 1994.
  * [42] N. Mehdizadeh and A. H. Wilman.  Myelin water fraction mapping from multiple echo spin echoes and an independent B1+map.  Magnetic Resonance in Medicine, 88(3):1380–1390, 2022.
  * [43] H. N. Mhaskar.  Introduction to the theory of weighted polynomial approximation.  Volume 56. World Scientific Singapore, 1996.
  * [44] H. N. Mhaskar.  Local approximation using Hermite functions.  In Progress in Approximation Theory and Applicable Complex Analysis, pages 341–362. Springer, 2017.
  * [45] H. N. Mhaskar and J. Prestin.  On the detection of singularities of a periodic function.  Advances in Computational Mathematics, 12(2-3):95–131, 2000.
  * [46] H. N. Mhaskar and J. Prestin.  On local smoothness classes of periodic functions.  Journal of Fourier Analysis and Applications, 11(3):353–373, 2005\.
  * [47] A. F. Oberhauser, P. E. Marszalek, H. P. Erickson, and J. M. Fernandez.  The molecular elasticity of the extracellular matrix protein tenascin.  Nature, 393(6681):181–5, 1998.  Oberhauser, A F Marszalek, P E Erickson, H P Fernandez, J M England 1998/05/29 Nature. 1998 May 14;393(6681):181-5. doi: 10.1038/30270.
  * [48] T. Parasram, R. Daoud, and D. Xiao.  T2 analysis using artificial neural networks.  Journal of Magnetic Resonance, 325:106930, 2021.
  * [49] R. Piessens.  A new numerical method for the inversion of the Laplace transform.  IMA Journal of Applied Mathematics, 10(2):185–192, 1972.
  * [50] R. Piessens.  A bibliography on numerical inversion of the Laplace transform and applications.  Journal of Computational and Applied Mathematics, 1(2):115–128, 1975.
  * [51] G. F. Piredda, T. Hilbert, J. Thiran, and T. Kober.  Probing myelin content of the human brain with MRI: A review.  Magnetic Resonance in Medicine, 85(2):627–652, 2021.
  * [52] D. Potts and M. Tasche.  Parameter estimation for exponential sums by approximate prony method.  Signal Processing, 90(5):1631–1642, 2010.
  * [53] R. d. Prony.  Essai éxperimental et analytique: sur les lois de la dilatabilité de uides élastique et sur celles de la force expansive de la vapeur de l’alkool, à di érentes températures.  J. Ec. Polytech, 1:24, 1795.
  * [54] S. W. Provencher.  A Fourier method for the analysis of exponential decay curves.  Biophysical Journal, 16(1):27–41, 1976.
  * [55] M. Rozowski, J. Palumbo, J. Bisen, C. Bi, M. Bouhrara, W. Czaja, and R. G. Spencer.  Input layer regularization for magnetic resonance relaxometry biexponential parameter estimation.  Magnetic Resonance in Chemistry, 60(11):1076–1086, 2022.
  * [56] B. W. Rust, D. P. O’Leary, and K. M. Mullen.  Exponential Data Fitting and Its Applications, chapter 8, pages 145–164.  Bentham Science Publishers, 2010.
  * [57] A. Sharafi, D. Xia, G. Chang, and R. R. Regatte.  Biexponential T1$\rho$ relaxation mapping of human knee cartilage in vivo at 3 T.  NMR in Biomedicine, 30(10):e3760, 2017.
  * [58] H. Shinmoto, K. Oshio, A. Tanimoto, N. Higuchi, S. Okuda, S. Kuribayashi, and R. V. Mulkern.  Biexponential apparent diffusion coefficients in prostate cancer.  Magnetic resonance imaging, 27(3):355–359, 2009.
  * [59] D. Sonnenschmidt.  Measurement of water flow velocity by analysis of short temperature steps.  Meas. Sci. Technol., (7):1536, 1996.
  * [60] R. G. Spencer and C. Bi.  A tutorial introduction to inverse problems in magnetic resonance.  NMR in Biomedicine, 33(12), 2020.
  * [61] G. Szegö.  Orthogonal polynomials.  In Colloquium Publications/American Mathematical Society. Volume 23. Providence, 1975.
  * [62] A. Townsend, T. Trogdon, and S. Olver.  Fast computation of Gauss quadrature nodes and weights on the whole real line.  IMA Journal of Numerical Analysis, 36(1):337–358, 2015.
  * [63] L. N. Trefethen.  Exactness of quadrature formulas.  Siam Review, 64(1):132–150, 2022.
  * [64] I. M. Vavasour, C. M. Clark, D. K. B. Li, and A. L. MacKay.  Reproducibility and reliability of MR measurements in white matter: Clinical implications.  NeuroImage, 32(2):637–642, 2006.
  * [65] V. Wiggermann, I. Vavasour, S. Kolind, A. L. MacKay, G. Helms, and A. Rauscher.  Non‐negative least squares computation for in vivo myelin mapping using simulated multi‐echo spin‐echo T2 decay data.  NMR in Biomedicine, 33(12):e4277, 2020.
  * [66] D. Xiao, Y. Chen, and D. D.-U. Li.  One-dimensional deep learning architecture for fast fluorescence lifetime imaging.  IEEE Journal of Selected Topics in Quantum Electronics, 27(4):1–10, 2021.