Bayesian uncertainty quantification for magnetic resonance fingerprinting

Magnetic resonance imaging (MRI) (Nishimura 1996) is a non-invasive medical imaging technique. It is applied in clinical routine for diagnosing diseases and monitoring therapies. While standard qualitative MRI produces images that show good contrasts, it does not yield quantitative results for tissue-related parameters like the spin relaxation times. It has been shown that such quantitative information could be beneficial for detecting deseases (Tofts 2005). Nevertheless, many methods for quantitative MRI (qMRI) (Schmitt et al 2004, Warntjes et al 2007, 2008), require long acquisition times which can be difficult to achieve in clinical practice.

In the past few years, Magnetic Resonance Fingerprinting (MRF) (Ma et al 2013) has evolved into a promising and fast method for qMRI. The main idea is to acquire and analyze a sequence of images from data that are highly undersampled in the spatial Fourier domain. Application of a pseudo-inverse Fourier transformation such as Greengard and Lee (2004) yields a sequence of images with undersampling (aliasing) artifacts. From these images individual time series of the dynamic magnetization are extracted for each voxel and subsequently correlated with a dictionary. The dictionary contains time series for different tissue types calculated for the applied pulse sequence using a physical model. For each voxel tissue-related parameters (T₁ and T₂ relaxation times) are then determined as those parameters of the physical model for which the corresponding dictionary entry shows highest correlation with the time series obtained from the measured data. While MRF has originally been used for qMRI of the brain (Ma et al 2013), recent applications include the measurement of blood flow (Flassbeck et al 2019), cerebral blood volume (Christen et al 2014) and cardiac tissue characterization (Hamilton et al 2017).

Several improvements have been proposed for MRF recently. For example, to accelerate the dictionary matching method a singular value decomposition of the dictionary has been used in McGivney et al (2014). In Hoppe et al (2017), Virtue et al (2017), Balsiger et al (2018), Cohen et al (2018), Fang et al (2019) machine learning has been applied to overcome the discrete nature of the dictionary. Another route of development has been to adopt statistical approaches that fit the MRF data directly in the Fourier domain (Zhao et al 2016, Wübbeler and Elster 2017, Metzner et al 2019). In this way aliasing effects caused by the pseudo-inverse Fourier transform used in the dictionary matching method are avoided. The drawback, however, is that fitting the data in the Fourier domain leads to a non-convex, large-scale optimization task since the tissue parameters of all voxels need to be adjusted at the same time. This optimization task is challenging and it is very time-consuming to solve. Furthermore, good starting values have to be assigned in order to reach the correct solution.

The dictionary matching method usually restricts itself to the calculation of estimates for the tissue parameters. As part of the matching process the correlation between the measured course of magnetization at a voxel and all entries of the chosen dictionary is calculated. The maximum correlation provides a measure of the quality of the fit. However, it does not yield an uncertainty of the tissue-related parameters. Uncertainties that reliably characterize the accuracy of the tissue-related parameter estimates (Polders et al 2012) are essential to assess the significance of observed differences in single qMRI results, which is particularly relevant in applications such as therapy monitoring. While the statistical approaches (Zhao et al 2016, Wübbeler and Elster 2017, Metzner et al 2019) that fit the data directly in the Fourier domain allow for such an uncertainty quantification, their practical applicability is challenged due to the involved large-scale optimization tasks.

The goal of this paper is to develop an uncertainty characterization of the estimates obtained by the dictionary matching method and to demonstrate its practical feasibility using MRF phantom measurements. The approach applies Bayesian statistics (Gelman et al 2013) and yields posterior probability distributions for the tissue-related parameters at all voxels. These probability distributions allow probability statements about the tissue-related parameters to be made conditional on the observed data, and to quantify the uncertainties associated with the estimates of the tissue-related parameters obtained by the dictionary matching method. The proposed Bayesian approach employs so-called noninformative priors (Kass and Wasserman 1996) combined with an informative prior for the relaxation times that utilizes as prior knowledge the ranges chosen for the relaxation times when calculating the dictionary. The calculation of results is facilitated through the analytical derivation of the marginal posterior for the relaxation times and can be carried out using numerical integration. Since these calculations can effectively be carried out using the pre-computed dictionary, they are also fast.

The practicability of the approach is validated by analyzing MRF measurements of a phantom that consists of several tubes. Relaxation times can be assumed constant within each tube, which allows for a consistency check by comparing the within-tube variability of their estimates obtained by the dictionary matching method with the proposed uncertainty quantification for these estimates. In addition, the estimates determined by the dictionary matching method are compared with estimates obtained by reference measurements using long acquisition times in terms of the proposed uncertainty quantification for the former. Validity of the proposed uncertainty evaluation for MRF is also assessed in terms of simulated data with known ground truth, and applicability to an in vivo measurement is demonstrated as well.

The proposed approach for uncertainty estimation in MRF fingerprinting is fast and can be applied in combination with any dictionary matching method. MATLAB^® source code will be shared to facilitate its application (Metzner and Wübbeler 2021).

The paper is organized as follows. Section 2 briefly describes the structure of MRF data and introduces the dictionary matching method. Then the proposed Bayesian uncertainty quantification is developed. In section 3 the experiment is described for recording MRF data. Results achieved by the proposed Bayesian uncertainty quantification for numerical simulations and phantom experiments are then presented in section 4. The results obtained for the phantom measurements are also compared with additional reference measurements of the relaxation times using long acquisition times. In addition, results for an in vivo measurement are presented. In section 5 then a discussion is provided and some conclusions are finally given in section 6.

This section briefly outlines the structure of MRF data and the dictionary matching method. Then the proposed Bayesian approach is developed. A description of the physical model used for compiling the dictionary as well as mathematical details of the Bayesian approach are given in the appendices.

2.1. MRF data

In MRI nuclear spins are excited in a selected slice by applying a high frequency pulse in the presence of a static magnetic field. The acquired data are samples in the Fourier domain of the image. The measured magnetization depends on the tissue parameters, namely relaxation times T₁, T₂ and proton spin density ρ, as well as on the MR sequence parameters (e.g. echo time, repetition time, RF-pulse design). In MRF a sequence of pulses with varying flip angles is applied, and measurements are conducted between subsequent pulses. These measurements do not cover the whole Fourier domain and can thus be carried out fast. The dynamics of the magnetization depends on the whole sequence of pulses and can be modeled through a three-dimensional linear discrete time system, see appendix A. Relevant parameters that are controlled in the design of the sequence of pulses are flip angles, repetition times and echo times.

The dependency of the dynamics of the magnetization on the flip angles and repetition times is not indicated in the following. To keep notations simple, we also suppress the presence of measurements using multiple receiver coils (as used in our experiments) and neglect the spatially varying received sensitivities of the multiple receiver coils in the following description. This simplification does not affect our treatment and complete relationships are given in appendix A.

Let N denote the number of considered voxels in the selected slice, and ${\theta }_{1}={({T}_{1}^{1},\,\ldots ,\,{T}_{1}^{N})}^{T}$, ${\theta }_{2}={({T}_{2}^{1},\,\ldots ,\,{T}_{2}^{N})}^{T}$, ${\theta }_{3}={({\rho }^{1},\,\ldots ,\,{\rho }^{N})}^{T}$ the sought vectors of relaxation times T₁, T₂ and proton spin density ρ, respectively, for N voxels. As is common in MRI, we model the proton spin density in terms of a complex number which can account for phase changes in the measured data (e.g. due to susceptibility artefacts). Let m_l = m_l (θ₁, θ₂, θ₃), l = 1, ..., L denote the N × 1 complex-valued vector of magnetizations at the N voxels after applying the lth excitation pulse, and let mⁱ denote the L × 1 complex-valued vector of magnetizations at the ith voxel after applying all L pulses, see figure 1. The sequence of magnetizations at the ith voxel, mⁱ , depends only on the unknown tissue parameters at that voxel (and the known parameters of the applied pulse sequence), i.e. ${m}^{i}={m}^{i}({T}_{1}^{i},{T}_{2}^{i},{\rho }^{i})$. Furthermore, the model m_l = m_l (θ₁, θ₂, θ₃) depends linearly on θ₃ (and hence mⁱ also depends linearly on ρⁱ ), which will be utilized in the proposed Bayesian treatment.

Zoom In Zoom Out Reset image size

The measured data consist of the complex-valued vectors

$$\begin{eqnarray}&&{z}_{l}={F}_{l}({m}_{l}),\,\,l=1,\,\ldots ,\,L.\end{eqnarray} \tag{ 2.1 }$$

F_l (m_l ) denotes the n_l × 1 vector of measurements in the Fourier domain after the lth pulse, where typically n_l ≪ N holds. Hence, the magnetization m_l cannot be reconstructed exactly. The mapping F_l is indexed by l as the locations of measurements in Fourier domain are generally different for the different pulses. Note that each entry of z_l depends on the tissue parameters of all voxels.

2.2. Dictionary matching MRF

The principles of the dictionary matching method according to Ma et al (2013) are briefly described in the following. In a first step, the magnetizations m_l , l = 1, ...L, are approximately reconstructed from (2.1) through

$$\begin{eqnarray}&&{\widehat{m}}_{l}={F}_{l}^{\dagger }({z}_{l}),l=1,\,\ldots ,\,L,\end{eqnarray} \tag{ 2.2 }$$

where ${F}_{l}^{\dagger }$ denotes a pseudo-inverse of the Fourier transform mapping F_l in (2.1). From the sequence of estimated magnetizations, ${\widehat{m}}_{1},\,\ldots ,\,{\widehat{m}}_{L}$, the estimated course of magnetization ${\widehat{m}}^{i}$ for each voxel i is then obtained and subsequently denoted by

$$\begin{eqnarray}&&{\widehat{m}}^{i}={\left({\left({\widehat{m}}_{1}\right)}_{i},\,\ldots ,\,{\left({\widehat{m}}_{L}\right)}_{i}\right)}^{T},\end{eqnarray} \tag{ 2.3 }$$

see also figure 1. These estimates are generally not accurate as they contain aliasing errors caused by the imperfect inversion ${F}_{l}^{\dagger }$. In our calculations, we employed the pseudo-inverse Fourier transformation from Greengard and Lee (2004). The basic idea of MRF is to choose flip angles and repetition times of the pulse sequence randomly in such a way that the sequence

$$\begin{eqnarray}&&{\epsilon }^{i}={\left({\left({\widehat{m}}^{i}-{m}^{i}\right)}_{1},{\left({\widehat{m}}^{i}-{m}^{i}\right)}_{2},\,\ldots ,\,{\left({\widehat{m}}^{i}-{m}^{i}\right)}_{L}\right)}^{T}\end{eqnarray} \tag{ 2.4 }$$

of (complex-valued) temporal errors in the reconstruction of the magnetizations (2.2) for each voxel i can be modeled as noise. The reconstruction errors ⁱ are usually dominated by aliasing effects caused by the imperfect inverse Fourier mapping, and these effects depend on the course of the magnetization itself and the undersampling scheme which both are determined through the choice of the pulse sequence parameters. Since the magnetization at each voxel depends on the tissue parameters of only that voxel, estimates of those tissue parameters can then be obtained by fitting the modeled dynamics of the magnetization of each voxel to the observed approximate dynamics (2.3) at that voxel. In Ma et al (2013) this fitting process is done by pre-computing a dictionary of (normalized) courses of magnetizations on a discrete set of possible values of the sought relaxation times and then selecting that dictionary entry which shows maximum correlation to the observed dynamics (2.2). The spin density, which enters linearly into the model for the magnetization, is determined subsequently such that a best fit to the data results.

Basically, this fitting can be achieved likewise without pre-computing a dictionary. We show in appendix B that the dictionary matching method yields estimates ${({({\widehat{\theta }}_{1})}_{i},{({\widehat{\theta }}_{2})}_{i},({\widehat{\theta }}_{3})}_{i})=({\widehat{T}}_{1}^{i},{\widehat{T}}_{2}^{i},{\widehat{\rho }}^{i})$ at the ith voxel which are given by minimizing

$$\begin{eqnarray}&&\parallel {\widehat{m}}^{i}-{m}^{i}({T}_{1}^{i},{T}_{2}^{i},{\rho }^{i}){\parallel }^{2}\end{eqnarray} \tag{ 2.5 }$$

with respect to ${T}_{1}^{i},{T}_{2}^{i},{\rho }^{i}$, provided that the resolution of the dictionary is made arbitrarily fine; in other words, the dictionary matching method basically fits the modeled dynamics of the magnetization independently at each voxel to the observed rough estimates (2.3) at that voxel. In (2.5) ∥ · ∥ stands for the usual 2-norm in ${{\mathbb{C}}}^{L}$.

2.3. Bayesian inference for dictionary matching method

A Bayesian inference is carried out to quantify the uncertainty associated with the tissue-related parameters calculated by the dictionary matching method. The Bayesian inference is based on a statistical model for the observed magnetizations in (2.3) and the choice of prior distributions. In order to improve readability, some additional notations are introduced first.

2.3.1. Notations

Let the 2L × 1 real-valued vector ${\hat{y}}^{i}$ denote real part and imaginary part of the observed estimates ${\widehat{m}}^{i}$ of the course of magnetizations at voxel i, i.e.

$$\begin{eqnarray}&&{\hat{y}}^{i}={\left({\mathfrak{R}}{\left({\widehat{m}}^{i}\right)}^{T},{\mathfrak{I}}{\left({\widehat{m}}^{i}\right)}^{T}\right)}^{T},\end{eqnarray} \tag{ 2.6 }$$

and the 2 × 1 real-valued vector μⁱ the corresponding proton spin density

$$\begin{eqnarray}&&{\mu }^{i}={\left({\mathfrak{R}}({\rho }^{i}),{\mathfrak{I}}({\rho }^{i})\right)}^{T},\end{eqnarray} \tag{ 2.7 }$$

where ${\mathfrak{R}}(\cdot )$ and ${\mathfrak{I}}(\cdot )$ stand for real and imaginary part, respectively. Furthermore, let the 2L × 2 real-valued matrix ${M}^{i}={M}^{i}({T}_{1}^{i},{T}_{2}^{i})$ be defined such that

$$\begin{eqnarray}&&{m}^{i}({T}_{1}^{i},{T}_{2}^{i},{\rho }^{i})={M}^{i}{\mu }^{i}\end{eqnarray} \tag{ 2.8 }$$

contains real and imaginary part of the modeled course of magnetizations ${m}^{i}({T}_{1}^{i},{T}_{2}^{i},{\rho }^{i})$ at voxel i, see appendix A for details.

In accordance with the dictionary matching method the subsequent statistical analysis is carried out independently at each voxel, and we will hence omit the reference to a particular voxel i in the following, e.g. we use M = M(T₁, T₂) to denote ${M}^{i}={M}^{i}({T}_{1}^{i},{T}_{2}^{i})$ for some considered voxel i.

Probability distributions will be loosely denoted through their arguments, e.g. π(T₁, T₂, μ) stands for a probability density function of T₁, T₂ and μ.

2.3.2. Statistical model

Following the rationale of the dictionary matching method, the errors in the observed course of magnetization in (2.3) are modeled as independently and identically distributed with zero mean and variance σ². In addition, we assume normality of this distribution. More precisely, we assume

$$\begin{eqnarray}&&y\sim N\left(M\mu ,{\sigma }^{2}{{I}}_{2L}\right),\end{eqnarray} \tag{ 2.9 }$$

where M = M(T₁, T₂) depends on T₁ and T₂ only, μ denotes real and imaginary part of the spin density ρ, and I_2L stands for the identity matrix of dimension 2L, i.e. the observed data $\hat{y}$ are modeled as a realization of a Gaussian distribution with mean vector M μ and covariance matrix σ² I_2L . Model (2.9) is a homoscedastic Gaussian model with a mean vector that depends linearly on the spin density μ and nonlinearly on the relaxation times T₁ and T₂. Note that the variance σ² in (2.9) is allowed to differ for different voxels.

2.3.3. Bayesian inference

A Bayesian inference combines the prior knowledge about the unknowns with the information contained in the data. The result is the posterior probability distribution which is obtained via Bayes theorem according to

$$\begin{eqnarray}&&\pi ({T}_{1},{T}_{2},\mu ,{\sigma }^{2}| y)\propto \pi ({T}_{1},{T}_{2},\mu ,{\sigma }^{2})l({T}_{1},{T}_{2},\mu ,{\sigma }^{2};y).\end{eqnarray} \tag{ 2.10 }$$

π(T₁, T₂, μ, σ²) denotes the prior beliefs about the unknowns T₁, T₂, μ, σ², and l(T₁, T₂, μ, σ²; y) is the likelihood function which is determined through the statistical model (2.9) for the data.

As a prior distribution we choose

$$\begin{eqnarray}&&\pi ({T}_{1},{T}_{2},\mu ,{\sigma }^{2})=\pi ({T}_{1},{T}_{2})\pi (\mu )\pi ({\sigma }^{2})\propto {{\bf{1}}}_{{\rm{\Omega }}}({T}_{1},{T}_{2})/{\sigma }^{2},\end{eqnarray} \tag{ 2.11 }$$

where Ω denotes the domain for T₁ and T₂ in the dictionary, and the function 1_Ω(T₁, T₂) equals 1 for (T₁, T₂) ∈ Ω and 0 otherwise. The prior π(T₁, T₂) ∝ 1_Ω(T₁, T₂) thus represents the prior belief about the range of possible values for T₁ and T₂, without a priori preferring any single value from this range. The spin density acts as a location parameter in the likelihood, and the prior π(μ) ∝ 1 is a standard noninformative prior in such a case, see Kass and Wasserman (1996). Similarly, π(σ²) ∝ 1/σ² is a standard noninformative prior for the scale parameter σ². Hence, we combine the available prior knowledge about the relaxation times with standard noninformative priors for the proton spin density and for the variance σ² in (2.9). Note that albeit the prior (2.11) is improper, i.e. cannot be normalized, propriety of the posterior (2.10) can be ensured, see appendix B.

From the joint posterior π(T₁, T₂, μ, σ²∣y) in (2.10) marginal posterior distributions such as, e.g. π(T₁∣y) are obtained by integrating out the remaining unknowns. Once the marginal posterior for one of the tissue parameters has been determined, the uncertainty u(T₁) associated with the estimate ${\widehat{T}}_{1}$ obtained by the dictionary matching method can be calculated as

$$\begin{eqnarray}&&{u}^{2}({\widehat{T}}_{1})=\int \pi ({T}_{1}| y){\left({T}_{1}-{\widehat{T}}_{1}\right)}^{2}{\rm{d}}{T}_{1}.\end{eqnarray} \tag{ 2.12 }$$

In addition, a 95% credible interval I can be determined which satisfies

$$\begin{eqnarray}&&0.95={\int }_{I}\pi ({T}_{1}| y){\rm{d}}{T}_{1}.\end{eqnarray} \tag{ 2.13 }$$

Relation (2.13) does not uniquely determine an interval I, and an additional condition needs to be posed, for example the interval of minimal length satisfying (2.13) or the choice of a probabilistically symmetric interval. One can also use the marginal posterior to make probability statements, for example to determine the probability that T₁ exceeds a certain limit, ${\overline{T}}_{1}$ say, according to

$$\begin{eqnarray}&&P({T}_{1}\gt {\overline{T}}_{1}| y)={\int }_{{\overline{T}}_{1}}^{\infty }\pi ({T}_{1}| y){\rm{d}}{T}_{1}.\end{eqnarray} \tag{ 2.14 }$$

In appendix B we derive the marginal joint posterior for T₁, T₂ as

$$\begin{eqnarray}&&\pi ({T}_{1},{T}_{2}| y)\propto \displaystyle \frac{1}{\sqrt{{\det }({M}^{T}M)}}\times \displaystyle \frac{{{\bf{1}}}_{{\rm{\Omega }}}({T}_{1},{T}_{2})}{{\left(\parallel y-M\mu ({T}_{1},{T}_{2}){\parallel }^{2}\right)}^{L-1}}\end{eqnarray} \tag{ 2.15 }$$

with

$$\begin{eqnarray}&&\mu ({T}_{1},{T}_{2})={\left[{M}^{T}M\right]}^{-1}{M}^{T}y\end{eqnarray} \tag{ 2.16 }$$

and M = M(T₁, T₂). Furthermore,

$$\begin{eqnarray}&&\mu | {T}_{1},{T}_{2},y\sim {t}_{2L-2}\left(\mu ({T}_{1},{T}_{2}),\displaystyle \frac{\parallel y-M\mu ({T}_{1},{T}_{2}){\parallel }^{2}}{2L-2}{\left[{M}^{T}M\right]}^{-1}\right),\end{eqnarray} \tag{ 2.17 }$$

i.e. the conditional posterior distribution for μ given T₁, T₂ is a scaled and shifted 2-dimensional t-distribution with 2L − 2 degrees of freedom. Since L ≫ 1, a very good approximation is given by the Gaussian distribution

$$\begin{eqnarray}&&\mu | {T}_{1},{T}_{2},y\sim N\left(\mu ({T}_{1},{T}_{2}),\displaystyle \frac{\parallel y-M\mu ({T}_{1},{T}_{2}){\parallel }^{2}}{2L-2}{\left[{M}^{T}M\right]}^{-1}\right),\end{eqnarray} \tag{ 2.18 }$$

where, again, M = M(T₁, T₂). The marginal posterior for μ is then obtained as

$$\begin{eqnarray}&&\pi (\mu | y)=\int \pi (\mu | {T}_{1},{T}_{2},y)\pi ({T}_{1},{T}_{2}| y){\rm{d}}{T}_{1}{\rm{d}}{T}_{2}.\end{eqnarray} \tag{ 2.19 }$$

The 2 × 2 covariance matrix $U(\widehat{\mu })$ associated with the estimate $\widehat{\mu }$ obtained by the dictionary matching method is simply given through

$$\begin{eqnarray}&&U(\widehat{\mu })=\int \left\{\left(\widehat{\mu }-\mu ({T}_{1},{T}_{2})\right){\left(\widehat{\mu }-\mu ({T}_{1},{T}_{2})\right)}^{T}+\displaystyle \frac{\parallel y-M\mu ({T}_{1},{T}_{2}){\parallel }^{2}}{2L-2}{\left[{M}^{T}M\right]}^{-1}\right\}\pi ({T}_{1},{T}_{2}| y){{dT}}_{1}{{dT}}_{2},\end{eqnarray} \tag{ 2.20 }$$

where again M = M(T₁, T₂).

Expressions (2.15), (2.18) and (2.20) allow for a practical calculation of results using 2-d approximate numerical integration, which is most efficiently carried out on the dictionary itself. Note that the columns of matrix M can be pre-calculated and also stored in the dictionary for the chosen set of values for T₁, T₂. Once the marginal posterior distributions for T₁, T₂, μ have been determined, uncertainties such as (2.12), or probability statements as in (2.14) can be determined again using numerical integration. Details are given in appendix C.

The shared MATLAB^® source code (Metzner and Wübbeler 2021) provides the proposed Bayesian uncertainty quantification using numerical integration. The required input consists of the employed dictionary, together with the observed sequence of magnetizations in (2.3) for each voxel as well as estimates for T₁ and T₂. Note that it is generally possible to use other estimates than the dictionary-based MRF estimates. Its output are the discretized marginal posteriors for T₁ and T₂ and the uncertainties for T₁ and T₂ calculated according to (2.12). Due to the use of the pre-computed dictionary, calculations are carried out fast.The computation time depends on the size of the dictionary and on the number of voxels of interest. The calculation time of the uncertainties shown in this paper was about 90 minutes, whereas calculating the estimates by the dictionary matching method took about 6 minutes. When the size of the dictionary is reduced by a factor of 4, also the computation time for the uncertainties is reduced by a factor of 4.

The phantom data was acquired on a 3T Verio MR Scanner (Siemens Healthineers, Erlangen, Germany) with a Fast Imaging with Steady State Precession (FISP) sequence starting with an inversion pulse using a radial readout, each with 384 data points per radial line. Repetition times (TR) were chosen to be 12.54 ms, the inversion time (TI) as 42 ms, echo times (TE) as 7 ms, and the flip angle follows the magnitude of a sinusoidal curve with maximum 70 degree as proposed in Flassbeck et al (2019). A total number of 1000 radial lines were acquired and for the parameter estimation one image per radial line was reconstructed using the pseudo-inverse Fourier transformation (Greengard and Lee 2004). Figure 2 shows two reconstructed magnetization images and the sum of all the reconstructions. Data were acquired with an 8-channel head coil, and the coil sensitivity maps required for image reconstruction were obtained from the data themselves. The inplane resolution of the MRF scan was 0.78 × 0.78 mm² with a slice thickness of 5.0 mm and a field of view (FoV) of 300 × 300 mm². The calculation of the dictionary was carried out using the Bloch equations (see appendix A) assuming 201 isochromats in each voxel.

Zoom In Zoom Out Reset image size

Reference values for T₁ were determined by a Cartesian inversion recovery spin echo sequence acquisition with 7 inversion times (TI = 25–2400 ms, TE/TR = 12/8000 ms). Reference values for T₂ were determined using a Cartesian spin echo sequence with multiple echo times (TE = 24–1000 ms, TR = 3000 ms). For both acquisitions, inplane resolution was 0.83 × 0.83 mm² with a slice thickness of 5.0 mm and a FoV of 159 ×144 mm². A 2-parameter nonlinear least squares fitting algorithm was used to determine the reference T₁ and T₂ values.

The simulation closely follows the phantom measurement. All scanning parameters (TR, TE, TI, flip angles) and k-space locations for the radial points are set to the values chosen for the phantom measurements. Again, a total number of 1000 radial lines were used, and for the dictionary matching method parameter estimation was carried out using one image per radial line. The same complex coil sensitivity maps were chosen as for the phantom scan. We assigned values for T₁, T₂ and the complex proton density that match those obtained in the analysis of the phantom data. Magnetization during the simulated data acquisition was modeled by the same FISP model used for the calculation of the dictionary (see appendix A). The modeled magnetization was then transformed into the Fourier domain using Fourier transformation (Greengard and Lee 2004), and Gaussian noise with similar signal-to-noise ratio as for the real scanner data was added with a standard deviation of 3.43 × 10⁻⁶ (the maximum value of the real and the imaginary part of the simulated magnetizations are of the order 10⁻³).

The in vivo data acquisition was performed on a 3T Verio MR Scanner (Siemens Healthineers, Erlangen, Germany) using a 32-channel cardiac coil (Siemens Healthineers, Erlangen, Germany) and the coil sensitivity maps required for image reconstruction were obtained from the data themselves. A cardiac scan was realized in a short-axis orientation utilizing a FISP sequence with a total acquisition time of 12 s (allowing for an acquisition during a single breathhold). A radial readout was applied, each with 320 data points per radial line. Repetition times (TR) were chosen to be 8.2 ms, the inversion time (TI) as 21 ms, echo times (TE) as 4 ms, and the flip angle follows the magnitude of a sinusoidal curve with maximum 70 degree. A total number of 1500 radial lines were acquired and for the parameter estimation one image per radial line was reconstructed using the pseudo-inverse Fourier transformation. The resolution of the MRF scan was 1.3 × 1.3 mm² with a FOV of 320 × 320 mm² and a slice thickness of 8 mm. To reduce the impact of the slice profile, a bandwidth time product of 8 was chosen (SINC pulse length = 1920 ms). ECG triggering was used to start the acquisition during diastole, but data were then acquired continuously without further triggering.

Real and simulated MRF measurements of the phantom (see figure 2) were conducted and analyzed by the proposed uncertainty quantification for the dictionary matching method. In addition, reference measurements for T₁ and T₂ were carried out to assess these results. The dictionary was calculated according to the physical model given in appendix A for T₁ = (200, 205, ..., 3000) ms and T₂ = (20, 22, ..., 300) ms. The phantom consists of 9 different tubes and the values for T₁ and T₂ within each tube can be assumed to be constant. The variation of T₁ and T₂ between the tubes almost covers their range in the chosen dictionary, see figure 4.

4.1. Phantom and in vivo measurements

Figure 3 illustrates the (estimated) residuals in (2.4) for a typical voxel. The residuals largely follow a Gaussian distribution, supporting the statistical model chosen for the Bayesian uncertainty quantification.

Zoom In Zoom Out Reset image size

In figure 4 calculated uncertainties are compared with the observed variation of estimated values for T₁ and T₂ within each tube. More precisely, the figure shows mean values of the estimates obtained by the dictionary matching method for T₁ and T₂ within the different tubes; in addition, standard deviations of the estimates within the tubes and within-tube means of uncertainties are displayed for the different tubes. Mean uncertainties reflect the within-tube variability of the estimates well, yet they are slightly larger than the observed variability of the estimates.

Zoom In Zoom Out Reset image size

In order to analyze the impact of the range of the values in the dictionary we also evaluated the phantom data with a larger and a smaller dictionary in T₂. When taking a larger dictionary with T₂ = (20, 22, ..., 500) ms, the within-tube standard deviations of the estimates differed from those obtained for the original dictionary by at most 10%. This is also true for the mean uncertainties, except for tube nine in which the mean uncertainty is 50% higher. Results for a smaller dictionary with T₂ = (20, 22, ..., 100) ms are shown in figure 5. In this case the T₂ estimates can no longer cover the whole range of the underlying T₂ values, and for tubes with a higher T₂ than allowed by the dictionary the corresponding estimates equal the maximum value of T₂ in the dictionary. Note that the resulting uncertainties indicate that the true T₂ values are beyond the upper bound induced by the (restricted) dictionary.

Zoom In Zoom Out Reset image size

In figure 6 the same quantities are plotted as in figure 4, but for data where 5 radial lines in k-space are acquired per image. That means, the amount of data is increased by a factor of 5. It can be seen that the mean uncertainties and the standard deviation of the estimates in each tube are smaller now (roughly a factor $\sqrt{5}$). The ratio of mean uncertainties and standard deviations of the estimates remains unchanged.

Zoom In Zoom Out Reset image size

Figure 7 compares the reference measurements with the results obtained by the dictionary matching method in terms of the uncertainties calculated for the latter. For T₁ estimation the calculated uncertainties appear to describe well the deviations between results of the dictionary matching method and the reference measurements. Regarding T₂ estimation this is true for the smaller values of T₂, while for larger values of T₂ a slight bias in the estimates produced by the dictionary matching method can be observed. Yet, the calculated uncertainties still characterize the size of the errors well in most cases.

Zoom In Zoom Out Reset image size

In figure 8 marginal distributions for T₁ and T₂ in two voxels of different tubes are shown. The dotted red lines illustrate the bounds for T₁ and T₂. Especially for T₂ some of the distributions are truncated at these bounds, i.e. the bounds are truly informative. The derived marginal distributions can be used, for example, to calculate the probability that T₁ is smaller than a certain limit. In figure 8 the probabilities for T₁ < 1600 ms of the two voxels indicated by the arrows are 0.54 and almost 1, respectively.

Zoom In Zoom Out Reset image size

Estimates for T₁ and T₂ as well as their uncertainties of a cardiac in vivo measurement are shown in figure 9. Note that in some areas the uncertainty increases towards the centre of the body, which is reasonable since the coil sensitivity is lower there.

Zoom In Zoom Out Reset image size

4.2. Simulated MRF

In figure 10 calculated uncertainties are compared with the observed variation of the estimates for T₁ and T₂ within each tube. The uncertainties reflect the different within tube variabilities of the estimates well. Similarly to the results for the real MRF measurements the calculated uncertainties appear to be slightly conservative.

Zoom In Zoom Out Reset image size

For higher T₁ and T₂ values the associated uncertainty is also higher. Nonetheless, the relative uncertainty is roughly constant. We simulated further MRF data by using the same simulation setup, changing only repetition times (TR was set to 11 ms) or flip angle (the maximum of the magnitude of the sinusoidal curve was set to 60 degree). In both cases, the relative uncertainties were again similar, the mean uncertainties in all nine tubes differed by a maximum of 10% when using the different repetition times, and by 20% when applying the different flip angle.

Figure 11 compares the estimates obtained by the dictionary matching method with the known ground truth in terms of the uncertainties calculated for the former. Here only the results for one tube are shown which are representative also for the other tubes. Calculated uncertainties describe well the size of the observed errors in the estimates obtained by the dictionary matching method for both T₁ and T₂. Again, the calculated uncertainties appear to be slightly conservative.

Zoom In Zoom Out Reset image size

The residuals in figure 3 obtained for the course of the reconstructed magnetizations from the phantom data in a single voxel appear to largely follow a Gaussian distribution. This finding supports the statistical model assumed for the Bayesian uncertainty quantification. Note that the results shown in figure 3 are typical and residuals in other voxels and other tubes look similarly. In figure 2 two single reconstructed magnetization images are shown with severe aliasing artifacts. When looking at the sum of all these individual reconstructions, the aliasing artifacts appear to be canceled out which supports the assumption that aliasing artifacts act like noise over time.

The simulated MRF data were constructed by adding in k-space Gaussian noise of realistic size. When carrying out the simulations without this additional measurement noise, results for the estimates of T₁ and T₂ essentially remain the same. This implies that the (unavoidable) aliasing errors in the course of the observed magnetizations at each voxel are the dominant uncertainty source. We have shown that the dictionary matching method is equivalent to fitting the modeled magnetizations to the observed ones, which can be recommended only under the assumption that the (aliasing) errors appear like random noise, ideally following a Gaussian distribution. The statistical model employed for the Bayesian uncertainty quantification is based on the very same assumption. Hence, when the design of the MRF data acquisition has been chosen properly to meet the assumptions underlying the the dictionary matching method, one can be confident that the assumptions of the proposed Bayesian uncertainty quantification are met as well.

The squared uncertainty has been taken as the expected quadratic loss of the estimate obtained by the dictionary matching method under the derived posterior. It describes the a posteriori belief about the size of the errors. We observed that the within-tube variances of the estimates are slightly smaller than these uncertainties for both the real MRF data (figure 4) and the simulated MRF data (figure 10). If the within-tube variances can be taken as a reference for the size of the errors, then this size is slightly overestimated by our Bayesian inference. One possible source of error is the presence of temporal and/or spatial correlations of the residuals. Although the assumption of a Gaussian distribution for the errors in the sequence of magnetization is in line with the assumption underlying the dictionary matching method, the statistical model will generally provide an approximation to the true distribution of errors only, and represents a limitation of the proposed approach. Future work could address more sophisticated statistical models. However, one benefit of the proposed method is that calculations can be carried out by numerical integration using only the dictionary and the observed data in an efficient way, which is hardly possible when more sophisticated statistical models, containing additional unknowns, are employed. Another possible source of error in the proposed approach is the calculation of the Bayesian uncertainty quantification itself. The calculation of integrals such as (2.12) are carried out on the dictionary and might be prone to numerical approximation errors. For the results shown in this paper this error source is not relevant since the results remain essentially unchanged when doubling the resolution of the dictionary. However, basically this is a potential source of error in the uncertainty calculation (as well as in the dictionary matching method itself), especially when dictionaries with low resolution are taken.

While for the simulated data the set of errors in the estimates produced by the dictionary matching method are well characterized by the calculated uncertainties (figure 11), this does not hold for all results obtained for the measured MRF data (figure 7). In the latter case a significant bias can be observed for the estimates obtained by the dictionary matching method for the larger T₂ values, albeit the size of single errors of the estimates are still largely captured by the calculated uncertainties. Since the reference measurements for T₂ can be trusted, and as their uncertainty is small compared to the observed discrepancies to the dictionary matching method, the results of the dictionary matching method appear to have some bias for large T₂ values. The reason likely is some error in the underlying dynamic model used to create the dictionary. The Bayesian uncertainty quantification is carried out assuming the employed physical model to be correct and therefore underrates calculated uncertainties in cases when a significant model error occurs. One reason for an error in the employed physical model are inhomogeneities of the B₁ field. The RF excitation in our physical model is assumed to be perfect (which would only be achieved with an infinite sinc pulse), but this is a simplification which can lead to errors in T₂.

The fact that the unavoidable (large) aliasing errors are the dominant uncertainty source for the dictionary matching method makes this method—as well as the proposed uncertainty quantification—robust. That is, as long as other error sources like measurement errors, incorrect k-space positioning, or errors in the dynamic model of magnetization, are small compared to the aliasing errors, both the dictionary matching method and the proposed uncertainty quantification can be expected to produce reasonable results.

We used an isochromat simulation with 201 isochromats per voxel to model the magnetization (see appendix A). It is also possible to apply the Extended Phase Graph formalism (Weigel 2015) for which no further parameters such as the number of isochromats are needed. When comparing results from both models, the mean values of the estimates of all nine tubes have a deviation less than 1%, and the same holds for within-tube standard deviations and for mean uncertainties.

The results of the estimates as well as of the uncertainties depend on the range of the dictionary. When taking a dictionary with a larger range of values for the relaxation times, the mean uncertainties (across single tubes) behave similarly compared to the (within-tube) standard deviations of the estimates. When taking a dictionary with a too small range for the relaxation times, biased estimates result and calculated uncertainties can no longer be expected to provide a reasonable characterization of the errors in the estimates. Hence, the range of the dictionary plays an important role both for the results achieved by the dictionary matching method and the proposed uncertainty quantification, and it should therefore be carefully chosen.

In Zhao et al (2018) the Cramér–Rao bound (CRB) was used to characterize the accuracy of estimates obtained in an MRF experiment, and to optimize its protocol (in (Zhao et al 2018) the design parameters are chosen as the sequence parameters repetition times and flip angle). The proposed Bayesian uncertainty quantification also depends on the chosen protocol, and it could likewise be used for protocol optimization. Future work could address this feature of our method and compare resulting protocols with those achieved by CRB minimization.

It was also possible to apply the proposed uncertainty calculation to in vivo data. Reasonable results for the estimates of the dictionary matching method as well as for their uncertainties could be achieved.

One important advantage of Bayesian inference is that a probability distribution is obtained. Basically, the whole distribution can serve as an uncertainty quantification. In some cases, a single summary derived from it such as the uncertainty as defined in (2.12) may be sufficient. In other cases, when the distributions deviate significantly from a Gaussian distribution (see the truncated posteriors in figure 8 for example), the whole distribution should be considered. The posterior distribution can be used to calculate probabilities about a hypothesis such as ${T}_{1}\gt {\underline{T}}_{1}$. In figure 12 a probability map of the hypothesis T₁ < 500 ms is shown for the phantom data. It can be seen that T₁ is almost certainly smaller than 500 ms in all voxels of the tube on the bottom right, while the opposite is true for, e.g. the tube on the top left. In other tubes a clear decision cannot be made, which is represented through probabilities for this hypothesis which neither close to 1 nor to 0. We believe that this sort of information can be useful in clinical follow-up measurements when trying to monitor possible changes in the properties of the tissue.

Zoom In Zoom Out Reset image size

Usually, Bayesian inferences require extensive numerical calculations, typically performed by Markov chain Monte Carlo methods (MCMC) (Robert and Casella 2013). Results achieved by MCMC methods are stochastic in nature, and need some form of convergence checks. An advantage of the proposed approach is that its results are deterministic and can be calculated efficiently. Furthermore, as input only the data (i.e. sequence of magnetization at every voxel) and the dictionary are required, which provides a convenient interface for applications.

Reliable quantification of uncertainties in quantitative magnetic imaging is beneficial as it allows a meaningful assessment of single results. A Bayesian uncertainty quantification has been proposed for estimates obtained by a dictionary matching method. The proposed uncertainty quantification has been tested through its application to measurements of a suitably designed phantom. Comparison with additional reference measurements, together with consistency checks, as well as results obtained for simulated data, show that the proposed approach yields a reasonable uncertainty quantification. The proposed uncertainty quantification has also been applied to an in vivo measurement. The Bayesian approach offers the possibility to also derive probability statements which can increase confidence in conclusions drawn from MRF results. While Bayesian inference typically requires extensive calculations, our analytical derivations allowed for a calculation of uncertainties based on the pre-calculated dictionary. As a consequence, uncertainty calculations can be carried out efficiently. We believe that applications of the dictionary-based MRF can benefit from the proposed uncertainty quantification scheme.

This work was supported by EMPIR project 18HLT05 QUIERO. This project has received funding from the EMPIR program co-financed by the Participating States and from the European Union's Horizon 2020 research and innovation program.

The dictionary used for the dictionary matching method is based on a physical model describing the magnetization process of the MRI sequence. We will use the Bloch equations to derive a three dimensional linear system of equation for the employed FISP sequence, for more details we refer to Hargreaves et al (2001). In the FISP sequence a dephasing gradient is applied in the middle of each repetition time TR. The signal is therefore an average of the spins in one voxel. We will simulate this by taking Nf different spins in one voxel. The following Matlab^® code is used to explain the implemented sequence. The flip angles will be denoted by flip, echo times by TE and the inversion time by TI.

The matrix Mⁱ containing the real and imaginary part of the modeled magnetization course ${m}^{i}({T}_{1}^{i},{T}_{2}^{i},{\rho }^{i})$ is introduced in (2.8). The modeled magnetization course is given as ${m}^{i}({T}_{1}^{i},{T}_{2}^{i},{\rho }^{i})=({\mathfrak{R}}({\rho }^{i})+i{\mathfrak{I}}({\rho }^{i}))({\mathfrak{R}}({\tilde{m}}^{i})\,+i{\mathfrak{I}}({\tilde{m}}^{i}))$ where ${\tilde{m}}^{i}$ is the part of the whole magnetization without the linear factor ρⁱ .

It follows that the real and imaginary part of mⁱ are ${\mathfrak{R}}({m}^{i}({T}_{1}^{i},{T}_{2}^{i},{\rho }^{i}))={\mathfrak{R}}({\rho }^{i}){\mathfrak{R}}({\tilde{m}}^{i})-{\mathfrak{I}}({\rho }^{i}){\mathfrak{I}}({\tilde{m}}^{i})$ and ${\mathfrak{I}}({m}^{i}({T}_{1}^{i},{T}_{2}^{i},{\rho }^{i}))={\mathfrak{I}}({\rho }^{i}){\mathfrak{R}}({\tilde{m}}^{i})+{\mathfrak{R}}({\rho }^{i}){\mathfrak{I}}({\tilde{m}}^{i})$. Mⁱ is defined as

$$\begin{eqnarray*}&&{\left({\mathfrak{R}}{\left({m}^{i}({T}_{1}^{i},{T}_{2}^{i},{\rho }^{i})\right)}^{T},{\mathfrak{I}}({m}^{i}{\left({T}_{1}^{i},{T}_{2}^{i},{\rho }^{i}\right)}^{T})\right)}^{T}={M}^{i}{\left({\mathfrak{R}}({\rho }^{i}),{\mathfrak{I}}({\rho }^{i})\right)}^{T}.\end{eqnarray*}$$

This leads to

$$\begin{eqnarray*}\begin{array}{c}{M}^{i}=\left(\begin{array}{cc}{\mathfrak{R}}({\tilde{m}}^{i}({T}_{1}^{i},{T}_{2}^{i})) & -{\mathfrak{I}}({\tilde{m}}^{i}({T}_{1}^{i},{T}_{2}^{i}))\\ {\mathfrak{I}}({\tilde{m}}^{i}({T}_{1}^{i},{T}_{2}^{i})) & {\mathfrak{R}}({\tilde{m}}^{i}({T}_{1}^{i},{T}_{2}^{i}))\end{array}\right).\end{array}\end{eqnarray*}$$

In section 2 the physical model is introduced. We highlighted that this model can easily be extended to an acquisition with several receiver coils. Let $\{{C}_{i}\}{}_{i=1}^{{n}_{c}}$ be the n_c different $\sqrt{N}\times \sqrt{N}$ coil sensitivity matrices. We define $\tilde{{C}_{i}}={\rm{diag}}(\{{\left({C}_{i}\right)}_{j}\}{}_{j=1}^{N})\in {{\mathbb{R}}}^{N\,\times \,N}$. Equation (2.1) will then change to

$$\begin{eqnarray}&&{z}_{{li}}={F}_{l}({\tilde{C}}_{i}{m}_{l}),\,\,l=1,\,\ldots ,\,L,\end{eqnarray} \tag{ A.1 }$$

and 2.2 will become

$$\begin{eqnarray}&&{\widehat{m}}_{l}=\displaystyle \frac{1}{{n}_{c}}\displaystyle \sum _{i=1}^{{n}_{c}}{\tilde{C}}_{i}^{-1}{F}_{l}^{\dagger }({z}_{{li}}),\,\,l=1,\,\ldots ,\,L.\end{eqnarray} \tag{ A.2 }$$

The entries in ${\tilde{C}}_{i}$ can be very close to zero or even be zero, therefore it is better to work with the complex conjugate which is denoted by $\overline{{\tilde{C}}_{i}}$. We will then get

$$\begin{eqnarray*}&&{\widehat{m}}_{l}={\left(\overline{{\tilde{C}}_{i}}{\tilde{C}}_{i}\right)}^{-1}\overline{{\tilde{C}}_{i}}{F}_{l}^{\dagger }({z}_{{li}}),\,\,l=1,\,\ldots ,\,L,i=1,\,\ldots ,\,{n}_{c}.\end{eqnarray*}$$

and

$$\begin{eqnarray*}&&(\overline{{\tilde{C}}_{i}}{\tilde{C}}_{i}){\widehat{m}}_{l}=\overline{{\tilde{C}}_{i}}{F}_{l}^{\dagger }({z}_{{li}}),\,\,l=1,\,\ldots ,\,L.i=1,\,\ldots ,\,{n}_{c}.\end{eqnarray*}$$

Equation (2.2) will then become

$$\begin{eqnarray}&&{\widehat{m}}_{l}={\left(\displaystyle \sum _{i=1}^{{n}_{c}}\overline{{\tilde{C}}_{i}}{\tilde{C}}_{i}\right)}^{-1}\displaystyle \sum _{i=1}^{{n}_{c}}\overline{{\tilde{C}}_{i}}{F}_{l}^{\dagger }({z}_{{li}}),\,\,l=1,\,\ldots ,\,L.\end{eqnarray} \tag{ A.3 }$$

All the other derivations follow easily.

The dictionary matching solution (see (Davies et al 2014)) for ${({\hat{T}}_{1})}_{i}$ and ${({\hat{T}}_{2})}_{i}$ in voxel i is chosen to be the entry ${\hat{k}}_{i}$ of the dictionary ${\{{D}_{k}\}}_{k=1}^{K}$ which maximizes

$$\begin{eqnarray*}&&{\hat{k}}_{i}={\mathrm{argmax}}_{k}\displaystyle \frac{| \lt {\hat{m}}^{i},{D}_{k}\gt | }{\parallel {D}_{k}\parallel }.\end{eqnarray*}$$

The complex proton density is then calculated as

$$\begin{eqnarray*}&&{\hat{\rho }}_{i}=\displaystyle \frac{\lt {\hat{m}}^{i},{D}_{{\hat{k}}_{i}}\gt }{\parallel {D}_{{\hat{k}}_{i}}{\parallel }^{2}}.\end{eqnarray*}$$

In the following, the reference to a specific voxel i is omitted as the treatment is the same for all voxels.

The dictionary is based on the physical model m(T₁, T₂, ρ) without the linear factor ρ, which we will denote by $\tilde{m}({T}_{1},{T}_{2})$. The above equations are therefore equivalent to

$$\begin{eqnarray*}\begin{array}{rcl}({\hat{T}}_{1},{\hat{T}}_{2}) & = & {\mathrm{argmax}}_{{T}_{1},{T}_{2}}\displaystyle \frac{| \lt \hat{m},\tilde{m}({T}_{1},{T}_{2})\gt | }{\parallel \tilde{m}({T}_{1},{T}_{2})\parallel }\\ \hat{\rho } & = & \displaystyle \frac{\lt \hat{m},\tilde{m}({\hat{T}}_{1},{\hat{T}}_{2})\gt }{\parallel \tilde{m}({\hat{T}}_{1},{\hat{T}}_{2}){\parallel }^{2}},\end{array}\end{eqnarray*}$$

provided that the resolution of dictionary is infinitely small.

Let ${\tilde{m}}_{R}$ and ${\tilde{m}}_{I}$ denote the real and imaginary part of $\tilde{m}({T}_{1},{T}_{2})$, and ${\hat{m}}_{R}$ and ${\hat{m}}_{I}$ real and imaginary part of $\hat{m}$. It follows

$$\begin{eqnarray*}\begin{array}{rcl}({\hat{T}}_{1},{\hat{T}}_{2}) & = & {\mathrm{argmax}}_{{T}_{1},{T}_{2}}\displaystyle \frac{\sqrt{{\left({\hat{m}}_{R}^{T}{\tilde{m}}_{R}+{\hat{m}}_{I}^{T}{\tilde{m}}_{I}\right)}^{2}+{\left(-{\hat{m}}_{R}^{T}{\tilde{m}}_{I}+{\hat{m}}_{I}^{T}{\tilde{m}}_{R}\right)}^{2}}}{\sqrt{{\tilde{m}}_{R}^{T}{\tilde{m}}_{R}+{\tilde{m}}_{I}^{T}{\tilde{m}}_{I}}}\\ & = & {\mathrm{argmax}}_{{T}_{1},{T}_{2}}\displaystyle \frac{\sqrt{h}}{\sqrt{{\tilde{m}}_{R}^{T}{\tilde{m}}_{R}+{\tilde{m}}_{I}^{T}{\tilde{m}}_{I}}}\\ & = & {\mathrm{argmax}}_{{T}_{1},{T}_{2}}\displaystyle \frac{h}{{\tilde{m}}_{R}^{T}{\tilde{m}}_{R}+{\tilde{m}}_{I}^{T}{\tilde{m}}_{I}},\end{array}\end{eqnarray*}$$

for $h={\left({\hat{m}}_{R}^{T}{\tilde{m}}_{R}\right)}^{2}+{\left({\hat{m}}_{I}^{T}{\tilde{m}}_{I}\right)}^{2}+{\left({\hat{m}}_{R}^{T}{\tilde{m}}_{I}\right)}^{2}+{\left({\hat{m}}_{I}^{T}{\tilde{m}}_{R}\right)}^{2}+2({\hat{m}}_{R}^{T}{\tilde{m}}_{R}{\hat{m}}_{I}^{T}{\tilde{m}}_{I})-2({\hat{m}}_{R}^{T}{\tilde{m}}_{I}{\hat{m}}_{I}^{T}{\tilde{m}}_{R})$ and

$$\begin{eqnarray}\begin{array}{rcl}\hat{\rho } & = & \displaystyle \frac{\lt {\hat{m}}_{R}+i{\hat{m}}_{I},{\tilde{m}}_{R}+i{\tilde{m}}_{I}\gt }{\parallel {\tilde{m}}_{R}+i{\tilde{m}}_{I}{\parallel }^{2}}\\ & = & \displaystyle \frac{{\hat{m}}_{R}^{T}{\tilde{m}}_{R}+{\hat{m}}_{I}^{T}{\tilde{m}}_{I}}{{\tilde{m}}_{R}^{T}{\tilde{m}}_{R}+{\tilde{m}}_{I}^{T}{\tilde{m}}_{I}}+i\displaystyle \frac{{\hat{m}}_{I}^{T}{\tilde{m}}_{R}-{\hat{m}}_{R}^{T}{\tilde{m}}_{I}}{{\tilde{m}}_{R}^{T}{\tilde{m}}_{R}+{\tilde{m}}_{I}^{T}{\tilde{m}}_{I}}.\end{array}\end{eqnarray} \tag{ B.1 }$$

We will now consider the minimization of (2.5). It is

$$\begin{eqnarray*}\begin{array}{rcl}\hat{\mu } & = & {\mathrm{argmin}}_{\mu }\parallel \hat{m}-M\mu {\parallel }^{2}\\ & = & {\mathrm{argmin}}_{\mu }-2{\hat{m}}^{T}M\mu +{\mu }^{T}{M}^{T}M\mu \\ & = & {\left({M}^{T}M\right)}^{-1}{M}^{T}\hat{m}\\ & = & {\left(\left(\begin{array}{cc}{\tilde{m}}_{R}^{T} & {\tilde{m}}_{I}^{T}\\ -{\tilde{m}}_{I}^{T} & {\tilde{m}}_{R}^{T}\end{array}\right)\left(\begin{array}{cc}{\tilde{m}}_{R} & -{\tilde{m}}_{I}\\ {\tilde{m}}_{I} & {\tilde{m}}_{R}\end{array}\right)\right)}^{-1}\left(\begin{array}{cc}{\tilde{m}}_{R}^{T} & {\tilde{m}}_{I}^{T}\\ -{\tilde{m}}_{I}^{T} & {\tilde{m}}_{R}^{T}\end{array}\right)\left(\begin{array}{c}{\hat{m}}_{R}\\ {\hat{m}}_{I}\end{array}\right)\\ & = & \displaystyle \frac{1}{{\tilde{m}}_{R}^{T}{\tilde{m}}_{R}+{\tilde{m}}_{I}^{T}{\tilde{m}}_{I}}\left(\begin{array}{c}{\hat{m}}_{R}^{T}{\tilde{m}}_{R}+{\hat{m}}_{I}^{T}{\tilde{m}}_{I}\\ -{\hat{m}}_{R}^{T}{\tilde{m}}_{I}+{\hat{m}}_{I}^{T}{\tilde{m}}_{R}\end{array}\right),\end{array}\end{eqnarray*}$$

which corresponds to the dictionary matching solution B.1.

Let $c={\left({\hat{m}}_{R}^{T}{\tilde{m}}_{R}+{\hat{m}}_{I}^{T}{\tilde{m}}_{I},\ -{\hat{m}}_{R}^{T}{\tilde{m}}_{I}+{\hat{m}}_{I}^{T}{\tilde{m}}_{R}\right)}^{T}$, then

$$\begin{eqnarray*}\begin{array}{rcl}{\mathrm{argmin}}_{{T}_{1},{T}_{2}}\parallel \hat{m}-M\hat{\mu }{\parallel }^{2} & = & {\mathrm{argmin}}_{{T}_{1},{T}_{2}}-2\lt \hat{m},M\hat{\mu }\gt +\lt M\hat{\mu },M\hat{\mu }\gt \\ & = & {\mathrm{argmax}}_{{T}_{1},{T}_{2}}2\lt \hat{m},M\hat{\mu }\gt -\lt M\hat{\mu },M\hat{\mu }\gt \\ & = & {\mathrm{argmax}}_{{T}_{1},{T}_{2}}2{\hat{m}}^{T}M\hat{\mu }-{\hat{\mu }}^{T}{M}^{T}M\hat{\mu }\\ & = & {\mathrm{argmax}}_{{T}_{1},{T}_{2}}\displaystyle \frac{2}{{\tilde{m}}_{R}^{T}{\tilde{m}}_{R}+{\tilde{m}}_{I}^{T}{\tilde{m}}_{I}}{c}^{T}c\\ & & -\,\displaystyle \frac{1}{{\left({\tilde{m}}_{R}^{T}{\tilde{m}}_{R}+{\tilde{m}}_{I}^{T}{\tilde{m}}_{I}\right)}^{2}}{c}^{T}\left(\begin{array}{cc}{\tilde{m}}_{R}^{T}{\tilde{m}}_{R}+{\tilde{m}}_{I}^{T}{\tilde{m}}_{I} & 0\\ 0 & {\tilde{m}}_{R}^{T}{\tilde{m}}_{R}+{\tilde{m}}_{I}^{T}{\tilde{m}}_{I}\end{array}\right)c\\ & = & {\mathrm{argmax}}_{{T}_{1},{T}_{2}}\displaystyle \frac{1}{{\tilde{m}}_{R}^{T}{\tilde{m}}_{R}+{\tilde{m}}_{I}^{T}{\tilde{m}}_{I}}{c}^{T}c\\ & = & {\mathrm{argmax}}_{{T}_{1},{T}_{2}}\displaystyle \frac{h}{{\tilde{m}}_{R}^{T}{\tilde{m}}_{R}+{\tilde{m}}_{I}^{T}{\tilde{m}}_{I}}.\end{array}\end{eqnarray*}$$

Hence, the dictionary matching method yields estimates which are given by minimizing (2.5).

In order to derive the posterior distribution in (2.15), we need two prior statements:

For $x\in {{\mathbb{R}}}^{p},\mu \in {{\mathbb{R}}}^{p},\nu \gt 1$ and symmetric, positive-definite ${\rm{\Sigma }}\in {{\mathbb{R}}}^{p\times p}$ we have

$$\begin{eqnarray*}&&\int \displaystyle \frac{1}{{\rm{\det }}{\left({\rm{\Sigma }}\right)}^{1/2}}{\left(1+\displaystyle \frac{1}{\nu }{\left(x-\mu \right)}^{T}{{\rm{\Sigma }}}^{-1}(x-\mu )\right)}^{-(\nu +p)/2}{dx}\propto 1\end{eqnarray*}$$

and therefore

$$\begin{eqnarray*}&&\int {\left(1+\displaystyle \frac{1}{\nu }{\left(x-\mu \right)}^{T}{{\rm{\Sigma }}}^{-1}(x-\mu )\right)}^{-(\nu +p)/2}{dx}\propto {\rm{\det }}{\left({\rm{\Sigma }}\right)}^{1/2}.\end{eqnarray*}$$

This follows from the PDF of the multivariate t-distribution.

For $A\in {{\mathbb{R}}}^{2\times 2}$ and $a\in {\mathbb{R}}$ it is

$$\begin{eqnarray*}\begin{array}{rcl}{\rm{\det }}{\left({\left(\displaystyle \frac{1}{{a}^{2}}A\right)}^{-1}\right)}^{1/2} & = & {\rm{\det }}{\left(\displaystyle \frac{1}{{a}^{2}}A\right)}^{-1/2}={\left({\left(\displaystyle \frac{1}{{a}^{2}}\right)}^{2}{\rm{\det }}(A)\right)}^{-1/2}\\ & = & {a}^{2}{\rm{\det }}{\left(A\right)}^{-1/2}.\end{array}\end{eqnarray*}$$

We then get for $\hat{\mu }=\mu ({T}_{1},{T}_{2})={\left[{M}^{T}M\right]}^{-1}{M}^{T}y$:

$$\begin{eqnarray}\begin{array}{rcl}\pi ({T}_{1},{T}_{2}| y) & \propto & {\displaystyle \int }_{-\infty }^{\infty }{\displaystyle \int }_{0}^{\infty }\pi ({T}_{1},{T}_{2},\mu ,{\sigma }^{2}| y){\rm{d}}{\sigma }^{2}{\rm{d}}\mu \\ & \propto & {\displaystyle \int }_{-\infty }^{\infty }{\displaystyle \int }_{0}^{\infty }{{\bf{1}}}_{{\rm{\Omega }}}({T}_{1},{T}_{2})\displaystyle \frac{1}{{\sigma }^{2L+2}}\exp \left(-\displaystyle \frac{1}{2{\sigma }^{2}}\parallel y-M\mu {\parallel }^{2}\right){\rm{d}}{\sigma }^{2}{\rm{d}}\mu \\ & = & {{\bf{1}}}_{{\rm{\Omega }}}({T}_{1},{T}_{2}){\displaystyle \int }_{-\infty }^{\infty }{\displaystyle \int }_{0}^{\infty }\displaystyle \frac{1}{{\tau }^{L+1}}\exp \left(-\displaystyle \frac{1}{2\tau }\parallel y-M\mu {\parallel }^{2}\right){\rm{d}}\tau {\rm{d}}\mu \\ & = & {{\bf{1}}}_{{\rm{\Omega }}}({T}_{1},{T}_{2}){\displaystyle \int }_{-\infty }^{\infty }{\displaystyle \int }_{0}^{\infty }\displaystyle \frac{1}{{\left(u\parallel y-M\mu {\parallel }^{2}\right)}^{L+1}}\parallel y-M\mu {\parallel }^{2}\exp \left(-\displaystyle \frac{1}{2u}\right){\rm{d}}u{\rm{d}}\mu \\ & = & {{\bf{1}}}_{{\rm{\Omega }}}({T}_{1},{T}_{2}){\displaystyle \int }_{-\infty }^{\infty }{\displaystyle \int }_{0}^{\infty }\displaystyle \frac{1}{({u}^{L+1})}{\left(\parallel y-M\mu {\parallel }^{2}\right)}^{-L}\exp \left(-\displaystyle \frac{1}{2u}\right){\rm{d}}u{\rm{d}}\mu \\ & \propto & {{\bf{1}}}_{{\rm{\Omega }}}({T}_{1},{T}_{2}){\displaystyle \int }_{-\infty }^{\infty }\displaystyle \frac{1}{{\left(\parallel y-M\mu {\parallel }^{2}\right)}^{L}}{\rm{d}}\mu \\ & = & {{\bf{1}}}_{{\rm{\Omega }}}({T}_{1},{T}_{2}){\displaystyle \int }_{-\infty }^{\infty }\displaystyle \frac{1}{{\left(\parallel y-M\hat{\mu }{\parallel }^{2}+(\mu -\hat{\mu }){M}^{T}M(\mu -\hat{\mu })\right)}^{L}}{\rm{d}}\mu \\ & = & {{\bf{1}}}_{{\rm{\Omega }}}({T}_{1},{T}_{2})\displaystyle \frac{1}{{\left(\parallel y-M\hat{\mu }{\parallel }^{2}\right)}^{L}}{\displaystyle \int }_{-\infty }^{\infty }\displaystyle \frac{1}{{\left(1+(\mu -\hat{\mu })\tfrac{{M}^{T}M}{\parallel y-M\hat{\mu }{\parallel }^{2}}(\mu -\hat{\mu })\right)}^{L}}{\rm{d}}\mu \end{array}\end{eqnarray} \tag{ B.2 }$$

$$\begin{eqnarray}\begin{array}{rcl} & \propto & {{\bf{1}}}_{{\rm{\Omega }}}({T}_{1},{T}_{2})\displaystyle \frac{1}{{\left(\parallel y-M\hat{\mu }{\parallel }^{2}\right)}^{L}}\parallel y-M\hat{\mu }{\parallel }^{2}\displaystyle \frac{1}{{\rm{\det }}{\left({M}^{T}M\right)}^{1/2}}\\ & = & {{\bf{1}}}_{{\rm{\Omega }}}({T}_{1},{T}_{2})\displaystyle \frac{1}{{\left(\parallel y-M\hat{\mu }{\parallel }^{2}\right)}^{L-1}}\displaystyle \frac{1}{{\rm{\det }}{\left({M}^{T}M\right)}^{1/2}}.\end{array}\end{eqnarray} \tag{ B.3 }$$

The posterior (B.3) exists if we make the following two assumptions: ${\rm{\det }}({M}^{T}M)\gt {\delta }_{1}\gt 0$ and ${\min }_{{T}_{1},{T}_{2}}\parallel y-M\hat{\mu }{\parallel }^{2}\gt {\delta }_{2}\gt 0$. Together with the compact set of T₁ and T₂, the propriety of the posterior is then guaranteed.

A similar derivation with marginalization w.r.t. σ² yields the expression over which is integrated in (B.2), from it follows that the conditional posterior for μ is that given in (2.17).

The joint posterior distribution for T₁ and T₂ is given in (2.15). This function is calculated on a grid of T₁ and T₂. By taking numerical integrals, the marginal posterior distribution for T₁ and T₂ can then be approximated. Note that if the stepsize of the grid is infinitely small, the integral is exactly calculated. We used a trapezoidal method with unit spacing given via

$$\begin{eqnarray*}\begin{array}{rcl}{\displaystyle \int }_{a}^{b}f(x){\rm{d}}x & \approx & \displaystyle \frac{b-a}{2N}\displaystyle \sum _{n=1}^{N}(f({x}_{n})+f({x}_{n+1}))\\ & = & \displaystyle \frac{b-a}{2N}\left(f({x}_{1})+2f({x}_{2})+\ldots +2f({x}_{N})+f({x}_{N+1})\right),\end{array}\end{eqnarray*}$$

where the scalar value $\tfrac{b-a}{N}$ is equal to one. This is implemented in the Matlab^® function trapz.