# Scalable magnetic resonance fingerprinting: Incremental inference of high dimensional elliptical mixtures from large data volumes

Geoffroy Oudoumanessah, Thomas Coudert, Carole Lartizien, Michel Dojat, Thomas Christen, Florence Forbes

► To cite this version:

Geoffroy Oudoumanessah, Thomas Coudert, Carole Lartizien, Michel Dojat, Thomas Christen, et al.. Scalable magnetic resonance fingerprinting: Incremental inference of high dimensional elliptical mixtures from large data volumes. Annals of Applied Statistics, In press. hal-04836864v2

HAL Id: hal-04836864

https://hal.science/hal-04836864v2

Submitted on 7 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
SCALABLE MAGNETIC RESONANCE FINGERPRINTING: INCREMENTAL
INFERENCE OF HIGH-DIMENSIONAL ELLIPTICAL MIXTURES FROM
LARGE DATA VOLUMES

BY GEOFFROY OUDOUMANESSAH<sup>1,3,2,a</sup> ID, THOMAS COUDERT<sup>3,d</sup> ID
CAROLE LARTIZIEN<sup>2,c</sup> ID MICHEL DOJAT<sup>1,3,e</sup> ID
THOMAS CHRISTEN<sup>3,f</sup> ID AND FLORENCE FORBES<sup>1,b</sup> ID

<sup>1</sup>Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, LJK, 38000 Grenoble, France, <sup>a</sup>geoffroy.oudoumanessah@inria.fr;
<sup>b</sup>florence.forbes@inria.fr

<sup>2</sup>Univ. Lyon, CNRS, Inserm, INSA Lyon, UCBL, CREATIS, UMR5220, U1294, F-69621, Villeurbanne, France,
<sup>c</sup>carole.lartizien@creatis.insa-lyon.fr

<sup>3</sup>Univ. Grenoble Alpes, Inserm U1216, CHU Grenoble Alpes, Grenoble Institut des Neurosciences, 38000 Grenoble, France,
<sup>d</sup>thomas.coudert@inserm.fr; <sup>e</sup>michel.dojat@inserm.fr; <sup>f</sup>thomas.christen@univ-grenoble-alpes.fr

Magnetic Resonance Fingerprinting (MRF) is an emerging technology with the potential to revolutionize radiology and medical diagnostics. In comparison to traditional magnetic resonance imaging (MRI), MRF enables the rapid, simultaneous, non-invasive acquisition and reconstruction of multiple tissue parameters, paving the way for novel diagnostic techniques. In the original *matching* approach, reconstruction is based on the search for the best matches between *in vivo* acquired signals and a dictionary of high-dimensional simulated signals (fingerprints) with known tissue properties. A critical and limiting challenge is that the size of the simulated dictionary increases exponentially with the number of parameters, leading to an extremely costly matching. In this work, we propose to address this scalability issue by considering probabilistic mixtures of high-dimensional elliptical distributions, to learn more efficient dictionary representations. Mixture components are modelled as flexible elliptical shapes in low dimensional subspaces. They are exploited to cluster similar signals and reduce their dimension locally cluster-wise limiting information loss. To estimate such a mixture model, we provide a new incremental algorithm capable of handling large numbers of signals, allowing us to go far beyond the hardware limitations encountered by standard implementations. We demonstrate, on simulated and real data, that our method effectively manages large volumes of MRF data with maintained accuracy. It offers a more efficient solution for accurate tissue characterization and significantly reduces the computational burden, making the clinical application of MRF more practical and accessible.

**1. Introduction.** Traditional Magnetic Resonance (MR) imaging relies on an analytical resolution of dynamical equations using conventional tuning of the MR hardware through sequences of pulses, each characterized by different values of parameters such as the flip angle and repetition time. Standard quantitative MRI (qMRI) methods are based on a single sequence for a single parameter measurement at a time. This leads to high scan times for multi-parametric protocols as each parameter estimate involves one MR sequence. A recent approach named Magnetic Resonance Fingerprinting (MRF, Ma et al. (2013)) has been developed to overcome these limitations. The MRF protocol involves fast undersampled acquisitions with time-varying parameters defining the MRF sequence that produces temporal signal evolutions (named *fingerprints*) in each voxel. In the original proposal, a dictionary

*Keywords and phrases:* Dimension reduction, Clustering, Incremental learning, High dimensional mixture models, Elliptical distributions, Expectation Maximization algorithm.
search approach is used to compare the *in vivo* fingerprints with millions of numerical simu-
lations of MR signals for which the associated parameters are known. These millions of sim-
ulated signals compose the so-called dictionary. The values of the parameters corresponding
to the closest simulated signals or *matches* are then assigned to the associated *in vivo* vox-
els, allowing the simultaneous reconstruction of multiple quantitative maps (images) from
extremely undersampled raw images, using only one single sequence, thus saving consid-
erable acquisition time (Poorman et al., 2020; McGivney et al., 2020). From this, standard
relaxometry MRF allows reconstructing parameter maps for relaxation times T₁ and T₂, over
the whole human brain (1 mm³ spatial resolution) in 3 min (Ye et al., 2017; Gu et al., 2018)
compared to 30 min for a standard T1/T2 exam. Moreover, the flexibility of the numerical
simulations enables correction of system imperfections as well as some patient motions by
including them in the model and post-processing pipelines (Bipin Mehta et al., 2019). Thus,
MRF could be a game changer for emergency patients who need to complete exams in a few
minutes. The power of the MRF approach is not limited to the estimation of relaxation times,
in theory, it allows the measurement of any parameter that influences nuclear magnetization
(e.g., microvascular networks), and could be added to the simulation model (Wang et al.,
2019a; Coudert et al., 2025a,b). However, increasing the number of estimated parameters,
even moderately, induces the design of more complex sequences and increased reconstruc-
tion times, from hours to days. This limits the clinical application of high-dimensional MRF
and necessitates the development of innovative processing methods. Consequently, a signif-
icant focus is on improving MRF reconstruction methods, as reviewed by Tippareddy et al.
(2021) and Monga et al. (2024). In this work, we propose to focus on reducing the reconstruc-
tion times of MRF when more than the main two parameters, T₁, T₂, are involved, including
the addition of δf the frequency offset, the sensitivity of the magnetic field B₁, cerebral blood
volume (CBV), and microvascular geometry (e.g., vessel radius denoted as R).
MRF reconstruction is first recast as an inverse problem that can be solved using differ-
ent approaches as recalled in the next section. All approaches make use, at some stage, of
a dictionary of simulated pairs (parameters, signal), which represent our knowledge of the
link between tissue parameters and MR time series, through a so-called direct or forward
model. The dictionary is then either used to learn an inversion operator, from signal to tissue
parameters, or to search for the best fits between observed signals and simulated ones. The
approaches' scalability relies thus greatly on their ability to extract efficiently the informa-
tion encoded in simulations. Efficiency has different aspects: for search-type methods, the
dictionary should not be too big, while for learning-based methods, the dictionary should
be informative enough. These potentially opposite requirements call for efficient representa-
tions of simulated data. We propose to explore a divide-and-conquer strategy by introducing
the framework of High Dimensional Mixtures of Elliptical Distributions (HD-MED). Proba-
bilistic mixtures of high-dimensional elliptical distributions allow us to learn more efficient
dictionary representations. Mixture components are modeled as flexible elliptical shapes in
low-dimensional subspaces. They are exploited to cluster similar signals and reduce their
dimension locally, at the cluster level, limiting information loss. To estimate such a mixture
model, we provide a new incremental algorithm capable of handling large numbers of signals,
allowing us to go far beyond the hardware limitations encountered with standard implemen-
tations. We demonstrate, on simulated and real data, that our method effectively manages
large volumes of MRF data with maintained accuracy. It offers a more efficient solution for
accurate tissue characterization and significantly reduces computational burden, making the
clinical application of MRF more practical and accessible.

MRF reconstruction is first recast as an inverse problem that can be solved using differ-
ent approaches as recalled in the next section. All approaches make use, at some stage, of
a dictionary of simulated pairs (parameters, signal), which represent our knowledge of the
link between tissue parameters and MR time series, through a so-called direct or forward
model. The dictionary is then either used to learn an inversion operator, from signal to tissue
parameters, or to search for the best fits between observed signals and simulated ones. The
approaches’ scalability relies thus greatly on their ability to extract efficiently the informa-
tion encoded in simulations. Efficiency has different aspects: for search-type methods, the
dictionary should not be too big, while for learning-based methods, the dictionary should
be informative enough. These potentially opposite requirements call for efficient representa-
tions of simulated data. We propose to explore a divide-and-conquer strategy by introducing
the framework of High Dimensional Mixtures of Elliptical Distributions (HD-MED). Proba-
bilistic mixtures of high-dimensional elliptical distributions allow us to learn more efficient
dictionary representations. Mixture components are modeled as flexible elliptical shapes in
low-dimensional subspaces. They are exploited to cluster similar signals and reduce their
dimension locally, at the cluster level, limiting information loss. To estimate such a mixture
model, we provide a new incremental algorithm capable of handling large numbers of signals,
allowing us to go far beyond the hardware limitations encountered with standard implemen-
tations. We demonstrate, on simulated and real data, that our method effectively manages
large volumes of MRF data with maintained accuracy. It offers a more efficient solution for
accurate tissue characterization and significantly reduces computational burden, making the
clinical application of MRF more practical and accessible.

In the rest of the paper, we first recall in Section 2 the two main types of approaches
that have been investigated in the literature on MRF reconstruction. We then review related
work on efficient dictionary representations and specify our contributions in Section 3. The
mixtures of elliptical distributions are presented in Section 4 with their use for dimension reduction in Section 5. The incremental algorithm proposed for the mixture estimation is presented in Section 6 and illustrated on MRF in Section 7.

**2. MRF reconstruction as an inverse problem.** An inverse problem refers to a situation where one aims to determine the causes of a phenomenon based on experimental observations. In MRF, the goal is to infer a set of tissue characteristics or parameter values from an observed signal, identifying those that best explain the signal. Compared to standard MRI, which relies on sequential measurements with fixed parameters, MRF acquires data using a sequence of pseudo-random parameters (e.g., varying repetition times or flip angles). This process generates tissue-dependent time-series signals, where each time-series y - known as a fingerprint - is assigned to a specific 3D location (voxel) in the brain. Such a resolution generally starts with modeling the phenomenon under consideration, which is called the "direct" or "forward problem." It is generally assumed that at least the numerical evaluation of the forward model is available because experts have designed equations that can be solved either analytically or numerically. The most common use of the forward model is via a simulator that allows the creation of a database $D_f$, usually referred to as a dictionary, of N signals $y_1, \dots, y_N$ with $y_i \in \mathcal{Y} \subset \mathbb{R}^M$, generated (stored or computed on the fly) by running the theoretical (physical) model f for many different tissue parameter values $t_1, \dots, t_N$ with $t_i \in \mathcal{T} \subset \mathbb{R}^L$ and $y_i = f(t_i)$. The typical dimension of each $y_i$ considered in this paper is $M = 260$. The generated tissue parameter values then only partially represent the full space $\mathcal{T}$ of possible values and correspond to a discrete grid in the space. In this context, we can distinguish two types of methods, referred to below as optimization and learning approaches.

**2.1. Optimization vs learning (or regression) approaches.** Optimization approaches include the most used method in MRF, namely dictionary matching, often called, in other domains, grid search, look-up-table or k-nearest neighbors. They consist of minimizing over parameters **t** a merit function **d** expressing the similarity between the observed signal $y_{obs}$ and simulated signal $f(\mathbf{t})$,

$$ (1) \qquad \hat{\mathbf{t}} \in \arg \min_{\mathbf{t} \in \mathcal{T}} d(\mathbf{y}_{obs}, f, \mathbf{t}). $$

Typically, $d(\mathbf{y}_{obs}, \mathbf{f}, \mathbf{t}) = d(\mathbf{y}_{obs}, f(\mathbf{t})) = \|\mathbf{y}_{obs} - f(\mathbf{t})\|^2$. Solutions are searched in the full $\mathcal{T}$ space but solutions could also be penalized as done in grid search methods. In grid search, the previous full search is replaced by a simpler look-up or matching operation making use of $D_f$, often created beforehand. The search space is significantly reduced from a continuous space $\mathcal{T}$ to a discrete and finite $D_f$. The speed gain is significant in comparison to traditional optimization methods as retrieving a value from memory is often faster than undergoing an expensive usually iterative computation. Their disadvantage is the instability of solutions. Many questions remain on how to choose the merit function, how many $y_n$ in the look-up table have to be kept to estimate parameters, how to choose the look-up table, etc. When the number of parameters is small, grid search is suitable and can provide very good predictions. However, for even moderate numbers of parameters, the required number of elements in the dictionary renders grid search either intractable or inaccurate. The technique is not amortized, for each new $y_{obs}$, we have to compute the matching score $d(\mathbf{y}_{obs}, \mathbf{y}_n)$ for all $y_n$ in the dictionary. When the dimension of $\mathbf{t}$ increases, the dictionary (and $N$) has to be larger too for better accuracy. The computation of $N$ matching scores can become time-consuming.

Regression or learning methods are more efficient in that sense and usually have better amortization properties. These methods have the advantage of providing tractable solutions in the case of massive inversions of high-dimensional data. The main principle is to transfer the computational cost and time from individual pointwise predictions to the learning of a
global inverse operator from $\mathcal{D}_f$. The advantage is that once the operator is learned, it can be
used, at negligible cost, for very large numbers of new signals. Then, the nature of the inverse
operator needs to be specified. Traditional learning or regression methods are not specifically
designed for high-dimensional data but there has been a large literature covering this case,
see e.g. Giraud (2014) for a review. In MRF, a popular approach from Cohen, Zhu and Rosen
(2018), uses a four-layer fully connected network (DRONE) to learn the dictionary signals
for reconstructing $T_1$ and $T_2$ parameters. Even though DRONE provides good results, this
perceptron-based model loses the temporal coherence of the signal. Recently, Cabini et al.
(2024) proposed a Recurrent Neural Network (RNN) with long-short term memory (LSTM)
blocks, which yields better reconstruction results for $T_1$ and $T_2$ with more robustness towards
noisy acquisitions. However, when more parameters need to be estimated, such as vessel
oxygenation or radius (Christen, Bolar and Zaharchuk, 2013; Christen et al., 2014; Coudert
et al., 2025b), the dimensionality of the signals dramatically increases. Barrier et al. (2024)
demonstrated that a simple RNN is prone to catastrophic forgetting, where the RNN well
estimates the beginning of the signals but learns the end less effectively. To mitigate this issue,
they proposed to use bidirectional LSTM (bi-LSTM) blocks within the RNN architecture,
ensuring that both the beginning and the end of the signal are efficiently learned by the new
bi-LSTM blocks. Despite the good results in reconstructing $T_1$, $T_2$, and vascular parameters,
bi-LSTM still faces challenges in learning the middle part of the signals.

Generally, the dictionary is seen as a collection of simulated signals with no particular
spatial correlation. Another way to simulate signals is to acquire real parametric maps from
different subjects and then simulate the MRF images. This approach was first proposed by
Soyak et al. (2021) and later improved by Gu et al. (2024). Both studies utilize a UNet
(Ronneberger, Fischer and Brox, 2015) to infer $T_1$ and $T_2$ directly using the entire MRF
image as an input to the network, preserving valuable spatial information. Given the high-
dimensionality of the signals, the authors proposed adding attention layers (Vaswani et al.,
2017) to focus on the most important dimensions of the signals. More recently, Li and Hu
(2024) highlighted the limitations of using CNNs, which have a restricted receptive field
and capture spatial information only locally. To overcome these limitations, the authors pro-
posed using a Local-Global vision Transformer to capture spatial information globally as
well. However, capturing spatial information has a cost. Indeed, one needs to acquire, from a
large group of subjects, multiple $T_1$ and $T_2$ maps, which takes about 30 min for a complete
exam, making the data acquisition process much more costly than using more conventional
dictionaries. Additionally, since $T_1$ and $T_2$ need to be acquired at two different times, this
method introduces more errors from the registration of the acquired maps.

Finally, Boux et al. (2021) also proposed to use GLLiM (Deleforge, Forbes and Horaud,
2015), a model that casts MRF reconstruction into a Bayesian inverse problem and then
solves it using a learning approach. This method takes into account the high-dimensional
property of MRF signals by defining the low-dimensional variables as the regressors, which
in our case are the tissue parameters. By doing so, they start learning the *low-to-high* regres-
sion model from which they can derive the forward model parameters and then the *high-to-
low* regression model, from MRF signals to tissue parameters, as desired.

**3. Related work and positioning.** In practice, acquired MRF acquisitions come as a 4D matrix, made of a time series of 3D MRI images where each voxel contains the acquired, potentially long, fingerprint signal. In this work, we focus on designing efficient representations of large highly precise grids of simulated signals counterparts.

3.1. *Parsimonious representations of dictionaries.* The curse of dimensionality goes with what is often called the bless of dimensionality, which refers to the fact that in high-dimensional data sets, useful information actually lives in much smaller-dimensional parts of
the data space. An approach to the MRF reconstruction problem is then to reduce the dimension of the dictionary beforehand to reduce the matching or learning cost. The first to propose such an approach for efficient matching were McGivney et al. (2014), who applied Singular Value Decomposition (SVD) to the dictionary of signals. Once the decomposition is learned, one can project any new acquisition into the SVD low-dimensional subspace. However, when the size of the dictionary increases, computing the SVD becomes costly as it requires loading the complete dictionary into fast-access memory (e.g. the RAM). To address this, Yang et al. (2018) proposed using randomized SVD (Halko, Martinsson and Tropp, 2011). Golbabaee et al. (2019) also suggested applying SVD to the dictionary before training a neural network. Other methods proposed a non-Euclidean analysis of the dictionary space projecting the signals into a lower-dimensional manifold (Li and Hu, 2023). Reducing the dimension of high-dimensional dictionaries assumes that most of the information in the signals can be captured and represented in a much lower-dimensional subspace. Classical techniques include principal component analysis (PCA, Jolliffe and Cadima (2016)), probabilistic principal component analysis (PPCA, Tipping and Bishop (1999a)), factor analyzers (FA) (Freitas Lopes and West, 2004), sparse models (Zou, Hastie and Tibshirani, 2006; d'Aspremont et al., 2007; Archambeau and Bach, 2008), and newer methods such as diffusion maps (Coifman and Lafon, 2006). More flexible approaches are based on mixtures of the previous ones, such as mixtures of factor analyzers (MFA, McLachlan, Peel and Bean (2003)) introduced by Ghahramani and Hinton (1997); Ghahramani and Beal (1999) and extended by McLachlan, Peel and Bean (2003); Baek and McLachlan (2011), and mixtures of PPCA (MPPCA, Tipping and Bishop (1999b); Xu, Balzano and Fessler (2023)) with recent generalizations (Hong et al., 2023; Xu, Balzano and Fessler, 2023). Another mixture approach is called HDDC in Bouveyron, Girard and Schmid (2007) and HD-GMM in Bouveyron and Brunet-Saumard (2014) for High Dimensional Gaussian Mixture Models, which encompass many forms of MFA and MPPCA and generalize them. In particular, HD-GMM can be used to obtain multiple low-dimensional subspaces of different dimensions. For a review on high-dimensional clustering via mixtures, see Bouveyron and Brunet-Saumard (2014). Such a divide-and-conquer strategy has been used in MRF by first creating a global clustering of the signals and then reducing the dimension. This method was applied by Cauley et al. (2015), who performed K-way partitioning of the dictionary before using multiple cluster-associated PCA to reduce the dimension in each part. When matching a new signal, one first determines which cluster it belongs to and then applies the projection learned by the associated PCA model. More recently, Ullah et al. (2023) proposed a simpler approach: applying a clustering algorithm to the dictionary, before utilizing GPUs for the dictionary matching, enabling fast matching without the need for dimensionality reduction. However, this method still faces issues when the number of parameters to estimate is larger than 3.

Most of these methods are designed for batch data and are thus sensitive to hardware
limits such as memory, restricting the amount of data they can process. For instance, some
dictionaries exceed terabyte (TB) sizes. A simple solution is to down-sample data sets before
processing, potentially losing useful information. Another approach is to design incremental,
also referred to as online, variants that handle data sequentially in smaller groups. Although
such approaches can be used for genuine streaming data that is not stored after being seen, in
this work we consider cases where past samples can possibly be revisited. A number of incre-
mental approaches exist for dimension reduction techniques, see the recent SHASTA-PCA
(Gilman et al., 2023) and references therein, or Balzano, Chi and Lu (2018) for a review. To
our knowledge, far fewer solutions exist for mixtures. Estimation of such models is generally
based on maximum likelihood estimation via the Expectation-Maximization (EM) algorithm
(McLachlan and Krishnan (2008)). A preliminary attempt for an incremental MPPCA can be
found in Bellas et al. (2013) but it is based on heuristic approximations of the EM steps.
### 3.2. Contribution.
Herein, we propose to explore this divide-and-conquer strategy by introducing the framework of High Dimensional Mixtures of Elliptical Distributions (HD-MED) to overcome the challenge of estimating more than 2 parameters, e.g. $T_1, T_2, \delta f, B_1, \text{CBV}, R$. Considering the reconstruction of 6 parameters makes the associated dictionaries both large in size (an order of 13 terabytes of signals to represent the tissue heterogeneity) and dimension, framing our proposal around two key components. Building on HD-MED, a generalization of HD-GMM, we show how they can be used to simultaneously compress and cluster large-scale high-dimensional MRF dictionaries, and more generally any dataset. We derive a new online algorithm, based on a principled EM framework, to learn such a model from very large data volumes. We demonstrate the effectiveness of our approach on MRF reconstruction, showing results comparable to the high-dimensional dictionary matching referred to as full-matching in the next sections. This approach allows us to exceed the resolution, size used in current implementations and to reconstruct a larger number of MR parameter maps with improved accuracy, thereby advancing the clinical feasibility of MRF.

### 4. High Dimensional Mixtures of Elliptical Distributions.
Due to time constraints and computational complexity, data may be acquired with errors and at high noise levels. In such cases, using Gaussian models, which are known to be sensitive to outliers, is not recommended. Elliptical distributions represent a family of distributions that contains Gaussian distributions but also heavy-tailed distributions, such as the Student or Laplace distributions. Elliptical distributions allow more flexible modeling, with additional parameters and better fit to data than Gaussian distributions, see e.g. (Lange, Little and Taylor, 1989; Lange and Sin-shheimer, 1993; Neri, Depalle and Badeau, 2021; Kotz, Kozubowski and Podgorski, 2001). This additional flexibility is often referred to as robustness because the heavier tails make estimation less sensitive to outlying samples. In this paper, we consider a sub-class of elliptical distributions, which can be expressed as infinite mixtures of Gaussian distributions.

#### 4.1. Gaussian scale mixture (GSM) distributions.
Scale mixtures of multivariate Gaussian distributions, also referred to as normal/independent distributions (Lange and Sin-shheimer, 1993), are an important subclass of elliptical distributions whose definition is recalled below. Scale mixtures of Gaussians share good properties with Gaussian distributions. They are tractable, lead to tractable inference procedures and provide more robust results, in contrast to Gaussian distributions that usually suffer from sensitivity to outliers.

**DEFINITION 4.1 (Elliptical Distributions (ED)).** A continuous random vector $\mathbf{Y} \in \mathbb{R}^M$ follows a multivariate elliptical symmetric distribution if its probability density function (pdf) $p(\mathbf{y})$ is of the following form (see Cambanis, Huang and Simons (1981) or Kelker (1970)),

$$ (2) \qquad p(\mathbf{y}) = C_{p,g} |\mathbf{\Sigma}|^{-1/2} g \left( (\mathbf{y} - \boldsymbol{\mu})^{\top} \mathbf{\Sigma}^{-1} (\mathbf{y} - \boldsymbol{\mu}) \right), $$

where $\mathbf{\Sigma} \in \mathbb{R}^{M \times M}$ is the scale matrix with determinant $|\mathbf{\Sigma}|$, $\boldsymbol{\mu} \in \mathbb{R}^M$ is the location or mean vector, $C_{p,g}$ is a normalizing constant such that the pdf $p(\mathbf{y})$ integrates to one. The non-negative function $g$ is called the density generator and determines the shape of the pdf. When $\mathbf{Y}$ has density (2), we write $\mathbf{Y} \sim \mathcal{E}_M(\boldsymbol{\mu}, \mathbf{\Sigma}, g)$.

Note that the scale matrix $\mathbf{\Sigma}$ is not necessarily the covariance matrix, $\mathbf{\Sigma}$ is proportional to the covariance matrix if the latter exists. The pdf of the multivariate normal distribution is a special case of ED with $g(u) = \exp(-u^2/2)$. Another member of the elliptical family is the multivariate Student distribution. This distribution is well-studied in the literature (Kotz and Nadarajah, 2004) and admits a useful representation as a Gaussian scale mixture. Denoting $\mathcal{N}_M(\mathbf{y}; \boldsymbol{\mu}, \mathbf{\Sigma})$ a $M$-variate Gaussian distribution with mean $\boldsymbol{\mu}$ and covariance matrix $\mathbf{\Sigma}$, a Gaussian scale mixture distribution is a distribution of the following form.
**DEFINITION 4.2 (Gaussian scale mixture distributions (GSM)).** If $\mu$ is an $M$-dimensional vector, $\Sigma$ is an $M \times M$ positive definite symmetric matrix and $f$ is a pdf of a univariate positive variable $W \in \mathbb{R}^+$, then the $M$-dimensional density given by

$$ (3) \qquad p(\mathbf{y}) = \int_{\mathbb{R}^+} \mathcal{N}_M \left( \mathbf{y}; \mu, \frac{\Sigma}{w} \right) f(w) dw $$

is said to be an infinite mixture of scaled Gaussians or Gaussian scale mixture (GSM) with mixing distribution function $f$. If vector $\mathbf{Y}$ has density (3), we still write $\mathbf{Y} \sim \mathcal{E}_M(\mu, \Sigma, f)$ and refer to $W$ as the mixing variable.

In practice, we will consider mixing distribution $f_\theta$ that depends on some parameter $\theta$ and also write $\mathbf{Y} \sim \mathcal{E}_M(\mu, \Sigma, \theta)$. As already mentioned, famous GSM distributions include the multivariate Student distribution (when $f_\theta$ is the pdf of a chi-squared variable), the Pearson type VII distribution (when $f_\theta$ is a gamma distribution), the generalized Gaussian (when a power of $W$ follows the gamma distribution) and the Laplace distribution (when $W^{-1}$ follows an exponential distribution). It is straightforward to see that GSM are elliptical distributions. However, not all elliptical distributions can be reduced to scale mixtures. For the previous reason, characterization and a way to represent elliptical distributions as GSM are very valuable. In Gómez-Sánchez-Manzano, Gómez-Villegas and Marín (2006), conditions are given under which elliptical distributions are GSMs. The issue of finding the corresponding mixing distribution is also addressed. An illustration of these results for generalized Gaussian distributions is given by Gomez, Gomez-Villegas and Marin (2008).

## 4.2. Mixtures of High Dimensional GSM.
We consider finite mixtures of $K$ GSM distributions, with a specific parameterization to handle potentially high-dimensional observations.

**DEFINITION 4.3 (Finite mixture).** $\mathbf{Y}$ follows a finite mixture model if its pdf writes as

$$ (4) \qquad p(\mathbf{y}) = \sum_{k=1}^{K} \pi_k f_k(\mathbf{y}), $$

with $\pi_k \in [0, 1]$ the mixing weights summing to one and $f_k$ the pdf of the $k$th component.

A mixture model where each $f_k$ is a GSM (3), $\mathcal{E}(\mu_k, \Sigma_k, \theta_k)$, is referred to as a mixture of ED (MED). The number of parameters in MED grows quadratically with the dimension $M$ due to the scale matrices $\Sigma_k$. For large $M$, this can be problematic for estimation. In the Gaussian case, to reduce the number of parameters, Bouveyron, Girard and Schmid (2007) proposed a family of parsimonious Gaussian mixtures, using the eigenvalue decomposition of the covariance matrices. We extend this idea to MED by reparameterizing $\Sigma_k$ as follows,

$$ (5) \qquad \Sigma_k = D_k A_k D_k^T, $$

where $D_k$ is an $M \times M$ orthogonal matrix, which contains the eigenvectors of $\Sigma_k$ and $A_k$ is an $M \times M$ diagonal matrix, which contains the associated eigenvalues in decreasing order. The key idea of Bouveyron, Girard and Schmid (2007) is to consider that each cluster lies in a low-dimensional subspace of dimension $d_k < M$, which can be expressed by assuming that

$$ (6) \qquad A_k = \mathrm{diag}(a_{k1}, \dots, a_{kd_k}, b_k, \dots, b_k), $$

where $a_{k1}, \dots, a_{kd_k}$ are the $d_k$ largest eigenvalues of $\Sigma_k$ and $b_k$ is a small negligible value. The $d_k$ eigenvectors associated to the first $d_k$ eigenvalues $\{a_{k1}, \dots, a_{kd_k}\}$ define a cluster-specific subspace $\mathbb{E}_k$, which captures the main cluster shape. The orthogonal subspace is
denoted by $\mathbb{E}_k^\perp$. Let $\tilde{D}_k$ consists of the $d_k$ first columns of $D_k$ supplemented by $(M - d_k)$ zero columns and $\overline{D}_k = (D_k - \tilde{D}_k)$. It follows that $P_k(y) = \tilde{D}_k \tilde{D}_k^T(y - \mu_k) + \mu_k$ and $P_k^\perp(y) = \overline{D}_k \overline{D}_k^T(y - \mu_k) + \mu_k$ are the projections of $y$ on $\mathbb{E}_k$ and $\mathbb{E}_k^\perp$ respectively. This parameterization allows us to handle high-dimensions in an efficient way. For instance, the quadratic form of the Mahalanobis distance, appearing in the generator $g$ in (2), writes

$$ (7) \qquad \begin{aligned} (y - \mu_k) D_k A_k^{-1} D_k^T (y - \mu_k) &= (y - \mu_k)^T \tilde{D}_k A_k^{-1} \tilde{D}_k^T (y - \mu_k) \\ &\quad + (y - \mu_k)^T \overline{D}_k A_k^{-1} \overline{D}_k^T (y - \mu_k) \\ &= \|\mu_k - P_k(y)\|_{\tilde{\Sigma}_k^{-1}}^2 + \frac{1}{b_k} \|y - P_k(y)\|^2, \end{aligned} $$

where $\|\cdot\|_{\tilde{\Sigma}_k^{-1}}^2$ is the norm defined by $\|y\|_{\tilde{\Sigma}_k^{-1}}^2 = y^T \tilde{\Sigma}_k^{-1} y$ with $\tilde{\Sigma}_k^{-1} = \tilde{D}_k A_k^{-1} \tilde{D}_k^T$. Equation (7) uses the definitions of $P_k$ and $P_k^\perp$ and $\|\mu_k - P_k^\perp(y)\|^2 = \|y - P_k(y)\|^2$. The gain comes from the fact that (7) does not depend on $P_k^\perp$ and thus does not require the computation of the $(M - d_k)$ latest columns of $D_k$, the eigenvectors associated to the smallest eigenvalues. Similarly, determinants can be efficiently computed as $\log(|\Sigma_k|) = (\sum_{m=1}^{d_k} \log(a_{km})) + (M - d)\log(b_k)$. This parameterization that now only depends on $\tilde{D}_k$ and not on the complete $D_k$, is indicated by denoting the corresponding ED as $\mathcal{H}\mathcal{E}_{Md_k}(\mu_k, \tilde{D}_k^*, a_k, b_k, \theta_k)$, where $\tilde{D}_k^*$ is matrix $\tilde{D}_k$ with the last zero $M - d_k$ columns omitted. We then refer to a MED coupled with this parameterization as a high-dimensional MED (HD-MED).

**DEFINITION 4.4 (HD-MED).** A random vector $Y \in \mathbb{R}^M$ follows a HD-MED distribution if for all $k \in [1:K]$, the pdf of the $k^{th}$ component $f_k$ is an ED $\mathcal{H}\mathcal{E}_{Md_k}(\mu_k, \tilde{D}_k^*, a_k, b_k, \theta_k)$ with reparameterization given by (5) and (6). We denote

$$ Y \sim \mathcal{M}\mathcal{H}\mathcal{E}_{Md}((\pi_k, \mu_k, \tilde{D}_k^*, a_k, b_k, \theta_k)_{k=1}^K), \text{ with } d = (d_1, \dots, d_k), a_k = (a_{k1}, \dots, a_{kd_k}). $$

## 5. Dimension reduction with HD-MED.

5.1. Latent variable dimension reduction. Standard PCA is defined without referring to a probabilistic model. Given a set of observations in $\mathbb{R}^M$, their $M \times M$ empirical covariance matrix is decomposed into eigenvalues and eigenvectors. For any $y \in \mathbb{R}^M$, a lower dimensional ($d << M$) representation can then be obtained by considering its projection to a lower dimensional subspace $\hat{y} = \hat{\Sigma}_d^T y$ where $\hat{\Sigma}_d$ is the matrix containing the $d$ first eigenvectors of the empirical covariance matrix. Its reconstruction in $\mathbb{R}^M$, optimal in the sense of the squared reconstruction error, can be obtained with $\tilde{y} = \hat{\Sigma}_d \hat{y}$. Alternatively, if $y$ is assumed to be a realization of a random vector $Y \sim \mathcal{H}\mathcal{E}_{Md}(\mu, \tilde{D}^*, a, b, \theta)$, a low dimensional representation of $y$ can be justified using the following latent variable model representation of $Y$. The following proposition generalizes the developments of Archambeau, Delannay and Verleysen (2006); Guo and Bondell (2023) for the Student distribution.

**PROPOSITION 5.1 (HD-ED latent variable model).** Let $d \le M - 1$, $Y \in \mathbb{R}^M$, $X \in \mathbb{R}^d$, $E \in \mathbb{R}^M$, $W \in \mathbb{R}^+$ be random variables, $V \in \mathbb{R}^{M \times d}$ a matrix of linearly independent columns, $\mu \in \mathbb{R}^M$ a vector and $f_\theta$ the pdf of a positive univariate random variable defined by some parameter $\theta$. Assume that $Y = VX + \mu + E$, $(X|W = w) \sim \mathcal{N}(0_d, w^{-1}I_d)$, $(E|W = w) \sim \mathcal{N}(0_M, b w^{-1}I_M)$ and $W \sim f_\theta$, then, $Y \sim \mathcal{H}\mathcal{E}_{Md}(\mu, \tilde{D}^*, a, b, \theta)$, with $DAD^T =$
$b\mathbf{I}_M + \mathbf{V}\mathbf{V}^T$ the eigenvalue decomposition of $b\mathbf{I}_M + \mathbf{V}\mathbf{V}^T$, $\mathbf{A} = \text{diag}(a_1, \dots, a_d, b, \dots, b)$ the ordered eigenvalues and $\tilde{\mathbf{D}}^*$ the matrix containing the first $d$ eigenvectors as columns.

Additionally, denoting by $\mathbf{U} = b\mathbf{I}_d + \mathbf{V}^T\mathbf{V}$, we have,

$$ (8) \qquad (\mathbf{X}|\mathbf{Y} = \mathbf{y}, W = w) \sim \mathcal{N}(\mathbf{U}^{-1}\mathbf{V}^T(\mathbf{y} - \boldsymbol{\mu}), w^{-1}b\mathbf{U}^{-1}). $$

It follows that $\mathbb{E}[\mathbf{Y}|\mathbf{X} = \mathbf{x}] = \mathbf{V}\mathbf{x} + \boldsymbol{\mu}$ and $\mathbb{E}[\mathbf{X}|\mathbf{Y} = \mathbf{y}] = \mathbf{U}^{-1}\mathbf{V}^T(\mathbf{y} - \boldsymbol{\mu})$.

PROOF. It comes from the first three assumptions that

$$ (\mathbf{Y}|W = w) \sim \mathcal{N}(\boldsymbol{\mu}, w^{-1}(\mathbf{V}\mathbf{V}^T + b\mathbf{I}_M)) $$

and from definition 4.2 that $\mathbf{Y} \sim \mathcal{E}_M(\boldsymbol{\mu}, \mathbf{V}\mathbf{V}^T + b\mathbf{I}_M, f_\theta)$. Since $\mathbf{V}\mathbf{V}^T$ is of rank $d$, $\mathbf{V}\mathbf{V}^T + b\mathbf{I}_M$ admits an eigenvalue decomposition $\mathbf{D}\mathbf{A}\mathbf{D}^T$ as stated. Distribution (8) follows from standard Gaussian vectors properties. Using the tower property, $\mathbb{E}[\mathbf{Y}|\mathbf{X} = \mathbf{x}] = \mathbb{E}[\mathbb{E}[\mathbf{Y}|\mathbf{X} = \mathbf{x}, W]] = \mathbf{V}\mathbf{x} + \boldsymbol{\mu}$ and $\mathbb{E}[\mathbf{X}|\mathbf{Y} = \mathbf{y}] = \mathbb{E}[\mathbb{E}[\mathbf{X}|\mathbf{Y} = \mathbf{y}, W]] = \mathbf{U}^{-1}\mathbf{V}^T(\mathbf{y} - \boldsymbol{\mu})$. $\square$

The previous proposition states that $\mathbf{y} \in \mathbb{R}^M$ from an HD-ED, can also be seen as originating from a model with a lower-dimensional latent variable $\mathbf{X} \in \mathbb{R}^d$. Hence, a natural alternative to the standard PCA projection, is the conditional mean $\hat{\mathbf{y}} = \mathbb{E}[\mathbf{X}|\mathbf{Y} = \mathbf{y}]$, that is

$$ (9) \qquad \hat{\mathbf{y}} = Q(\mathbf{y}) = \mathbf{U}^{-1}\mathbf{V}^T(\mathbf{y} - \boldsymbol{\mu}), $$

and as a reconstruction or an approximation of the original information, the conditional mean $\hat{\mathbf{y}} = \mathbb{E}[\mathbf{Y}|\mathbf{X} = \hat{\mathbf{y}}]$, that is $\hat{\mathbf{y}} = \mathbf{V}\hat{\mathbf{y}} + \boldsymbol{\mu}$. Using the previous formulas, lower dimensional representations can thus be obtained using $\mathbf{U}$ and $\mathbf{V}$ but when estimating the parameters of the HD-ED, we get estimates for $(a_1, \dots, a_d)$, $b$ and $\tilde{\mathbf{D}}^*$. However, using that $\mathbf{V}$ is of rank $d$ and $\mathbf{D}\mathbf{A}\mathbf{D}^T = b\mathbf{I}_M + \mathbf{V}\mathbf{V}^T$, we can set

$$ (10) \qquad \mathbf{V} = \tilde{\mathbf{D}}^* \sqrt{\text{diag}(a_1, \dots, a_d) - b\mathbf{I}_d}, $$

and deduce $\mathbf{U}$ straightforwardly. Interestingly, the projection and reconstruction formulas do not directly depend on $f_\theta$. The dependence is only implicit via the impact of $f_\theta$ on the estimation of the other parameters.

5.2. Cluster globally, reduce locally. The proposed dimension reduction and reconstruction, using the HD-ED model, generalizes to GSM, PPCA (Tipping and Bishop, 1999c) and robust PPCA (Archambeau, Delannay and Verleysen, 2006), the former using Gaussian and the latter Student distributions. Both have been extended to account for potential heterogeneity in data, considering mixtures, with the MPPCA model (Tipping and Bishop, 1999d) and its Student-based robust version (Archambeau, Delannay and Verleysen, 2008). In these models, all clusters are assumed to live in subspaces of the same dimension $d$. In the Gaussian case, an extension, allowing varying reduced dimensions $d_k$ across clusters, have been proposed by Bouveyron, Girard and Schmid (2007) with their high-dimensional Gaussian mixture model (HD-GMM). The proposed HD-GMM parameterization allows to handle high-dimensional data in a computationally efficient way. However, it does not provide an actual lower dimensional representation of the data. While such a reduced-dimensional representation may often not be needed, it may be crucial to deal with hardware or software limitations. Originally, HD-GMM have not been designed for dimension reduction or compression but rather for clustering and density estimation in high-dimensional heterogeneous settings. Our HD-MED model generalizes the covariance matrix decomposition to the scale matrix of a GSM distribution. In addition, we describe how it can be further exploited as a dimension reduction technique. As mixture models, HD-MED can be used for clustering data into $K$ clusters. For any possible observation $\mathbf{y}$, a HD-MED model provides a probability $r_k(\mathbf{y})$ that
y is assigned to cluster $k$ for each $k = 1:K$. Using (9), a reduced-dimension representation $\hat{y}_k$ of y, for each of the K different subspaces, is denoted by $\hat{y}_k = Q_k(y)$ and given by,

$$ (11) \qquad \hat{y}_k = Q_k(y) = U_k^{-1} V_k^T (y - \mu_k), $$

while its reconstruction $\tilde{y}_k$ in the original space is given by $\tilde{y}_k = V_k \hat{y}_k + \mu_k$. In practice, it is reasonable to use as a reduced-dimension representation of y only the one corresponding to the most probable cluster $k$, i.e. with the highest $r_k(y)$. In this setting, HD-MED acts as a divide-and-conquer paradigm by initially clustering the data into K clusters and then performing cluster-specific data reduction. The divide step allows a much more effective reduction than if a single subspace was considered, while in the conquer step, little information is lost, as for any new observation y, cluster assignment probabilities $r_k(y)$ can be straightforwardly computed to decide on the best reduced representation to be used. However, it is important to keep track of clustering information for each observation. The reduced representations cannot be pooled back altogether, as they are likely to become impossible to distinguish across clusters. Also, as summarized in Algorithm 1 and illustrated in Section 7, each reduced cluster may have to be processed separately but this additional cost is negligible compared to the hardware and software gain of a more efficient representation.

The Expectation-Maximization (EM) algorithm (Dempster, Laird and Rubin (1977)) that iteratively maximizes the conditional expectation of the complete-data log-likelihood, is commonly used to estimate finite mixture models. Gradient algorithms can also be used but have not been yet very popular, see a more detailed discussion in Appendix Section 2. The number of clusters K, and their respective inner dimension $d_k$ are hyperparameters that need to be tuned prior to the EM steps. As a simple solution, we use the Bayesian Information Criterion (BIC) to tune K and the $d_k$'s at the same time. In contrast to MPPCA solutions which assume the same subspace dimension d for all clusters, the possibility to handle different $d_k$'s and to allow non Gaussian cluster shapes is important for the target applications involving datasets that are very large. Using different dimensions across clusters is likely to yield a more efficient reduced representation of the data as illustrated in Section 7.

The standard, or batch EM algorithm needs all the dataset to be loaded in a fast access memory (e.g the RAM) which is often limited as this kind of memory is expensive. In the case of large-scale dataset, the RAM is often overloaded and supported by a slow-memory which makes EM iterations very slow. Batch sizes are then limited by resource constraints, so that very large data sets need either to be downsampled or to be handled in an incremental manner. Incremental versions of EM exist and can be adapted to our setting. In section 6, we provide a way to deal with this large-scale case by using an online version of EM.

# 6. Online Learning of High Dimensional Mixtures of Elliptical Distributions.

## 6.1. Online EM, main assumptions.
When the data volume is too large, EM becomes slow because of multiple data transfer between the RAM and the computer storage. A way to handle large volumes is to use online learning. Online learning refers to procedures dealing with data acquired sequentially. Online variants of EM, among others, are described in Cappé and Moulines (2009); Maire, Moulines and Lefebvre (2017); Karimi et al. (2019a,b); Fort, Moulines and Wai (2020); Kuhn, Matias and Rebafka (2020); Nguyen, Forbes and McLachlan (2020); Dieuleveut et al. (2023). As an archetype of such algorithms, we consider the online EM of Cappé and Moulines (2009) which belongs to the family of stochastic approximation algorithms (Borkar (2009); Dieuleveut et al. (2023)). This algorithm has been theoretically studied and extended. However, it is designed only for distributions that admit a data augmentation scheme yielding a complete likelihood of the exponential family form, see (12) below. This case is already very broad, including e.g. Gaussian, gamma, Student
distributions and mixtures of those. We recall the main assumptions required and the online EM iteration, based on a latent variable formulation.

Assume $(\mathbf{Y}_i)_{i=1}^N$ is a sequence of $N$ i.i.d. replicates of a random variable $\mathbf{Y} \in \mathcal{Y} \subset \mathbb{R}^M$, observed one at a time. Extension to successive mini-batches is straightforward (Nguyen, Forbes and McLachlan (2020)). In addition, $\mathbf{Y}$ is assumed to be the visible part of $(\mathbf{Y}, \mathbf{Z})$, where $\mathbf{Z} \in \mathbb{R}^l$ is a latent variable, e.g. the unknown component label in a mixture or a mixing weight in a GSM formulation, and $l \in \mathbb{Z}_+$. For $i \in [1:N]$ then, each $\mathbf{Y}_i$ is the visible part of $(\mathbf{Y}_i, \mathbf{Z}_i)$. Suppose $\mathbf{Y}$ arises from some data generating process (DGP) characterised by a probability density function $f(\mathbf{y}; \mathbf{\Theta}_0)$, with unknown parameters $\mathbf{\Theta}_0 \in \mathbb{T} \subseteq \mathbb{R}^p$, for $p \in \mathbb{Z}_+$.

Using the sequence $(\mathbf{Y}_i)_{i=1}^N$, the method of Cappé and Moulines (2009) allows to sequentially estimate $\mathbf{\Theta}_0$ provided the following assumptions are met:

(A1) The complete-data likelihood for $(\mathbf{Y}, \mathbf{Z})$ is of the exponential family form:

$$ (12) \qquad f_c(\mathbf{y}, \mathbf{z}; \mathbf{\Theta}) = h(\mathbf{y}, \mathbf{z}) \exp \{[\mathbf{s}(\mathbf{y}, \mathbf{z})]^{\top} \phi(\mathbf{\Theta}) - \psi(\mathbf{\Theta})\}, $$

with $h: \mathbb{R}^{M+l} \to [0, \infty)$, $\psi: \mathbb{R}^p \to \mathbb{R}$, $\mathbf{s}: \mathbb{R}^{M+l} \to \mathbb{R}^q$, $\phi: \mathbb{R}^p \to \mathbb{R}^q$, for $q \in \mathbb{Z}_+$.

(A2) The function

$$ (13) \qquad \bar{s}(\mathbf{y}; \mathbf{\Theta}) = \mathbb{E}[\mathbf{s}(\mathbf{Y}, \mathbf{Z}) | \mathbf{Y} = \mathbf{y}; \mathbf{\Theta}] $$

is well-defined for all $\mathbf{y}$ and $\mathbf{\Theta} \in \mathbb{T}$, where $\mathbb{E}[|\mathbf{Y} = \mathbf{y}; \mathbf{\Theta}|]$ is the conditional expectation when $\mathbf{X}$ arises from the DGP characterised by $\mathbf{\Theta}$.

(A3) There is a convex $\mathbb{S} \subseteq \mathbb{R}^q$, satisfying: (i) for all $\gamma \in (0, 1)$, $\mathbf{s} \in \mathbb{S}$, $\mathbf{y} \in \mathcal{Y}$, and $\mathbf{\Theta} \in \mathbb{T}$, $(1-\gamma)\mathbf{s} + \gamma\bar{s}(\mathbf{y}; \mathbf{\Theta}) \in \mathbb{S}$; and (ii) for any $\mathbf{s} \in \mathbb{S}$, the function $Q(\mathbf{s}; \mathbf{\Theta}) = \mathbf{s}^{\top}\phi(\mathbf{\Theta}) - \psi(\mathbf{\Theta})$ has a unique global maximizer on $\mathbb{T}$ denoted by

$$ (14) \qquad \mathbf{\Theta}(\mathbf{s}) = \underset{\mathbf{\Theta} \in \mathbb{T}}{\arg \max} Q(\mathbf{s}; \mathbf{\Theta}). $$

Let $(\gamma_i)_{i=1}^N$ be a sequence of learning rates in $(0, 1)$ and $\mathbf{\Theta}^{(0)} \in \mathbb{T}$ an initial estimate of $\mathbf{\Theta}_0$. For each $i \in [1:N]$, the online EM of Cappé and Moulines (2009) proceeds by computing

$$ (15) \qquad \mathbf{s}^{(i)} = \gamma_i \bar{s}(\mathbf{y}_i; \mathbf{\Theta}^{(i-1)}) + (1-\gamma_i) \mathbf{s}^{(i-1)}, $$

$$ (16) \qquad \text{and} \qquad \mathbf{\Theta}^{(i)} = \bar{\mathbf{\Theta}}(\mathbf{s}^{(i)}), $$

where $\mathbf{s}^{(0)} = \bar{s}(\mathbf{y}_1; \mathbf{\Theta}^{(0)})$. Theorem 1 of Cappé and Moulines (2009) shows that the sequence $(\mathbf{\Theta}^{(i)})_{i=1:N}$ of estimators of $\mathbf{\Theta}_0$ satisfies a convergence result to stationary points of the likelihood (cf. Cappé and Moulines (2009) for a more precise statement).

## 6.2. Online EM for HD-ED.
We now derive the online EM for HD-ED. The extension to HD-MED mixtures is straightforward and is detailed in Nguyen and Forbes (2022) or in the supplementary Section 3. The weight distribution $f$ in the GSM formulation (3) is assumed to belong to the exponential family. This case may seem restrictive but it encompasses a number of ED as the Gaussian, Student, Normal Inverse Gamma with no skewness etc. distributions.

**PROPOSITION 6.1 (HD-ED exponential form).** Let $\mathbf{Y}$ be a HD-ED distributed variable, $\mathbf{Y} \sim \mathcal{H}\mathcal{E}_{\mathrm{MD}}(\boldsymbol{\mu}, \tilde{\mathbf{D}}^*, \mathbf{a}, \mathbf{b}, \boldsymbol{\theta})$, and $\mathbf{W}$ a weight variable with pdf $f_{\boldsymbol{\theta}}$. The set of parameters is denoted by $\mathbf{\Theta} = (\boldsymbol{\mu}, \tilde{\mathbf{D}}^*, \mathbf{a}, \mathbf{b}, \boldsymbol{\theta})$, with $\tilde{\mathbf{D}}^*$ defined by its column vectors $\tilde{\mathbf{D}}^* = [\mathbf{d}_1, \dots, \mathbf{d}_d]$. If $\mathbf{W}$ belongs to the exponential family, i.e. $f_{\boldsymbol{\theta}}(\mathbf{w}) = h_w(\mathbf{w}) \exp [\mathbf{s}_w(\mathbf{w})^T \phi_w(\boldsymbol{\theta}) - \psi_w(\boldsymbol{\theta})]$, the complete data likelihood

$$ (17) \qquad f_c(\mathbf{y}, \mathbf{w}; \mathbf{\Theta}) = f_{\boldsymbol{\theta}}(\mathbf{w}) \mathcal{N}_M(\mathbf{y}; \boldsymbol{\mu}, \mathbf{w}^{-1} \mathbf{D} \mathbf{A} \mathbf{D}^T), $$
can be expressed in an exponential family form (12) with

$$
(18) \qquad \mathbf{s}(\mathbf{y}, w) = \begin{bmatrix} w\mathbf{y} \\ w \operatorname{vec}(\mathbf{y}\mathbf{y}^T) \\ w\mathbf{y}^T\mathbf{y} \\ w \\ \mathbf{s}_w(w) \end{bmatrix}, \quad \phi(\boldsymbol{\Theta}) = \begin{bmatrix} \sum_{m=1}^{d} \left(\frac{1}{a_m} - \frac{1}{b}\right) \mathbf{d}_m \mathbf{d}_m^T \boldsymbol{\mu} + \frac{1}{b} \boldsymbol{\mu} \\ \frac{1}{2} \sum_{m=1}^{d} \left(\frac{1}{b} - \frac{1}{a_m}\right) \operatorname{vec}(\mathbf{d}_m \mathbf{d}_m^T) \\ -\frac{1}{2b} \\ \frac{1}{2} \sum_{m=1}^{d} \left(\frac{1}{b} - \frac{1}{a_m}\right) \boldsymbol{\mu}^T \mathbf{d}_m \mathbf{d}_m^T \boldsymbol{\mu} - \frac{1}{2b} \boldsymbol{\mu}^T \boldsymbol{\mu} \\ \phi_w(\boldsymbol{\theta}) \end{bmatrix},
$$
$$
\psi(\boldsymbol{\Theta}) = \frac{1}{2} \sum_{m=1}^{d} \log a_m + \frac{M-d}{2} \log b + \psi_w(\boldsymbol{\theta}),
$$

where vec denotes the vectorization operator (Schott (2016)).

PROOF. The proof is detailed in supplementary Section 4.

<figure><img src="image_2.png" /></figure>

The online EM algorithm (OEM) consists, as the batch EM, of two steps, the computation of the sufficient statistics (13) and the maximization of the likelihood (14). For the first step, we need to compute $\bar{s}(\mathbf{y}; \mathbf{\Theta}) = \mathbb{E}[\mathbf{s}(\mathbf{Y}, \mathbf{Z})|\mathbf{Y}=\mathbf{y}; \mathbf{\Theta}]$. This quantity requires to compute the following expectations $\mathbb{E}[W|\mathbf{Y}=\mathbf{y}; \mathbf{\Theta}]$, and $\mathbb{E}[s_w(W)|\mathbf{Y}=\mathbf{y}; \mathbf{\Theta}]$.

PROPOSITION 6.2 (Expectations of sufficient statistics). Let $\mathbf{Y} \sim \mathcal{H}\mathcal{E}_{M_d}(\boldsymbol{\mu}, \tilde{\mathbf{D}}^*, \mathbf{a}, \mathbf{b}, \boldsymbol{\theta})$, $\mathbf{\Theta} = (\boldsymbol{\mu}, \tilde{\mathbf{D}}^*, \mathbf{a}, \mathbf{b}, \boldsymbol{\theta})$, and $W$ the mixing variable $W \sim f_{\boldsymbol{\theta}}$. Then $\mathbf{Y}$ has density (2) with generator $g$ and Mahalanobis distance $u = (\mathbf{y} - \boldsymbol{\mu})^T \mathbf{D} \mathbf{A}^{-1} \mathbf{D}^T (\mathbf{y} - \boldsymbol{\mu})$ defined as in (7). With $g'$ denoting the derivative of $g$, it follows that

$$
(19) \qquad \mathbb{E}[W | \mathbf{Y} = \mathbf{y}; \mathbf{\Theta}] = - \frac{2}{(2\pi)^{M/2}} \frac{g'(u)}{g(u)}.
$$

PROOF. Equation (19) results from equations (2) and (3), which lead to

$$
\mathbb{E}[W | \mathbf{Y} = \mathbf{y}; \mathbf{\Theta}] = \frac{1}{(2\pi)^{M/2} g(u)} \int_{\mathbb{R}^{+}} w f_{\theta}(w) \exp(-\frac{w}{2}u) dw = - \frac{2}{(2\pi)^{M/2}} \frac{g'(u)}{g(u)}.
$$

<figure><img src="image_3.png" /></figure>

In contrast, there is no general formula for the expectation of $s_w(W)$, which depends on
the mixing distribution $f_{\theta}$. Once the expectation of the sufficient statistics in (13) is com-
puted, we can update it following (15).

The next OEM step is the maximization step described in (14), which gives an estimation
of the parameters at each iteration. The solution for $\theta$ varies with $f_{\theta}$, but solutions for $\mu$, $\mathbf{A}$,
and $\tilde{\mathbf{D}}^*$ can be derived as follows. Let $\tilde{\mathbf{\Theta}}(\mathbf{s})$ be defined as the unique maximizer of function
$Q(\mathbf{s}, \mathbf{\Theta}) = \mathbf{s}^T \phi(\mathbf{\Theta}) - \psi(\mathbf{\Theta})$ with $\mathbf{s}$ a vector that matches the definition and dimension of
$\phi(\mathbf{\Theta})$ in (18), and can be conveniently written as

$$
(20) \qquad \mathbf{s} = \begin{bmatrix} s_1 \\ \operatorname{vec}(S_2) \\ s_3 \\ s_4 \\ s_5 \end{bmatrix},
$$
with $s_1$ an M-dimensional vector, $S_2$ an $M \times M$ matrix, and $s_3, s_4, s_5$ three scalar values. Parameters are updated by maximizing $Q$ with respect to $\Theta$. $\Theta(s)$ is defined as the root of the first-order condition involving $J_{\phi}(\Theta) = \frac{\partial \phi}{\partial \Theta}$ the Jacobian of $\phi$,

$$ (21) \qquad J_{\phi}(\Theta)s - \frac{\partial \psi}{\partial \Theta}(\Theta) = 0. $$

Computing gradients leads to $\tilde{\Theta}(s) = (\bar{\mu}(s), \tilde{D}^*(s), \bar{A}(s), \bar{\theta}(s))$, where $\bar{\mu}(s)$, and $\bar{A}(s)$ are closed form, and $\tilde{D}^*(s)$ can be found using Riemannian optimization.

**PROPOSITION 6.3 (Maxima).** We use an ECM-like procedure (Meng and Rubin (1993)) by optimizing parameters separately and incorporating them during the optimization of each other parameter. $(\mu, A)$ are easily optimized using (21) and computing the gradients gives

$$ (22) \qquad \bar{\mu} = \frac{s_1}{s_4}, $$

$$ (23) \qquad \bar{a}_m = d_m^T (S_2 + s_4 \bar{\mu}\bar{\mu}^T - 2\bar{\mu}s_1^T) d_m \quad \text{for } m \in [1: d], $$

$$ (24) \qquad \bar{b} = \frac{1}{M-d} \left( s_4 \bar{\mu}^T \bar{\mu} + s_3 - 2 \bar{\mu}^T s_1 - \sum_{m=1}^{d} \bar{a}_m \right). $$

Maximizing over $\tilde{D}^*$ has to take into account that $\tilde{D}^* \in St(M, d)$, the Stiefel manifold of the $M \times d$ matrices. Using the expressions of $\mu$ and $A$, omitting parts on $\theta$ we have,

$$ (25) \qquad \tilde{D}^* = \arg\max_{\tilde{D}^* \in St(M,d)} \sum_{m=1}^{d} \left( \frac{1}{\bar{a}_m} - \frac{1}{b} \right) d_m^T (2\bar{\mu}s_1^T - S_2 - s_4\bar{\mu}\bar{\mu}^T) d_m. $$

For $\bar{\theta}$, a general closed-form expression is not available, but if there are no particular constraints it results from solving the following equation

$$ (26) \qquad s_5 \frac{\partial \phi_w}{\partial \theta}(\bar{\theta}) - \frac{\partial \psi_w}{\partial \theta}(\bar{\theta}) = 0. $$

**PROOF.** We compute the gradients and use (21) for $\mu$ and $A$, however for the $\tilde{D}^*$ we only plug-in the optimized values $\bar{\mu}$, and $\bar{A}$ in (17) and solve

$$ (27) \qquad \arg\max_{\tilde{D}^* \in St(M,d)} s^T \phi (\bar{\mu}, \tilde{D}^*, \bar{A}, \bar{\theta}) - \psi (\bar{\mu}, \tilde{D}^*, \bar{A}, \bar{\theta}), $$

which does not depend on the value of $\bar{\theta}$. $\square$

In practice (25) is solved using a Riemannian optimization framework in the setting where $M \gg d$ (Wen and Yin, 2013).

**7. Application to magnetic resonance fingerprinting (MRF) reconstruction.** MRF provides multiple tissue parameter maps from shorter acquisition times, thanks to the simultaneous application of transient states excitation and highly undersampled k-space read-outs. These two aspects have a combined impact on acquisition times and image reconstruction accuracy. More undersampling allows more parameter estimations but is also responsible for larger undersampling errors, noise and artifacts, reducing map reconstruction accuracy.

In some earlier work (Oudoumanessah et al., 2025), the Gaussian version of our procedure, referred to as HD-GMM, was evaluated, for the reconstruction, in an ideal setting, of fully
**Algorithm 1** Divide & Conquer high-dimensional matching for MRF reconstruction

**Input** Dictionary of (signal, parameters) pairs $\mathcal{D}_f = \{\mathbf{y}_i, \mathbf{t}_i\}_{i=1:N}, N \gg 1$, $\mathbf{t}_i \in \mathbb{R}^L, \mathbf{y}_i \in \mathbb{R}^M, M \gg 1$. In vivo acquired signals $\{\tilde{\mathbf{y}}_j\}_{j=1:N}, \tilde{\mathbf{y}}_j \in \mathbb{R}^M$.

1. **Reduced dimension representation of the dictionary:** $\{\hat{\mathbf{y}}_i, \mathbf{t}_i, r(\mathbf{y}_i)\}_{i=1:N}$

1.1 **Online HD-MED inference from** $\{\mathbf{y}_i\}_{i=1:N}: K$ clusters, $d_k < M$ for $k=1:N$ $\implies$ cluster assignment probabilities and cluster-wise projections $(\mathbf{r}, \mathbf{Q}) = \{r_k(\cdot), Q_k(\cdot)\}_{k=1:N}$

1.2 **Cluster-wise fingerprint reductions:** $\{\mathbf{y}_i\}_{i=1:N}, \mathbf{r}, \mathbf{Q} \implies \{\hat{\mathbf{y}}_i = Q_k(\mathbf{y}_i), i \in I_k\}$ with $I_k = \{i, \text{s.t } k = \arg\max_{\ell} r_{\ell}(\mathbf{y}_i)\}, \text{ for } k=1:N$

2. **Cluster-wise matching of acquired signals:**

2.1 **Cluster-wise invivo signal reductions:** Use learned $(\mathbf{r}, \mathbf{Q})$ from step 1.1 to obtain $\{Q_k(\tilde{\mathbf{y}}_j), j \in \tilde{I}_k\}$ with $\tilde{I}_k = \{j, \text{s.t } k = \arg\max_{\ell} r_{\ell}(\tilde{\mathbf{y}}_j)\}, \text{ for } k=1:N$

2.2 **Matching:** For $k=1:N$, for $j \in \tilde{I}_k$, determine $i(\tilde{\mathbf{y}}_j) = \arg\min_{i \in I_k} d(Q_k(\tilde{\mathbf{y}}_j), \hat{\mathbf{y}}_i)$ and set $\tilde{\mathbf{t}}_j = \mathbf{t}_i(\tilde{\mathbf{y}}_j)$

**Return** Matched tissue properties $\{\tilde{\mathbf{t}}_j\}_{j=1:N}, \tilde{\mathbf{t}}_j \in \mathbb{R}^L$

sampled acquisitions, targeting $L=3$ parameters. In this scenario, it was reported that HD-GMM, coupled with the online EM algorithm, achieved results comparable to full dictionary matching while significantly reducing reconstruction times. However, as we report in this section, HD-GMM performance degrades when dealing with largely undersampled acquisitions and becomes inadequate for estimating six parameters. Performance can be maintained by considering an elliptical version of our procedure (HD-MED) less sensitive to outliers.

In this section, we thus show how we can accurately reconstruct six parameter maps from *in vivo* undersampled acquisitions (Section 7.1) by leveraging an extensive high-resolution dictionary (Section 7.2) and the HD-MED model. Algorithm 1 provides a schematic summary of our procedure. Figure 1 provides an illustration of our matching strategy or step 2 in Algorithm 1, once the HD-MED model has been estimated. Each acquired signal is first assigned to one of the learned clusters, and reduced accordingly to be then matched to the best dictionary signal in the corresponding cluster. It leads to improvements in both memory management complexity and reconstruction speed for parameter maps when compared to traditional dictionary matching, see Section 7.5.

We compare two instances of HD-MED, namely HD-GMM and HD-STM, where the mixture components are respectively set to Gaussian and Student distributions. The computations for the corresponding OEM are detailed in Section 5 of the Appendix. Additional computations and results for mixtures of Laplace distributions (HD-LM) are given in Appendix Section 6. All experiments are performed with a Python code using the JAX library (Bradbury et al., 2018) with Nvidia H100-80gb GPU for training, and V100-32gb GPU for inference except for the dictionary generation part that is done with a mix of Matlab and Python.

7.1. *Undersampled MRF acquisitions.* *In vivo* acquisitions were conducted on 6 healthy volunteers (28±5.5 years old, 3 males and 3 females) using a 32-channel head receiver array on a Philips 3T Achieva dStream MRI at the IRMaGe facility (MAP-IRMaGe protocol, NCT05036629). This study was approved by the local medical ethics committee and informed consent was obtained from each volunteer prior to image acquisition. The imaging pulse sequence was based on an IR-bSSFP acquisition. 260 repetitions were acquired following the parameters proposed in Coudert et al. (2025a). The acquisitions were performed using quadratic variable density spiral sampling (12 interleaves out of 13), matrix
<figure><img src="image_4.png" /><figcaption>FIG 1. HD-MED-based low-dimensional signal matching step (step 2 of Algorithm 1) Green block: each acquired signal is assigned to one of the K HD-MED learned clusters and reduced accordingly using the assigned cluster projection. Blue block: reduced acquired signals are then matched with their closest reduced simulated counterpart in the cluster.</figcaption></figure>

size=192x192x(4-5), voxel size=1.04x1.04x3.00 mm³ for a total scan duration of 2 minutes per slice. While the acquisition time may appear long for an MRF context compared to Gu et al. (2018), where full T1 and T2 maps are generated in 2 minutes, the longer sequences in Coudert et al. (2025a) account for vascular parameters without requiring contrast agents.

Spiral sampling in MRI involves acquiring data in a spiral trajectory through k-space, covering the center first and gradually moving outward, which allows for faster data acquisition and more efficient use of the scanner's gradients. However, this method leads to undersampling noise, manifesting as artifacts, because it captures less information about the image, reducing the ability to accurately reconstruct fine details (Körzdörfer et al., 2019). It results in a number $\tilde{N}_s \approx 140,000 - 180,000$ of in vivo MR signals per subject $s$, each of dimension $M = 260$, and a total number $\tilde{N}$ of approximately $\tilde{N} = 960,000$ signals to be matched for the reconstruction of $L = 6$ parameter maps for 6 subjects.

### 7.2. MRF dictionary.
The MRF sequence used here is a bSSFP-derived sequence. This type of sequence has been preferred because of its sensitivity to local frequency distributions related to microvascular structures in the imaging voxel. The MRF dictionary is generated using the approach described in Coudert et al. (2025a). First, a base dictionary $D_{f_0}$ is simulated using Bloch equations for different combinations of ($T_1, T_2, B_1, \delta f$) sampled over a 4-dimensional regular grid. Simulations were made at a magnetic field strength of 3.0T, on a regular parameter grid made of 20 $T_1$ values (from 200 to 3500ms), 20 $T_2$ values (from 10 to 600ms), 10 $B_1$ values (from 0.7 to 1.2) and 100 frequency offset $\delta f$ values (from -50 to 49 Hz with an increment of 1 Hz), keeping only signals for which $T_1>T_2$, resulting in an initial 390,000 entries dictionary. These Bloch simulations are obtained using an in-house mix of Python and Matlab code based on initial implementation by B. Hargreaves (Hargreaves). The
<figure><img src="image_5.png" /><figcaption>FIG 2. Illustration of the base dictionary expansion. The blue signal corresponding to a non-blood voxel with associated $T_1$, $T_2$, $B_1$ and $\delta f$ values. By convoluting this signal along the $\delta f$ dimension with two different frequency distributions corresponding to different values of CBV and R, the resulting orange and green signals are then obtained for voxels with the same $T_1$, $T_2$, $B_1$ and $\delta f$ values as the blue signal but with different CBV and R values.</figcaption></figure>

relaxometric parameters vary within the dictionary to allow their estimation. $B_1$ is varied to ensure the realism of the dictionary, given the sequence's sensitivity to this parameter, which could otherwise bias the estimation of the other parameters. Finally, varying $\delta f$ values are necessary to compute microvascular contributions, as detailed below.

To build signals that capture additional microvascular network's blood volume (CBV) and mean vessel radius (R) information, we follow the construction of Coudert et al. (2025a). Microvascular network segmentations are used to pre-compute 2500 different intra-voxel frequency distributions centred at $\delta f$ values. A CBV and R value characterize each distribution. Using these pre-computed distributions, new signals are obtained by summing signals from the base dictionary weighted according to each frequency distribution. The result is an expanded 6-dimensional dictionary $\mathcal{D}_f$ of almost $N = 400,000,000$ signals of dimension $M = 260$, encoding for $(T_1, T_2, B_1, \delta f, \text{CBV}, R)$. Figure 2 illustrates how this operation changes a signal in the base dictionary into 2 new signals depending on CBV and R values.

## 7.3. Model selection and initialization.
In practice, the online EM algorithm is run in batches (Nguyen, Forbes and McLachlan, 2020) of size $B = 2048$ signals at each iteration $r$ where $r \in [1:R]$ for typically $R = \lfloor N/B \rfloor$. When $N$ is huge, the $N$ signals are read only once in a single epoch, while for more moderate $N$, several epochs, e.g. 10, can be run with shuffling, making use of the fact that samples are available and can be revisited. Regarding the learning rates, $(\gamma_r)_{r=1}^R$, they are set to $\gamma_r = (1-\epsilon)r^{-0.5-\epsilon'}$ with $\epsilon = 10^{-10}, \epsilon' = 10^{-1}$ and Polyak-Ruppert averaging, following theoretical recommendations. In particular, for this choice, the MSE converges at the optimal rate and the covariance in the CLT is minimal, see Lauand and Meyn (2024) for additional results and a recent literature survey.

The only two other hyperparameters that need to be set are the number $K$ of components and the desired reduced dimension vector $d$. Once initialized, these quantities are considered as fixed and not updated during the online procedure. For mixture models, $K$ can be selected using the Bayesian Information Criterion (BIC) (Schwarz, 1978), which requires running multiple models with varying $K$ values and finding the elbow or the minimum of the BIC curve. The desired dimension vector $d$ is automatically set during an initialization step, which makes $K$ the only parameter that the user needs to set prior to running the OEM algorithm.

The initialization step includes determining the reduced dimension vector $d$. We randomly choose a subset of the dictionary that fits into memory and run a batch full covariance EM algorithm to determine a first estimation of the $\Sigma_k$'s. In practice, a subset of 10,000 signals
was used. The covariance estimations are then decomposed into eigenvalues and eigenvectors to determine $D_k$'s and $A_k$'s. For each $k$, $d_k$ is then determined by applying a scree plot to the eigenvalues using the kneadle algorithm (Satopaa et al., 2011), $a_k$ is then initialized to the first $d_k$ eigenvalues, and the last eigenvalues are averaged to initialize $b_k$. This kind of initialization of $D_k$ and $A_k$ is called spectral initialization and is similar to the one proposed by Hong et al. (2021), proving to be relatively stable compared to other types of initializations.

In Figure 3 (left), we show the different BIC curves for HD-GMM and HD-STM models trained on the dictionary of signals with $K$ varying from 5 to 80. The elbows of the curves are found at $K = 30$ for HD-GMM, and at $K = 25$ for HD-STM. Figures 3 (middle and right) show the different vectors $d$ obtained for varying $K$. The larger the points, the higher the number of components having the reduced dimension indicated on the y-axis. Most of these dimensions are between 5 and 40, meaning that the original $M = 260$ signal dimension can be reduced by a factor of approximately 7 to 10 leading to a reduction from 13 TB to approximately 1.5 TB of signals stored in memory. To check, how much is lost in this reduction, we compute the root mean square error (RMSE) between the dictionary $\mathbf{y}_i$ signals and their reconstructions $\tilde{\mathbf{y}}_i$ from their reduced representations (11), $RMSE = \sqrt{\frac{1}{N}\sum_{i=1}^{N}(\mathbf{y}_i - \tilde{\mathbf{y}}_i)^2}$.

<figure><img src="image_6.png" /><figcaption>FIG 3. Model selection. BIC scores (left) and clusters dimensions, with respect to the number of components $K$, for HD-GMM (middle) and HD-STM (right) applied to signals of dimension $M = 260$. Points sizes reflect the proportion of clusters with a given dimension (y-axis).</figcaption></figure>

Figure 4 shows that the RMSE decreases when $K$ increases and is smaller for HD-STM in particular for small $K$. For the selected $K$ values it represents about 2.5% of the signal.

<figure><img src="image_7.png" /><figcaption>FIG 4. Dictionary reduction loss measured by RMSE over the dictionary signals, for HD-GMM and HD-STM, with respect to the number of components $K$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIG 5. Parameter maps comparison. Average SSIM in [-1,1], the higher the better, and standard deviation over subjects, between HD-GMM, (resp. HD-STM, DRONE) and full matching maps.</figcaption></figure>

Similar plots are shown in Figure 1 of Appendix Section 6 for the Laplace model HD-LM, leading to select $K = 30$ and similar cluster dimension values.
7.4. *Tissue and sensitivity parameters reconstruction.* For the six subjects, we compare the parameter maps obtained with our Algorithm 1 using HD-GMM, HD-STM and with a MRF state-of-the-art neural network approach DRONE (Cohen, Zhu and Rosen, 2018) tailored for $T_1$ and $T_2$, which we adapted for additional parameters, to that obtained with traditional matching, referred to as *full matching* (i.e., matching with uncompressed signals). While *full matching* cannot be considered as the ground truth, it remains the reference method due to its robustness to the highly undersampled MRF acquisitions used in this study, as demonstrated in (Coudert et al., 2025a), despite its high computational and memory demands. Conventionally, the distance to compare two signals for matching is the Euclidean distance.

Table 1 reports for each of the 6 parameters, the distances to full matching parameter estimations, showing the average Mean Absolute Errors (MAE) across voxels for all $n_{subject} = 6$ subjects and all slices. For another comparison that takes into account the image structure, Figure 5 shows the reconstruction quality as measured by the structural similarity index measure (SSIM) (Wang et al., 2004) when compared to full matching, for each parameter and averaged across subjects. The SSIM is a decimal value between -1 and 1, where 1 indicates perfect similarity, 0 indicates no similarity, and -1 indicates perfect anti-correlation.

$$MAE = \frac{1}{n_{\text{subject}}} \sum_{s=1}^{n_{\text{subject}}} \frac{1}{\tilde{N}_s} \sum_{j=1}^{\tilde{N}_s} |\tilde{t}_{j,s} - t_{j,s}^{\text{match}}|,$$
$$SSIM = \frac{1}{n_{\text{subject}}} \sum_{s=1}^{n_{\text{subject}}} SSIM((\tilde{t}_{j,s})_{1 \le j \le \tilde{N}_s}, (t_{j,s}^{\text{match}}))_{1 \le j \le \tilde{N}_s}).$$

<table><caption>TABLE 1
Parameter maps reconstruction. Average MAE and standard deviation, over voxels and subjects, for HD-GMM ($K = 30$) HD-STM ($K = 25$) and DRONE with respect to full matching. Best values in bold.</caption><thead><tr><th>Parameter</th><th><i>T</i><sub>1</sub> (ms)</th><th><i>T</i><sub>2</sub> (ms)</th><th>δ<i>f</i> (Hz)</th><th><i>B</i><sub>1</sub> sensitivity (10<sup>-3</sup>)</th><th>CBV (%)</th><th><i>R</i> (μm)</th></tr></thead><tbody><tr><td>HD-GMM</td><td>462 ± 29</td><td>325 ± 17</td><td>10 ± 1</td><td>92 ± 18</td><td>5 ± 0.8</td><td>2 ± 0.03</td></tr><tr><td>HD-STM</td><td><b>78 ± 14</b></td><td><b>5 ± 2</b></td><td><b>6 ± 1</b></td><td><b>40 ± 5</b></td><td><b>1.5 ± 0.5</b></td><td><b>1.6 ± 0.06</b></td></tr><tr><td>DRONE</td><td>103 ± 8</td><td>351 ± 10</td><td>7 ± 1</td><td>60 ± 4</td><td>7 ± 1</td><td>4 ± 0.2</td></tr></tbody></table>

HD-STM consistently outperforms HD-GMM and DRONE for all parameters, as evidenced by Figures 5 and Table 1. HD-STM achieves superior SSIM and MAE values, indicating better structural correspondence. This advantage is particularly pronounced for parameters such as $T_1$ and $T_2$, where HD-STM demonstrates robust reconstruction capabilities, whereas HD-GMM, and to a lesser extent DRONE, fail to accommodate too much undersampling noise. However, for parameters like CBV and $R$, HD-STM's performance decreases, suggesting that these parameters are inherently more challenging to model accurately. In these cases, DRONE and HD-GMM yield the lowest SSIM values and the highest MAE values, further underscoring their limitations. Overall, while HD-STM proves to be the more reliable method, the significant variability across subjects—particularly for $\delta f$ and CBV—and the observed SSIM drop for $R$ highlight the need for further refinements.

For another quality assessment, Table 2 presents the mean and standard deviations of $T_1$, $T_2$, CBV, and $R$, computed over voxels in white and grey matter ROIs delineated on $T_1$ maps obtained from full matching. Compared to the ranges for healthy subjects reported by Wansapura et al. (1999); Bjørnerud and Emblem (2010); Delphin et al. (2023), HD-STM
produces values that are both more consistent with expected ranges and closer to those obtained using full matching than to those derived from HD-GMM or DRONE. In particular, with HD-GMM (resp. DRONE), $T_2$ (resp. $R$) values significantly depart from the literature reference, as also visible in Figure 6.

The conclusions drawn from Figure 5 and Tables 1–2 are further illustrated in Figure 6. This figure highlights HD-GMM and DRONE’s failure to reconstruct parameters accurately, even when some estimates fall within healthy ranges reported in the literature from Table 2. In contrast, HD-STM consistently delivers parameter estimates for $T_1$, $T_2$, $\delta f$, and $B_1$ that are nearly artifact-free and closely aligned with expected values. Additionally, HD-STM effectively captures primary structures and generates homogeneous maps for vascular parameters such as CBV and $R$ even though minor residual *shim* artifacts remain, likely resulting from significant undersampling in the acquired images. Nevertheless, the presence of these artifacts indicates that HD-STM preserves the physical characteristics of the image, rather than suppressing them. Reconstruction results for HD-LM are shown in Appendix Section 6.1. HD-LM performs much better than HD-GMM and DRONE but worse than HD-STM.

<figure><img src="image_9.png" /><figcaption>FIG 6. $T_1$, $T_2$, $\delta f$, $B_1$ sensitivity, CBV, $R$ maps (columns) with different methods (lines) - full matching, HD-STM ($K = 25$), HD-GMM ($K = 30$), and DRONE. All other slices of this subject are shown in Appendix Section 8, Figure 3.</figcaption></figure>

7.5. *Computation times.* Training is performed on an Nvidia H100-80GB GPU, requiring 26 hours for HD-GMM, 28 hours for HD-STM, and 12 hours for DRONE. For HD models, these times do not include initialization, which consists of computing BIC for different $K$ on a subset of signals. The different $K$ values can be run in parallel so that initialization time mainly depends on the maximum $K$ tested and the number of sub-sampled signals. As an illustration, for $K = 80$ and 10,000 signals, initialization takes 15 minutes. Reconstruction is executed on an Nvidia V100-32GB GPU, taking 3 minutes $\pm$ 2 seconds per slice for both HD-GMM and HD-STM, significantly faster than the 29 minutes $\pm$ 46 seconds required for full matching. Due to memory constraints, as the full uncompressed dictionary exceeds the 32GB GPU memory capacity, an incremental matching approach was implemented. Slice
<table><caption>TABLE 2
Mean $T_1$, $T_2$, CBV, R and standard deviations, in white (WM) and grey (GM) matter ROIs. Values significantly departing from literature (last column) are in red. 1: Cohen, Zhu and Rosen (2018). 2: Wang et al. (2019b)</caption><thead><tr><th>Parameter</th><th>ROI</th><th>Full matching</th><th>HD-GMM</th><th>HD-STM</th><th>DRONE<sup>1</sup></th><th>Literature<sup>2</sup></th></tr></thead><tbody><tr><td rowspan="2">$T_1$(ms)</td><td>WM</td><td>891 ± 6</td><td>650 ± 23</td><td>869 ± 19</td><td>696 ± 14</td><td>690-1100</td></tr><tr><td>GM</td><td>1566 ± 12</td><td>1650 ± 57</td><td>1631 ± 56</td><td>1660 ± 53</td><td>1286-1393</td></tr><tr><td rowspan="2">$T_2$(ms)</td><td>WM</td><td>52 ± 0.0</td><td>400 ± 15</td><td>48 ± 2</td><td>125 ± 5</td><td>56-80</td></tr><tr><td>GM</td><td>95 ± 4</td><td>370 ± 43</td><td>96 ± 12</td><td>135 ± 24</td><td>78-117</td></tr><tr><td rowspan="2">CBV(%)</td><td>WM</td><td>5.4 ± 0.6</td><td>2.0 ± 0.1</td><td>4.8 ± 0.8</td><td>11.5 ± 0.3</td><td>1.7 - 3.6</td></tr><tr><td>GM</td><td>9.2 ± 1.1</td><td>4.2 ± 0.8</td><td>7.6 ± 1</td><td>14.6 ± 0.6</td><td>3.0 - 8.0</td></tr><tr><td rowspan="2">$R$ (μs)</td><td>WM</td><td>5.7 ± 0.1</td><td>5.6 ± 0.1</td><td>5.7 ± 0.1</td><td>0.5 ± 0.01</td><td>6.8 ± 0.3</td></tr><tr><td>GM</td><td>6.2 ± 0.1</td><td>5.4 ± 0.4</td><td>6.1 ± 0.1</td><td>0.6 ± 0.1</td><td>7.3 ± 0.3</td></tr></tbody></table>

reconstruction using DRONE on an Nvidia V100-32GB GPU is notably rapid, completing in under 1 minute. This highlights the efficiency of neural network-based regression methods. However, DRONE exhibits limitations in accurately estimating vascular parameters (Figure 6). In more realistic scenarios, calculations are performed locally by a medical practitioner on a CPU. On an Apple M2 Pro CPU, full matching takes days and is prone to many memory issues that need to be solved, while the HD-GMM and HD-STM variants require only 4.5 hours. This is a strength of HD models over complex neural network approaches, which cannot be easily adapted to CPU and whose running times would be prohibitive.

**8. Conclusion.** In this work, we combine robust latent variable representations, clustering and incremental learning to propose a tractable and accurate way to represent and handle large volumes of potentially heterogeneous high-dimensional data. To our knowledge, this combination and its use as a data compression strategy is novel. The clustering structure of HD-MED allows to address data heterogeneity, for a greater dimensionality reduction without increasing information loss. Incremental learning allows to handle large volumes, resulting in significant reduction in both computational costs and information losses. In terms of implementation, the procedure is flexible and easy to interpret. It depends on two main hyperparameters, the number of clusters *K* and the vector of reduced dimensions *d*, that can be set using conventional model selection criteria or chosen by the user to meet resource constraints, e.g. increasing *K* to reduce the size of clusters or decreasing *d* to increase compression.

As an illustration, we focus on a Magnetic Resonance Fingerprinting (MRF) application. The proposed method drastically reduces the computational time required on standard hardware, such as CPUs, which are commonly used in clinical environments. Beyond this crucial processing time gain, dimensionality reduction has also an interesting impact on patient data acquisition time, as it helps mitigate the effect of noise in fast-acquired *in vivo* signals. This reduction is thus not only a technical achievement but also a significant step forward in making advanced MRF techniques more accessible and usable in daily clinical practice. Currently, the sequence used to acquire these images is in the pre-clinical testing phase on pathological animal, including those with tumors and stroke. With respect to cross-vendor data, since MRF provides quantitative measurements, there are fewer harmonization issues between scanners from different vendors. However, a potential limitation is the noise, which can differ between scanners; our simulated dictionary originally contains denoised signals, which we artificially noised with different values of the signal-to-noise ratio (SNR), but we did not adjust the noise to exactly match that of the acquisition data, also suggesting a potentially good generalization between acquisitions from different scanners.
In addition, although illustrated with a simple matching procedure, which is the current reference in the target MRF application, our procedure can be coupled with other simulation-based inference (SBI) approaches (Cranmer, Brehmer and Louppe, 2020), such as approximate Bayesian computation (ABC) or other neural techniques (Barrier et al., 2024). Like MRF, simulation-based inference has to face two opposite requirements, which are the need for large and high-dimensional simulated data sets to accurately capture information on the physics under study and the issue of handling such large volumes due to computational resource constraints. ABC uses distances between observations and simulations, and inference is based on these distances. A reduced data representation can typically be coupled with ABC techniques proposing automatic summary statistics selection such as in Forbes et al. (2022); Fearnhead and Prangle (2012). It can also be used with SBI neural approaches, e.g. Hägström et al. (2024), when coupled with a regression model for high-dimensional data, see Algorithm 1 in supplementary material. Compared to matching and standard ABC, regression-based procedures are more amortized, which further reduces the reconstruction time. In terms of MRF and other medical imaging applications, future work will thus involve combining HD-MED with a regression model or a neural network, such as in Deleforge, Forbes and Horaud (2015) or Golbabaee et al. (2019). More generally, the possibility to incrementally learn HD-MED models, from high-dimensional heterogeneous data in large volumes, opens the way to a broad range of applications. A first type of applications relates to anomaly detection using a likelihood score, see Lagogiannis et al. (2024) for a recent review. Anomaly detection can be described as an outlier detection problem, where the aim is to discriminate between in- and out-of-distribution samples of a normative distribution learned from a dataset of normal samples. These normal samples can typically be high-dimensional feature vectors summarizing a possibly heterogeneous set of subjects. Second, inverse problems in an SBI context provide numerous applications other than MRF, such as remote sensing tasks in planetary science, extending, for instance, the work by Kugler, Forbes and Douté (2022). Third, HD-MED could also simply be used as an alternative to standard PCA to achieve a better low-dimensional representation of datasets of interest, or simply as a clustering technique (Bouveyron and Brunet-Saumard, 2014).

Lastly, none of the mentioned methods makes use of spatial information, to perform map reconstruction. The idea would be to account for voxel proximity to improve or accelerate prediction. A natural extension of mixture-based methods is to account for spatial information by adding a Markov dependence on the clusters, as done by Deleforge et al. (2015).

**Acknowledgments.** The authors would like to thank the referees, an Associate Editor and the Editor for their constructive comments that improved the quality of this paper.

**Funding.** GO was supported by the AURA region and was granted access to the HPC resources of IDRIS under the allocation 2022-AD011013867 made by GENCI. The MRI facility IRMaGe is partly funded by the program "Investissement d'avenir" run by the French National Research Agency, grant "Infrastructure d'avenir en Biologie et Santé" [ANR-11-INBS-0006]. The project is supported by the French National Research Agency [ANR-20-CE19-0030 MRFUSE].

## SUPPLEMENTARY MATERIAL

Supplementary material includes an additional discussion on the practical benefits of the online EM algorithm. Details on its application to finite mixture models, and in particular to mixtures of Student and Laplace distributions in a high-dimensional setting are also provided. More illustrations on the fingerprinting application can also be found. The code to run the proposed algorithm and reproduce key results is also accessible on
GitHub (https://github.com/geoffroyO/onlineHDEM). The repository includes a notebook, *tutorial_hd_models.ipynb*, which provides a comprehensive tutorial on how to run the different models with simulated data. However, for regulatory and patient privacy reasons, the data used in this study cannot be made publicly available.

## REFERENCES

* ARCHAMBEAU, C. and BACH, F. R. (2008). Sparse probabilistic projections. In *21st International Conference on Neural Information Processing Systems* 73–80.

* ARCHAMBEAU, C., DELANNAY, N. and VERLEYSEN, M. (2006). Robust probabilistic projections. In *Proceedings of the 23rd International conference on machine learning* 33–40.

* ARCHAMBEAU, C., DELANNAY, N. and VERLEYSEN, M. (2008). Mixtures of robust probabilistic principal component analyzers. *Neurocomputing* 71 1274–1282.

* BAEK, J. and MCLACHLAN, G. J. (2011). Mixtures of common t-factor analyzers for clustering high-dimensional microarray data. *Bioinformatics* 27 1269–1276.

* BALZANO, L., CHI, Y. and LU, Y. M. (2018). Streaming PCA and Subspace Tracking: The Missing Data Case. In *Proceedings of the IEEE* 106 1293–1310.

* BARRIER, A., COUDERT, T., DELPHIN, A., LEMASSON, B. and CHRISTEN, T. (2024). MARVEL: MR Fingerprinting with Additional micRoVascular Estimates using bidirectional LSTMs. In *International Conference on Medical Image Computing and Computer-Assisted Intervention* 259–269. Springer.

* BELLAS, A., BOUVEYRON, C., COTTRELL, M. and LACAILLE, J. (2013). Model-based Clustering of High-dimensional Data Streams with Online Mixture of Probabilistic PCA. *Advances in Data Analysis and Classification* 7 281–300.

* BIPIN MEHTA, B., COPPO, S., FRANCES MCGIVNEY, D., IAN HAMILTON, J., CHEN, Y., JIANG, Y., MA, D., SEIBERLICH, N., GULANI, V. and ALAN GRISWOLD, M. (2019). Magnetic resonance fingerprinting: a technical review. *Magnetic resonance in medicine* 81 25–46.

* BJØRNERUD, A. and EMBLEM, K. E. (2010). A fully automated method for quantitative cerebral hemodynamic analysis using DSC–MRI. *Journal of Cerebral Blood Flow & Metabolism* 30 1066–1078.

* BORKAR, V. S. (2009). *Stochastic approximation: a dynamical systems viewpoint* 48. Springer, Gurgaon.

* BOUVEYRON, C. and BRUNET-SAUMARD, C. (2014). Model-based clustering of high-dimensional data: A review. *Computational Statistics & Data Analysis* 71 52–78.

* BOUVEYRON, C., GIRARD, S. and SCHMID, C. (2007). High-dimensional data clustering. *Computational statistics & data analysis* 52 502–519.

* BOUX, F., FORBES, F., ARBEL, J., LEMASSON, B. and BARBIER, E. L. (2021). Bayesian inverse regression for vascular magnetic resonance fingerprinting. *IEEE Transactions on Medical Imaging* 40 1827–1837.

* BRADBURY, J., FROSTIG, R., HAWKINS, P., JOHNSON, M., LEARY, C., MACLAURIN, D., NECULA, G., PASZKE, A., VANDERPLAS, J., WANDERMAN-MILNE, S. and ZHANG, Q. (2018). JAX: composable transformations of Python+NumPy programs. http://github.com/google/jax.

* CABINI, R. F., BARZAGHI, L., CICOLARI, D., AROSIO, P., CARRAZZA, S., FIGINI, S., FILIBIAN, M., GAZZANO, A., KRAUSE, R., MARIANI, M. et al. (2024). Fast deep learning reconstruction techniques for preclinical magnetic resonance fingerprinting. *NMR in Biomedicine* 37 e5028.

* CAMBANIS, S., HUANG, S. and SIMONS, G. (1981). On the theory of elliptically contoured distributions. *Journal of Multivariate Analysis* 11 368–385.

* CAPPÉ, O. and MOULINES, E. (2009). On-line expectation–maximization algorithm for latent data models. *Journal of the Royal Statistical Society Series B: Statistical Methodology* 71 593–613.

* CAULEY, S. F., SETSOMPOP, K., MA, D., JIANG, Y., YE, H., ADALSTEINSSON, E., GRISWOLD, M. A. and WALD, L. L. (2015). Fast group matching for MR fingerprinting reconstruction. *Magnetic resonance in medicine* 74 523–528.

* CHRISTEN, T., BOLAR, D. S. and ZAHARCHUK, G. (2013). Imaging brain oxygenation with MRI using blood oxygenation approaches: methods, validation, and clinical applications. *American journal of neuroradiology* 34 1113–1123.

* CHRISTEN, T., PANNETIER, N., NI, W. W., QIU, D., MOSELEY, M. E., SCHUFF, N. and ZAHARCHUK, G. (2014). MR vascular fingerprinting: A new approach to compute cerebral blood volume, mean vessel radius, and oxygenation maps in the human brain. *Neuroimage* 89 262–270.

* COHEN, O., ZHU, B. and ROSEN, M. S. (2018). MR fingerprinting deep reconstruction network (DRONE). *Magnetic resonance in medicine* 80 885–894.

* COIFMAN, R. R. and LAFON, S. (2006). Diffusion maps. *Applied and Computational Harmonic Analysis* 21 5–30.
COUDERT, T., DELPHIN, A., BARRIER, A., LEGRIS, L., WARNKING, J. M., LAMALLE, L., DONEVA, M., LEMASSON, B., BARBIER, E. L. and CHRISTEN, T. (2025a). Relaxometry and contrast-free cerebral microvascular quantification using balanced steady-state free precession MR fingerprinting. *Magnetic Resonance in Medicine* **94** 302-316.

COUDERT, T., DELPHIN, A., BARRIER, A., BARBIER, E., LEMASSON, B., WARNKING, J. and CHRISTEN, T. (2025b). MR Fingerprinting for Imaging Brain Hemodynamics and Oxygenation. *Journal of Magnetic Resonance Imaging*.

CRANMER, K., BREHMER, J. and LOUPPE, G. (2020). The frontier of simulation-based inference. *Proceedings of the National Academy of Sciences* **117** 201912789.

D'ASPREMONT, A., ELGHAOUI, L., JORDAN, M. I. and LANCKRIET, G. (2007). A direct formulation for sparse PCA using semidefinite programming. *SIAM Review* **49** 434–48.

DELEFORGE, A., FORBES, F. and HORAUD, R. (2015). High-Dimensional Regression with Gaussian Mixtures and Partially-Latent Response Variables. *Statistics & Computing* **25** 893-911.

DELEFORGE, A., FORBES, F., BA, S. O. and HORAUD, R. (2015). Hyper-Spectral Image Analysis With Partially Latent Regression and Spatial Markov Dependencies. *IEEE J. Sel. Top. Signal Process.* **9** 1037–1048.

DELPHIN, A., COUDERT, T., FAN, A., MOSELEY, M. E., ZAHARCHUK, G. and CHRISTEN, T. (2023). MR Vascular Fingerprinting with 3D realistic blood vessel structures and machine learning to assess oxygenation changes in human volunteers. In *2023 ISMRM & ISMRT Annual Meeting & Exhibition*.

DEMPSTER, A. P., LAIRD, N. M. and RUBIN, D. B. (1977). Maximum likelihood from incomplete data via the EM algorithm. *Journal of the royal statistical society: series B (methodological)* **39** 1–22.

DIEULEVEUT, A., FORT, G., MOULINES, E. and WAI, H.-T. (2023). Stochastic Approximation Beyond Gradient for Signal Processing and Machine Learning. *IEEE Transactions on Signal Processing* **71** 3117-3148.

FEARNHEAD, P. and PRANGLE, D. (2012). Constructing summary statistics for approximate Bayesian computation: semi-automatic approximate Bayesian computation. *Journal of the Royal Statistical Society Series B: Statistical Methodology* **74** 419–474.

FORBES, F., NGUYEN, H. D., NGUYEN, T. and ARBEL, J. (2022). Summary statistics and discrepancy measures for approximate Bayesian computation via surrogate posteriors. *Statistics and Computing* **32** 85.

FORT, G., MOULINES, E. and WAI, H.-T. (2020). A stochastic path-integrated differential estimator expectation maximization algorithm. In *Proceedings of the 34th Conference on Neural Information Processing Systems (NeurIPS)*.

FREITAS LOPES, H. and WEST, M. (2004). Bayesian Model Assessment in Factor Analysis. *Statistica Sinica* **14** 41-67.

GHAHRAMANI, Z. and BEAL, M. J. (1999). Variational inference for Bayesian mixtures of factor analysers. In *Proceedings of the 13th International Conference on Neural Information Processing Systems. NIPS'99* 449–455. MIT Press, Cambridge, MA, USA.

GHAHRAMANI, Z. and HINTON, G. E. (1997). The EM algorithm for mixtures of factor analyzers Technical Report, University of Toronto.

GILMAN, K., HONG, D., FESSLER, J. A. and BALZANO, L. (2023). Streaming Probabilistic PCA for Missing Data with Heteroscedastic Noise. *arXiv preprint arXiv:2310.06277*.

GIRAUD, C. (2014). *Introduction to high-dimensional statistics*. Chapman and Hall/CRC.

GOLBABAEE, M., CHEN, D., GÓMEZ, P. A., MENZEL, M. I. and DAVIES, M. E. (2019). Geometry of deep learning for magnetic resonance fingerprinting. In *ICASSP 2019-2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)* 7825-7829. IEEE.

GOMEZ, E., GOMEZ-VILLEGAS, M. and MARIN, J. (2008). Multivariate Exponential Power Distributions as Mixtures of Normal Distributions with Bayesian Applications. *Communications in Statistics-theory and Methods* **37**.

GÓMEZ-SÁNCHEZ-MANZANO, E., GÓMEZ-VILLEGAS, M. and MARÍN, J. (2006). Sequences of elliptical distributions and mixtures of normal distributions. *Journal of multivariate analysis* **97** 295–310.

GU, Y., WANG, C. Y., ANDERSON, C. E., LIU, Y., HU, H., JOHANSEN, M. L., MA, D., JIANG, Y., RAMOS-ESTEBANEZ, C., BRADY-KALNAY, S. et al. (2018). Fast magnetic resonance fingerprinting for dynamic contrast-enhanced studies in mice. *Magnetic resonance in medicine* **80** 2681–2690.

GU, Y., PAN, Y., FANG, Z., MA, L., ZHU, Y., ANDROJNA, C., ZHONG, K., YU, X. and SHEN, D. (2024). Deep learning-assisted preclinical MR fingerprinting for sub-millimeter T1 and T2 mapping of entire macaque brain. *Magnetic Resonance in Medicine* **91** 1149–1164.

GUO, Y. and BONDELL, H. (2023). On robust probabilistic principal component analysis using multivariate t-distributions. *Communications in Statistics - Theory and Methods* **52** 8261-8279.

HÄGGSTRÖM, H., RODRIGUES, P. L., OUDOUMANESSAH, G., FORBES, F. and PICCHINI, U. (2024). Fast, accurate and lightweight sequential simulation-based inference using Gaussian locally linear mappings. *arXiv preprint arXiv:2403.07454*.
HALKO, N., MARTINSSON, P.-G. and TROPP, J. A. (2011). Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions. *SIAM review* **53** 217–288.

HARGREAVES, B. Bloch Equation Simulator. http://www-mrsrl.stanford.edu/brian/blochsim/.

HONG, D., GILMAN, K., BALZANO, L. and FESSLER, J. A. (2021). HePPCAT: Probabilistic PCA for data with heteroscedastic noise. *IEEE Transactions on Signal Processing* **69** 4819–4834.

HONG, D., YANG, F., FESSLER, J. A. and BALZANO, L. (2023). Optimally Weighted PCA for High-Dimensional Heteroscedastic Data. *SIAM Journal on Mathematics of Data Science* **5** 222–250.

JOLLIFFE, I. T. and CADIMA, J. (2016). Principal component analysis: a review and recent developments. *Philosophical transactions of the royal society A: Mathematical, Physical and Engineering Sciences* **374** 20150202.

KARIMI, B., MIASOJEDOW, B., MOULINES, E. and WAI, H.-T. (2019a). Non-asymptotic analysis of biased stochastic approximation scheme. *Proc. Mach. Learn. Res.* **99** 1–31.

KARIMI, B., WAI, H.-T., MOULINES, R. and LAVIELLE, M. (2019b). On the global convergence of (fast) incremental expectation maximization methods. In *Proceedings of the 33rd Conference on Neural Information Processing Systems (NeurIPS)*.

KELKER, D. (1970). Distribution theory of spherical distributions and a location-scale parameter generalization. *Sankhya: The Indian Journal of Statistics, Series A* **32** 419–430.

KÖRZDÖRFER, G., PFEUFFER, J., KLUGE, T., GEBHARDT, M., HENSEL, B., MEYER, C. H. and NITTKA, M. (2019). Effect of spiral undersampling patterns on FISP MRF parameter maps. *Magnetic resonance imaging* **62** 174–180.

KOTZ, S., KOZUBOWSKI, T. and PODGORSKI, K. (2001). *The Laplace Distribution and Generalizations: A Visit with Applications to Communications, Economics, Engineering, and Finance. Progress in Mathematics*. Birkhäuser Boston.

KOTZ, S. and NADARAJAH, S. (2004). *Multivariate t-distributions and their applications*. Cambridge University Press.

KUGLER, B., FORBES, F. and DOUTÉ, S. (2022). Fast Bayesian inversion for high dimensional inverse problems. *Statistics and Computing* **32**.

KUHN, E., MATIAS, C. and REBAFAK, T. (2020). Properties of the stochastic approximation EM algorithm with mini-batch sampling. *Statistics and Computing* **30** 1725–1739.

LAGOGIANNIS, I., MEISSEN, F., KAISSIS, G. and RUECKERT, D. (2024). Unsupervised Pathology Detection: A Deep Dive Into the State of the Art. *IEEE Transactions on Medical Imaging* **43** 241–252.

LANGE, K. L., LITTLE, R. J. A. and TAYLOR, J. M. G. (1989). Robust Statistical Modeling Using the t Distribution. *Journal of the American Statistical Association* **84** 881–896.

LANGE, K. and SINSHEIMER, J. S. (1993). Normal/Independent Distributions and Their Applications in Robust Regression. *Journal of Computational and Graphical Statistics* **2** 175–198.

LAUAND, C. K. and MEYN, S. (2024). Revisiting Step-Size Assumptions in Stochastic Approximation.

LI, P. and HU, Y. (2023). Learned Tensor Low-CP-Rank and Bloch Response Manifold Priors for Non-Cartesian MRF Reconstruction. *IEEE Transactions on Medical Imaging*.

LI, P. and HU, Y. (2024). Deep magnetic resonance fingerprinting based on Local and Global Vision Transformer. *Medical Image Analysis* **95** 103198.

MA, D., GULANI, V., SEIBERLICH, N., LIU, K., SUNSHINE, J. L., DUERK, J. L. and GRISWOLD, M. A. (2013). Magnetic resonance fingerprinting. *Nature* **495** 187–192.

MAIRE, F., MOULINES, E. and LEFEBVRE, S. (2017). Online EM for functional data. *Computational Statistics and Data Analysis* **111** 27–47.

MCGIVNEY, D. F., PIERRE, E., MA, D., JIANG, Y., SAYBASILI, H., GULANI, V. and GRISWOLD, M. A. (2014). SVD compression for magnetic resonance fingerprinting in the time domain. *IEEE transactions on medical imaging* **33** 2311–2322.

MCGIVNEY, D. F., BOYACIOĞLU, R., JIANG, Y., POORMAN, M. E., SEIBERLICH, N., GULANI, V., KEENAN, K. E., GRISWOLD, M. A. and MA, D. (2020). Magnetic resonance fingerprinting review part 2: Technique and directions. *Journal of Magnetic Resonance Imaging* **51** 993–1007.

MCLACHLAN, G. J. and KRISHNAN, T. (2008). *The EM Algorithm and Extensions*. Wiley 2nd edition.

MCLACHLAN, G. J., PEEL, D. and BEAN, R. (2003). Modelling high-dimensional data by mixtures of factor analyzers. *Computational Statistics & Data Analysis* **41** 379–388.

MENG, X.-L. and RUBIN, D. B. (1993). Maximum likelihood estimation via the ECM algorithm: A general framework. *Biometrika* **80** 267–278.

MONGA, A., SINGH, D., DE MOURA, H. L., ZHANG, X., ZIBETTI, M. V. and REGATTE, R. R. (2024). Emerging Trends in Magnetic Resonance Fingerprinting for Quantitative Biomedical Imaging Applications: A Review. *Bioengineering* **11** 236.

NERI, J., DEPALLE, P. and BADEAU, R. (2021). Approximate Inference and Learning of State Space Models With Laplace Noise. *IEEE Transactions on Signal Processing* **69** 3176–3189.
NGUYEN, H. D., FORBES, F. and McLACHLAN, G. J. (2020). Mini-batch learning of exponential family finite mixture models. *Statistics and Computing* **30** 731-748.

NGUYEN, H. and FORBES, F. (2022). Global implicit function theorems and the online expectation-maximisation algorithm. *Australian & New Zealand Journal of Statistics* **64**.

OUDOUMANESSAH, G., COUDERT, T., MEYER, L., DELPHIN, A., CHRISTEN, T., DOJAT, M., LARTIZIEN, C. and FORBES, F. (2025). Cluster globally, Reduce locally: Scalable efficient dictionary compression for magnetic resonance fingerprinting. In *International Symposium in Biomedical Imaging* 1-5.

POORMAN, M. E., MARTIN, M. N., MA, D., MCGIVNEY, D. F., GULANI, V., GRISWOLD, M. A. and KEENAN, K. E. (2020). Magnetic resonance fingerprinting Part 1: Potential uses, current challenges, and recommendations. *Journal of Magnetic Resonance Imaging* **51** 675-692.

RONNEBERGER, O., FISCHER, P. and BROX, T. (2015). U-net: Convolutional networks for biomedical image segmentation. In *Medical image computing and computer-assisted intervention–MICCAI 2015: 18th international conference, Munich, Germany, October 5-9, 2015, proceedings, part III* **18** 234–241. Springer.

SATOPAA, V., ALBRECHT, J., IRWIN, D. and RAGHAVAN, B. (2011). Finding a" kneedle" in a haystack: Detecting knee points in system behavior. In *2011 31st international conference on distributed computing systems workshops* 166–171. IEEE.

SCHOTT, J. R. (2016). *Matrix analysis for statistics*. Wiley probability & statistics, Hoboken.

SCHWARZ, G. (1978). Estimating the dimension of a model. *The Annals of Statistics* 461–464.

SOYAK, R., NAVRUZ, E., ERSOY, E. O., CRUZ, G., PRIETO, C., KING, A. P., UNAY, D. and OKSUZ, I. (2021). Channel attention networks for robust MR fingerprint matching. *IEEE Transactions on Biomedical Engineering* **69** 1398–1405.

TIPPAREDDY, C., ZHAO, W., SUNSHINE, J. L., GRISWOLD, M., MA, D. and BADVE, C. (2021). Magnetic resonance fingerprinting: an overview. *European Journal of Nuclear Medicine and Molecular Imaging* **48** 4189–4200.

TIPPING, M. E. and BISHOP, C. M. (1999a). Probabilistic principal component analysis. *Journal of the Royal Statistical Society B: Statistical Methodology* **61** 611–622.

TIPPING, M. E. and BISHOP, C. M. (1999b). Mixtures of probabilistic principal component analyzers. *Neural computation* **11** 443–482.

TIPPING, M. E. and BISHOP, C. M. (1999c). Probabilistic principal component analysis. *Journal of the Royal Statistical Society Series B: Statistical Methodology* **61** 611–622.

TIPPING, M. E. and BISHOP, C. M. (1999d). Mixtures of probabilistic principal component analyzers. *Neural computation* **11** 443–482.

ULLAH, I., HASSAN, A. M., SAAD, R. M. and OMER, H. (2023). GPU accelerated grouped magnetic resonance fingerprinting using clustering techniques. *Magnetic Resonance Imaging* **97** 13–23.

VASWANI, A., SHAZEER, N., PARMAR, N., USZKOREIT, J., JONES, L., GOMEZ, A. N., KAISER, Ł. and POLOSUKHIN, I. (2017). Attention is all you need. *Advances in neural information processing systems* **30**.

WANG, Z., BOVIK, A. C., SHEIKH, H. R. and SIMONCELLI, E. P. (2004). Image Quality Assessment: From Error Visibility to Structural Similarity. *IEEE Transactions on Image Processing* **13** 600-612.

WANG, Z., ZHANG, J., CUI, D., XIE, J., LYU, M., HUI, E. S. and WU, E. X. (2019a). Magnetic Resonance Fingerprinting Using a Fast Dictionary Searching Algorithm: MRF-ZOOM. *IEEE trans. bio-medical engineering* **66** 1526-1535.

WANG, C. Y., COPPO, S., MEHTA, B. B., SEIBERLICH, N., YU, X. and GRISWOLD, M. A. (2019b). Magnetic resonance fingerprinting with quadratic RF phase for measurement of T2* simultaneously with δf, T1, and T2. *Magnetic resonance in medicine* **81** 1849–1862.

WANSAPURA, J. P., HOLLAND, S. K., DUNN, R. S. and BALL JR, W. S. (1999). NMR relaxation times in the human brain at 3.0 tesla. *Journal of Magnetic Resonance Imaging: An Official Journal of the International Society for Magnetic Resonance in Medicine* **9** 531–538.

WEN, Z. and YIN, W. (2013). A feasible method for optimization with orthogonality constraints. *Mathematical Programming* **142** 397–434.

XU, A. S., BALZANO, L. and FESSLER, J. A. (2023). HeMPPCAT: mixtures of probabilistic principal component analysers for data with heteroscedastic noise. In *IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)* 1–5.

YANG, M., MA, D., JIANG, Y., HAMILTON, J., SEIBERLICH, N., GRISWOLD, M. A. and MCGIVNEY, D. (2018). Low rank approximation methods for MR fingerprinting with large scale dictionaries. *Magnetic resonance in medicine* **79** 2392–2400.

YE, H., CAULEY, S. F., GAGOSKI, B., BILGIC, B., MA, D., JIANG, Y., DU, Y. P., GRISWOLD, M. A., WALD, L. L. and SETSOMPOP, K. (2017). Simultaneous multislice magnetic resonance fingerprinting (SMS-MRF) with direct-spiral slice-GRAPPA (ds-SG) reconstruction. *Magnetic resonance in medicine* **77** 1966–1974.

ZOU, H., HASTIE, T. and TIBSHIRANI, R. (2006). Sparse Principal Component Analysis. *Journal of Computational and Graphical Statistics* **15** 265–286.
# SUPPLEMENTARY MATERIAL
## SCALABLE MAGNETIC RESONANCE FINGERPRINTING: INCREMENTAL INFERENCE OF HIGH DIMENSIONAL ELLIPTICAL MIXTURES FROM LARGE DATA VOLUMES

**1. Practical benefit of online EM vs batch EM.** Dictionary-based techniques aim at being as amortized as possible, meaning that most of the computational burden is performed offline, ahead of the task of interest, which can then be performed straightforwardly for each new patient. Thus, in principle, the standard batch EM could be executed once, using powerful cloud computing resources, without resorting to a more frugal incremental EM. However, there are several reasons why this may not be the case in practice:

* For large dictionaries, the resources needed for a batch EM are huge, so that not all physicians / hospitals / researchers can afford the resource and time costs. A useful dictionary exceeds in size several terabytes in memory, which requires an excessive amount of RAM to process efficiently. On typical computer clusters, it is hard to find a single node with sufficient RAM, necessitating the use of multiple nodes and thus, resorting to implementing a distributed EM algorithm. The incremental EM setting makes the implementation simpler.

* In addition, building a dictionary can be an experimental design task. For practitioners, the batch EM cost can be detrimental and prevent the necessary repetitive trials to design an appropriate dictionary for a given task. Note that the dictionary considered in the paper is already 400 million signals for a size of 13 terabytes, which, for 6 parameters, represents approximately a grid made with only 27 values for each parameter. Although the hope is that a given dictionary be used for several patients, depending on the study, other parameter values may need to be included or changed, so that the dictionary may have to be changed or be completed with more signals. In the latter case, the incremental EM allows either to re-learn a dictionary representation from scratch, without excessive cost, or to pursue the previous learning with the newly added dictionary elements. In other words, incremental EM enables to refine a dictionary on the fly by adding more simulated signals—such as transitioning to a finer grid—very easily. We do not need to re-run the complete batch EM algorithm but can simply pass mini-batches of new simulated signals incrementally to the model.

* Finally, beyond computational costs, incremental/stochastic methods often enable better generalization. This is particularly the case when Stochastic Gradient Descent (SGD) is compared to Gradient Descent (GD) in deep learning (Amir, Koren and Livni, 2021). In the case of the EM algorithm, experimental results also show that this can hold. For instance, Figure 1 of Gilman et al. (2023) demonstrates that when only 50% of the dataset is observed, the incremental EM achieves better results in terms of log-likelihood compared to the batch version of EM (HePPCAT).

To summarize, compared to standard EM, the proposed incremental EM allows 1) much more flexible experimental design and investigations, and 2) better learning performance. So doing it broadens the applicability and accuracy of MRF in general, and in particular beyond the settings studied in the paper.

**2. Online EM vs Stochastic Gradient Descent (SGD).** SGD can be used to estimate mixtures but this is not yet very popular. Both EM and GD approaches have strengths. Gradient approaches do not need to compute E and M steps and can thus easily be adapted from one model to another. However, although automatic differentiation is a powerful tool, its use when
there are constraints on the parameters is not straightforward and requires reparametrization and sometimes over-reparametrization, typically when parameters are proportions that sum to 1, covariance matrices or orthogonal matrices. See, for instance, an elegant reparametrization used by Hosseini and Sra (2015) that handle covariance matrices by switching to a Riemannian manifold. Handling orthogonality constraints, like for matrix $D$ is likely to be more involved. In contrast, incremental EM, which works on sufficient statistics, does not need smoothness assumptions and gradient computations, see also Dieuleveut et al. (2023) for a recent analysis of stochastic non-gradient algorithms.

In addition, GD/SGD is extremely sensitive to initialization, even more than EM for which initialization can also be important. EM avoids certain poor local minima that GD/SGD might fall into, as shown in analyses of two-component mixtures (Zhang, Poupart and Trimponias, 2020). EM convergence is generally faster. EM can be seen as a preconditioned version of GD, see e.g. Xu and Jordan (1996).

Another feature is that SGD is more memory - and time - intensive than incremental EM, primarily because automatic differentiation tools like autodiff construct and store computational graphs to compute gradients, increasing overhead for large-scale, high-dimensional problems. In contrast, incremental EM can often leverage closed-form updates for the E- and M-steps, avoiding the need for such graphs.

EM-based methods, including incremental versions, are then frequently more robust and efficient for mixture estimation. However, Geppert and Pfülb (2021) compare incremental EM with SGD. To allow SGD to escape local optima, they propose an annealing procedure. Their conclusion is that both SGD and incremental EM perform similarly but SGD can outperform incremental EM in high dimensions. In our work, high dimensionality is handled via our HD-MED model and we chose to benefit from the other strengths of EM.

**3. Online EM for finite mixture models.** This section completes Section 6.1 of the main paper. We demonstrate that the online EM (OEM) algorithm for a finite mixture model can be derived from the OEM algorithm for a single component. We denote the mixture parameters as $\Theta_M$, which include the pairs $(\pi_k, \Theta_k)$ for $k = 1, \dots, K$, where $\Theta_k$ are the component-wise parameters of the finite mixture model. Furthermore, we consider the pair $(\mathbf{Y}^T, Z)$, where $Z = 1:K$ is a latent variable such that $\mathbb{P}(Z = k) = \pi_k$.

We assume that the complete-data likelihood for each component can be expressed in the exponential family form, as given by Equation (12) in the main paper, for each $f_c(\mathbf{y}, z; \Theta_k)$. Following Nguyen and Forbes (2022), the complete-data likelihood for the mixture model $f_c(\mathbf{y}, z, \Theta_M)$ can be written as:

$$ (1) \qquad f_c(\mathbf{y}, z, \Theta_M) = h(\mathbf{y}, z) \exp \left[ \sum_{\xi=1}^{K} \left[ \log \pi_{\xi} + [\mathbf{s}(\mathbf{y}, z)]^T \phi(\Theta_{\xi}) - \psi(\Theta_{\xi}) \right] \right] \\ = h_M(\mathbf{y}, z) \exp \left[ [\mathbf{s}_M(\mathbf{y}, z)]^T \phi_M(\Theta_M) - \psi_M(\Theta_M) \right], $$
where $h_M = h, \psi_M = 0$,

$$ (2) \qquad \mathbf{s}_M(\mathbf{y}, z) = \begin{bmatrix} \mathbf{1}_{z=1} \\ \mathbf{1}_{z=1} \mathbf{s}(\mathbf{y}, z) \\ \vdots \\ \mathbf{1}_{z=K} \\ \mathbf{1}_{z=K} \mathbf{s}(\mathbf{y}, z) \end{bmatrix}, \quad \text{and} \quad \phi_M(\theta_M) = \begin{bmatrix} \log \pi_1 - \psi(\Theta_1) \\ \phi(\Theta_1) \\ \vdots \\ \log \pi_K - \psi(\Theta_K) \\ \phi(\Theta_K) \end{bmatrix}. $$

We can conveniently write the vector $\mathbf{s}_M$ that matches the dimension of $\mathbf{s}_M(\mathbf{y}, z)$ as

$$ (3) \qquad \mathbf{s}_M^T = (s_{01}, \mathbf{s}_{M1}^T, \dots, s_{0K}, \mathbf{s}_{MK}^T), $$
where $s_{Mz}$ corresponds to one component of the exponential family form of the mixture model distribution. Via Equation (2), the objective function $Q_M$ for the mixture complete-data likelihood as defined in Equation (14) of the main paper can be written as

$$ (4) \qquad Q_M(s_M, \Theta_M) = s_M^T \phi_M(\Theta_M) = \sum_{z=1}^{K} s_{0z} (\log \pi_z - \psi(\Theta_z)) + s_{Mz}^T \phi(\Theta_z). $$

Whatever the form of the component p.d.f, the maximisation with respect to $\pi_z$ yields

$$ (5) \qquad \bar{\pi}_z(s_M) = \frac{s_{0z}}{\sum_{\xi=1}^{K} s_{0\xi}}. $$

Then, for each $z \in \{1...K\}$,

$$ \begin{align} \frac{\partial Q_M}{\partial \Theta_z}(s_M, \theta_M) &= -s_{0z} \frac{\partial \psi}{\partial \Theta_z}(\Theta_z) + \mathbf{J}_{\phi}(\Theta_z)s_{Mz} \nonumber \\ &= s_{0z} \left( \mathbf{J}_{\phi}(\Theta_z) \left[ \frac{s_{Mz}}{s_{0z}} \right] - \frac{\partial \psi}{\partial \Theta_z} \right) \nonumber \\ (6) \qquad &= s_{0z} \frac{\partial Q}{\partial \Theta_z} \left( \left[ \frac{s_{Mz}}{s_{0z}} \right], \Theta_z \right), \end{align} $$

where $Q$ is the objective function of form Equation (14) from main paper corresponding to a single component p.d.f. Since $s_{0z} > 0$, for all $z \in \{1...K\}$, it follows that the maximization of $Q_M$ can be conducted by solving

$$ (7) \qquad \frac{\partial Q}{\partial \Theta_z} \left( \left[ \frac{s_{Mz}}{s_{0z}} \right], \Theta_z \right) = 0, $$

with respect to $\Theta_z$ for each $z$. Therefore, it is enough to show that for each component $z$ of the finite mixture model that there exists a root of the equation above denoted as $\bar{\Theta}_z(s_{Mz}/s_{0z})$ in order to find a solution for the maximizer of the mixture objective $Q_M$. Then we can set

$$ (8) \qquad \bar{\Theta}_M(s_M) = \begin{bmatrix} \bar{\pi}_1(s_M) \\ \bar{\Theta}_1(s_{M1}/s_{01}) \\ \vdots \\ \bar{\pi}_K(s_M) \\ \bar{\Theta}_K(s_{MK}/s_{0K}) \end{bmatrix}, $$

where the form of $\bar{\pi}_K(s_M)$ is given in (5).

To complete the OEM algorithm, we need to specify the quantity $\bar{s}_M(y; \Theta_M) = E[s_M(y, Z) | y = y; \Theta_M]$ as defined in Equation (13) of the main paper. Using the definition of $s_M(y, z)$ in Equation (2), we compute,

$$ (9) \qquad r_k = E[\mathbf{1}_{\{Z=k\}} | \mathbf{Y} = \mathbf{y}; \mathbf{\Theta}_M] = \frac{\pi_k f_k(\mathbf{y})}{\sum_{\xi=1}^{K} \pi_k f_k(\mathbf{y})}, $$

and

$$ (10) \qquad E[\mathbf{1}_{\{Z=k\}} s(\mathbf{Y}, Z) | \mathbf{Y} = \mathbf{y}; \mathbf{\Theta}_M] = r_k E[s(\mathbf{Y}) | \mathbf{Y} = \mathbf{y}, Z = k; \mathbf{\Theta}_M], $$
where $f_k$ is defined in Equation (4) of the main paper as the conditional pdf of the $k^{th}$ mixture component.

Thus, in the finite mixture framework, deriving the OEM algorithm for a single component suffices to deduce the OEM algorithm for the entire mixture. In the next section, we derive the OEM algorithm for a single component of HD-MED. For simplicity, we omit further mention of the form of $h$ in the exponential family representation, as it does not affect the OEM algorithm derivation.

### 4. HD-ED exponential family form: Proof of Proposition 6.1.

**PROOF.** Recall that $\tilde{D}$ denotes the $d$ first columns of $D$ completed with zero columns and $\overline{D} = D - \tilde{D}$. The following multivariate normal distribution can be decomposed as follows,

$$ \mathcal{N}_M \left( \mathbf{y}; \boldsymbol{\mu}, \frac{DAD^T}{w} \right) = \prod_{m=1}^{d} \mathcal{N}_1 \left( [\tilde{D}^T(\mathbf{y} - \boldsymbol{\mu})]_m; 0, \frac{a_m}{w} \right) \prod_{m=d+1}^{M} \mathcal{N}_1 \left( [\overline{D}^T(\mathbf{y} - \boldsymbol{\mu})]_m; 0, \frac{b}{w} \right), $$

where $[\cdot]_m$ denotes the $m^{th}$ element of the input vector. For $m \le d$ it comes,

$$ [\tilde{D}^T(\mathbf{y} - \boldsymbol{\mu})]_m^2 = \mathrm{vec}(\mathbf{d}_m \mathbf{d}_m^T)^T \mathrm{vec}(\mathbf{y}\mathbf{y}^T) + \mathrm{vec}(\mathbf{d}_m \mathbf{d}_m^T)^T \mathrm{vec}(\boldsymbol{\mu}\boldsymbol{\mu}^T) - 2\boldsymbol{\mu}^T \mathbf{d}_m \mathbf{d}_m^T \mathbf{y}. $$

And for $d < m \le M$ we have:

$$ \prod_{m=d+1}^{M} \mathcal{N}_1 \left( [\overline{D}^T(\mathbf{y} - \boldsymbol{\mu})]_m; 0, \frac{b}{w} \right) = \frac{1}{(2\pi b)^{\frac{M-d}{2}}} \exp \left( -\frac{w}{2b} \sum_{m=d+1}^{M} \overline{D}^T(\mathbf{y} - \boldsymbol{\mu}) \right) $$

It is then possible to remove the dependence on $\overline{D}$ by noting that

$$ \sum_{m=d+1}^{M} [\overline{D}^T(\mathbf{y} - \boldsymbol{\mu})]_m^2 = \|\mathbf{y} - P(\mathbf{y})\|^2 \\ = \mathbf{y}^T \mathbf{y} + \boldsymbol{\mu}^T \boldsymbol{\mu} - 2\mathbf{y}^T \boldsymbol{\mu} - \sum_{m=1}^{d} [\tilde{D}^T(\mathbf{y} - \boldsymbol{\mu})]_m^2, $$

using that $P(\mathbf{y}) = \tilde{D}\tilde{D}^T(\mathbf{y} - \boldsymbol{\mu}) + \boldsymbol{\mu}$. Rewriting expressions with the vec operator and combining both parts where $m \le d$ and $d \le m \le M$ gives the exponential family form for the complete data likelihood (18) in main paper. □

### 5. Gamma-distributed mixing variable.
In Section 6 of the main paper we proposed an online E-step with incomplete calculations that depend on the exponential family form of mixing variable $W$. In this section we propose to specify these calculations when $W$ follows a Gamma distribution $f_{\theta}(w) = G(w; \alpha, \beta)$ with the two parameters $\theta = (\alpha, \beta)$ in $\mathbb{R}^{+2}$. More specifically, we need to compute the expectation (19), $E[W|\mathbf{Y}; \mathbf{\Theta}]$, from the main paper, and $E[s_w(W)|\mathbf{Y}; \mathbf{\Theta}]$. Then we need to derive the root equations for $\alpha$ and $\beta$. This Gamma-distributed mixing $W$ encompasses several distributions as indicated in Table 1.

The Gamma distribution has density $G(w; \alpha, \beta) = w^{\alpha-1} \frac{\beta^{\alpha} e^{-\beta w}}{\Gamma(\alpha)}$ with $\Gamma$ the Gamma function. This distribution is easily put in the exponential family form with

$$ \begin{array}{l} (11) \qquad \mathbf{s}_w(w) = \begin{bmatrix} w \\ \log w \end{bmatrix}, \quad \phi_w(\alpha, \beta) = \begin{bmatrix} -\beta \\ \alpha \end{bmatrix}, \\ \qquad \psi_w(\alpha, \beta) = \log(\Gamma(\alpha)) - \alpha \log(\beta). \end{array} $$

The next proposition gives the necessary formulas to derive the online E-step for an ED with a Gamma-distributed mixing variable.
PROPOSITION 5.1 (Online E-step with Gamma mixing variable).

$$ (12) \qquad \mathbb{E}[W | \mathbf{Y} = \mathbf{y}] = \frac{2\alpha + M}{u + 2\beta} $$

$$ (13) \qquad \mathbb{E}[\log(W) | \mathbf{Y} = \mathbf{y}] = \Psi^{(0)} \left( \frac{2\alpha + M}{2} \right) - \log \left( \frac{u + 2\beta}{2} \right), $$

with $u = (\mathbf{y} - \boldsymbol{\mu})^T \mathbf{D} \mathbf{A}^{-1} \mathbf{D}^T (\mathbf{y} - \boldsymbol{\mu})$ and $\Psi^{(0)}$ the diGamma function defined as the derivative of $\log \Gamma(\cdot)$.

PROOF. We have $(\mathbf{Y}|W=w) \sim \mathcal{N}_M(\boldsymbol{\mu}, \boldsymbol{\Sigma}/w)$ and $W \sim \mathcal{G}(\alpha, \beta)$, then

$$ \begin{aligned} p(w | \mathbf{Y} = \mathbf{y}) &\propto w^{\alpha + \frac{M}{2} - 1} \exp \left( -w \left( \frac{u}{2} + \beta \right) \right) \\ &\propto \mathcal{G} \left( w; \alpha + \frac{M}{2}, \frac{u}{2} + \beta \right), \end{aligned} $$

which provides the values for both required expectations. $\square$

<figure><img src="image_10.png" /></figure>

For the online maximization step, we need to find the solution of the root equation (26) in the main paper. By conveniently writing $s_w^T = [s_{1w} \ s_{2w}]$ and differentiating with respect to $\alpha$ and $\beta$, it comes the following system of equations,

$$ (14) \qquad s_{2w} - \Psi^{(0)}(\alpha) + \log \beta = 0 $$

$$ (15) \qquad -s_{1w} + \frac{\alpha}{\beta} = 0. $$

Then substituting $\alpha/s_{1w}$ to $\beta$ in (14) and under the condition that $s_{1w} \neq 0$, (14) leads to an equation with $\alpha$ only.

5.1. *Mixture of HD Student-t distributions.* As recalled in Table 1, the Student distribution is an ED with a Gamma-distributed mixing variable with $\alpha = \beta = \frac{\nu}{2}$ where $\nu$ is usually referred to as the degrees of freedom. Previous formulas apply and the root equation for $\nu$ simplifies into,

$$ s_{2w} - s_{1w} - \Psi^{(0)}\left(\frac{\nu}{2}\right) + \log\left(\frac{\nu}{2}\right) + 1 = 0. $$

<table><caption>TABLE 1
Some EDs with Gamma mixing variable.</caption><thead><tr><th>ED</th><th>Mixing variable</th></tr></thead><tbody><tr><td>Student</td><td>W ~ G(ν/2, ν/2)</td></tr><tr><td>Pearson Type VII</td><td>W ~ G(α - p/2, 1/2)</td></tr><tr><td>Generalized Student</td><td>W ~ G(ν/2, α/2)</td></tr><tr><td>Cauchy</td><td>W ~ G(1/2, 1/2)</td></tr></tbody></table>
**6. Inverse Exponential mixing variable and Laplace distribution.** Another example of widely used distribution is the symmetric multivariate Laplace distribution whose pdf in dimension M is

$$p(\mathbf{y}; \boldsymbol{\mu}, \boldsymbol{\Sigma}) = \frac{2}{(2\pi)^{M/2}(\det \boldsymbol{\Sigma})^{1/2}} \left(\sqrt{\frac{u}{2}}\right)^{1-\frac{M}{2}} K_{1-\frac{M}{2}}(\sqrt{2u})$$

where $\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$ are location and scale parameters, $u = (\mathbf{y} - \boldsymbol{\mu})^T \mathbf{D} \mathbf{A}^{-1} \mathbf{D}^T (\mathbf{y} - \boldsymbol{\mu})$ is the Mahalanobis distance and $K$ is the modified Bessel function of the second kind, see e.g. Appendix of Kotz, Kozubowski and Podgorski (2001). This is equivalent to the definition used by Eltoft, Kim and Lee (2006) by setting the parameter $\lambda$ therein to 1 and removing the constraint $\det \boldsymbol{\Sigma} = 1$, which for non fixed $\lambda$ is necessary to guarantee identifiability.

This multivariate distribution, can be written as a GSM with a mixing variable $W$ that follows an inverse exponential distribution, i.e. $W^{-1} = Z \sim \mathcal{E}(1)$. In our notation, the mixing variable $W$ is the inverse of the variable $Z$ used by Eltoft, Kim and Lee (2006). It follows that the quantity required for our online EM, namely $\mathbb{E}[W|\mathbf{Y} = \mathbf{y}]$ can be directly deduced from formula (20) in (Eltoft, Kim and Lee, 2006),

$$(16) \qquad \mathbb{E}[W|\mathbf{Y} = \mathbf{y}] = \sqrt{\frac{2}{u} \frac{K_{-\frac{M}{2}}(\sqrt{2u})}{K_{-\frac{M}{2}+1}(\sqrt{2u})}}.$$

We can then update the sufficient statistics following equation (15) in the main text. Applying Proposition 6.3, estimations for $\boldsymbol{\mu}$, $\mathbf{A}$ and $\tilde{\mathbf{D}}^*$ follow straightforwardly.

Note that the above moment can also be obtained by starting from the expression of the multivariate Laplace distribution as a GSM

$$\begin{aligned} p(\mathbf{y}; \boldsymbol{\mu}, \boldsymbol{\Sigma}) &= \int_{\mathbb{R}^+} \mathcal{N}_M (\mathbf{y}; \boldsymbol{\mu}; z\boldsymbol{\Sigma}) e^{-z} dz \\ &= \int_{\mathbb{R}^+} \mathcal{N}_M \left( \mathbf{y}; \boldsymbol{\mu}; \frac{\boldsymbol{\Sigma}}{w} \right) w^{-2} e^{-\frac{1}{w}} dw, \end{aligned}$$

and deriving that

$$p(w|\mathbf{Y} = \mathbf{y}) \propto w^{\frac{M}{2}-2} e^{-\frac{1}{2}wu - \frac{1}{w}} \propto \text{GIG}\left(w; \frac{M}{2} - 1, u, 2\right),$$

where GIG is the pdf of the Generalized Inverse Gaussian distribution, using the parametrization of Jørgensen (1982). The moments of the GIG are known and given, for $r \in \mathbb{R}$, by (Jørgensen, 1982),

$$\mathbb{E}[W^r|\mathbf{Y} = \mathbf{y}] = \left(\sqrt{\frac{2}{u}}\right)^r \frac{K_{\frac{M}{2}-1+r}(\sqrt{2u})}{K_{\frac{M}{2}-1}(\sqrt{2u})}.$$

For $r = 1$ it comes (16), using the fact that the modified Bessel functions satisfy $K_p = K_{-p}$. In practice we use Takekawa (2022) to make efficient computation of the ratio $\frac{K_{\cdot+1}(\cdot)}{K_{\cdot}(\cdot)}$ in Python-Jax.

**6.1. MRF reconstruction with High-Dimensional Mixture of Multivariate Laplace distributions (HD-LM).** We run the experiment to reconstruct the MRF signals using the same data and procedure, from Section 7 of the main paper. Training the model multiple times with different numbers of components and computing BIC gives an optimal number of classes of
$K = 30$. Figure 1 shows the BIC values and the different dimensions of the sub-spaces, with respect to the number of components. Table 2 shows the average MAE and SSIM over the 6 subjects, when compared to full matching. The results show clearly that the Laplace mixture model HD-LM performs much better than HD-GMM and DRONE but less than HD-STM.

<figure><img src="image_11.png" /><figcaption>FIG 1. Model selection. BIC scores (left) and clusters dimensions (right), with respect to the number of components $K$, for HD-LM applied to signals of dimension $M = 260$.</figcaption></figure>

<table><caption>TABLE 2
Parameter maps reconstruction with HD-LM. Average MAE, SSIM and standard deviations, over voxels and subjects with respect to full-matching.</caption><thead><tr><th></th><th><i>T</i><sub>1</sub> (ms)</th><th><i>T</i><sub>2</sub> (ms)</th><th>δ<i>f</i> (Hz)</th><th><i>B</i><sub>1</sub> sensitivity (10<sup>-3</sup>)</th><th>CBV (%)</th><th>R (μm)</th></tr></thead><tbody><tr><td>MAE</td><td>84 ± 18</td><td>12 ± 3</td><td>7 ± 1</td><td>45 ± 7</td><td>2.3 ± 1</td><td>1.9 ± 0, 2</td></tr><tr><td>SSIM</td><td>0.8 ± 0.02</td><td>0.76 ± 0.12</td><td>0.4 ± 0.05</td><td>0.53 ± 0.02</td><td>0.38 ± 0.07</td><td>0.25 ± 0.09</td></tr></tbody></table>

Finally, Figure 2 shows HD-LM reconstructions in comparison to full matching, for the different parameters and one slice of a subject. We observe some noise at the top of the $T_2$ map.

<figure><img src="image_12.png" /><figcaption>FIG 2. $T_1$, $T_2$, $\delta f$, $B_1$ sensitivity, CBV, R maps (columns) with full matching (top) and HD-LM with $K = 30$ (bottom)</figcaption></figure>

**7. Divide & Conquer high dimensional inverse regression of large data volumes.** As mentioned in the main paper conclusion, the presented data compression procedure could also be combined with a regression model in place of the matching step. Algorithm 1 below is a variant of Algorithm 1 in the paper, illustrating this possibility. Step 1 remains the same, Step 2 corresponds to the estimation of regression functions for each cluster and is
new, Step 3 is the cluster-wise prediction of new signals using the K estimated regression functions. More specifically, step 3.1 consists of projecting an observation $\tilde{\mathbf{y}}_j$ on each of the K subspaces and of computing the probability of $\tilde{\mathbf{y}}_j$ to be assigned to each cluster. Step 3.2 then corresponds to a cluster-wise prediction when only the regression function for the most probable cluster is used, i.e. $\tilde{\mathbf{t}}_j = \tilde{\mathbf{t}}_j^k$ with $k = \arg \max_l r_l(\tilde{\mathbf{y}}_j)$. Step 3.2 is written in a slightly more general manner using $\mathcal{L}$ to denote a function that would combine, in a way to be chosen, all predictions and all assignment probabilities.

**Algorithm 1** Divide & Conquer high dimensional inverse regression of large data volumes

**Input** Data set of simulated (signal, parameters) pairs $\mathcal{D}_f = \{\mathbf{y}_i, \mathbf{t}_i\}_{i=1:N}, N \gg 1, \mathbf{t}_i \in \mathbb{R}^L, \mathbf{y}_i \in \mathbb{R}^M, M \gg 1$. Observed signals $\{\tilde{\mathbf{y}}_j\}_{j=1:N}, \tilde{\mathbf{y}}_j \in \mathbb{R}^M$.

1. Reduced dimension representation of simulated data: $\{\hat{\mathbf{y}}_i, \mathbf{t}_i, r(\mathbf{y}_i)\}_{i=1:N}$.

1.1 Online HD-MED inference from $\{\mathbf{y}_i\}_{i=1:N}: K$ clusters, $d_k < M$ for $k = 1:K$ $\Rightarrow$ cluster assignment probabilities and cluster-wise projections $(\mathbf{r}, \mathbf{Q}) = \{r_k(\cdot), Q_k(\cdot)\}_{k=1:N}$.

1.2 Cluster-wise signal reduction: $\{\mathbf{y}_i\}_{i=1:N}, \mathbf{r}, \mathbf{Q} \Rightarrow \{\hat{\mathbf{y}}_i = Q_k(\mathbf{y}_i), i \in I_k\}$ with $I_k = \{i, \text{s.t } k = \arg \max_l r_l(\mathbf{y}_i)\}$, for $k = 1:K$.

2. Cluster-wise regressions: for $k = 1:K$, learn a regression function $l_k: \mathbb{R}^{d_k} \to \mathbb{R}^L$ from $\{(\mathbf{t}_i, \hat{\mathbf{y}}_i), i \in I_k\}$.

3. Cluster-wise predictions for observed signals:

3.1 Cluster-wise observed signal reduction: Use learned $(\mathbf{r}, \mathbf{Q})$ from step 1.1 to obtain $\{Q_k(\tilde{\mathbf{y}}_j), r_k(\tilde{\mathbf{y}}_j), j = 1:N, k = 1:K\}$

3.2 Cluster-wise prediction: For $j = 1:N, k = 1:K$, compute $\tilde{\mathbf{t}}_j^k = l_k(Q_k(\tilde{\mathbf{y}}_j)) \Rightarrow \text{predict } \tilde{\mathbf{t}}_j = \mathcal{L}(\tilde{\mathbf{t}}_j^1, r_1(\tilde{\mathbf{y}}_j), \dots, \tilde{\mathbf{t}}_j^K, r_K(\tilde{\mathbf{y}}_j)) \quad \text{e.g. } \tilde{\mathbf{t}}_j = \tilde{\mathbf{t}}_j^k \text{ with } k = \arg \max_l r_l(\tilde{\mathbf{y}}_j)$.

**Return** Predicted parameters $\{\tilde{\mathbf{t}}_j\}_{j=1:N}, \tilde{\mathbf{t}}_j \in \mathbb{R}^L$

**8. Application to MRF reconstruction.** We provide additional illustration by showing in Figure 3 the parameter maps obtained with full matching and HD-STM, for all the slices of the subject presented in main Figure 6. Both methods, HD-STM and full matching provide visually similar reconstructions.

**9. Generalized scale mixtures of Gaussian distributions.** In the main text, we focused on the case where the covariance matrix of each Gaussian component is scaled by a factor of $\frac{1}{w}$. However, as proposed in Lee and McLachlan (2014), it is possible to generalize Gaussian scale mixtures by using a more flexible link function $k(\cdot)$, leading to the formulation:

$$ (17) \qquad p(\mathbf{y}) = \int_{\mathbb{R}^+} \mathcal{N}_M(\mathbf{y}; \mu, k(w)\boldsymbol{\Sigma})f(w) dw. $$

This generalization encompasses a wider class of distributions, as summarized in Table 1 of Lee and McLachlan (2014). A possible continuation of our work is to derive properties from the main paper for this more general case, while still applying the proposed parcimonious parameterization of high dimensional scale matrices. It is noteworthy that we retain the ability to define a latent variable model for this extended class of distributions. Using the same notation, the following proposition holds for this new class, denoted by $\mathcal{HE}$ as in the main paper:
**PROPOSITION 9.1.** Let $d \le M-1$, $\mathbf{Y} \in \mathbb{R}^M$, $\mathbf{X} \in \mathbb{R}^d$, $\mathbf{E} \in \mathbb{R}^M$, $W \in \mathbb{R}^+$ be random variables, $\mathbf{V} \in \mathbb{R}^{M \times d}$ a matrix of linearly independent columns, $\mu \in \mathbb{R}^M$ a vector, a link function $k: \mathbb{R} \to \mathbb{R}$ and $f_\theta$ the pdf of a positive univariate random variable defined by some parameter $\theta$. Assume that

$$ 
\begin{align*}
\mathbf{Y} &= \mathbf{V}\mathbf{X} + \mu + \mathbf{E} \\
(\mathbf{X}|W=w) &\sim \mathcal{N}(\mathbf{0}_d, k(w)\mathbf{I}_d) \\
(\mathbf{E}|W=w) &\sim \mathcal{N}(\mathbf{0}_M, b k(w)\mathbf{I}_M) \\
W &\sim f_\theta
\end{align*}
$$

then,

$$ \mathbf{Y} \sim \mathcal{H}\mathcal{E}_{Md}(\mu, \tilde{\mathbf{D}}^*, \mathbf{a}, b, \theta), $$

with $\mathbf{I}_M + \mathbf{V}\mathbf{V}^T = \mathbf{D}\mathbf{A}\mathbf{D}^T$ and $\mathbf{A} = \text{diag}(a_1, \dots, a_d, b, \dots, b)$.

Additionally, denoting by $\mathbf{U} = b\mathbf{I}_d + \mathbf{V}^T\mathbf{V}$, we have,

$$ (18) \qquad (\mathbf{X}|\mathbf{Y}=\mathbf{y}, W=w) \sim \mathcal{N}(\mathbf{U}^{-1}\mathbf{V}^T(\mathbf{y}-\mu), k(w)b\mathbf{U}^{-1}). $$

It follows that $\mathbb{E}[\mathbf{Y}|\mathbf{X}=\mathbf{x}] = \mathbf{V}\mathbf{x} + \mu$ and $\mathbb{E}[\mathbf{X}|\mathbf{Y}=\mathbf{y}] = \mathbf{U}^{-1}\mathbf{V}^T(\mathbf{y}-\mu)$.

**PROOF.** The proof follows the same structure as Proposition 5.1 in the main paper, with the scaling factor $\frac{1}{w}$ replaced by the link function $k(w)$. $\square$

Regarding the OEM part, we can similarly derive propositions for this new class of distributions, assuming that $\mathbf{W}$ follows a distribution in the exponential family. For all link functions $k(w)$, most calculations involve replacing the scaling factor $\frac{1}{w}$ with $k(w)$. This opens the possibility of studying high-dimensional finite mixtures of distributions such as the slash, variance gamma, and Laplace distributions.

## REFERENCES

* AMIR, I., KOREN, T. and LIVNI, R. (2021). SGD Generalizes Better Than GD (And Regularization Doesn't Help). In *Proceedings of Thirty Fourth Conference on Learning Theory* (M. BELKIN and S. KLOTUFE, eds.). *Proceedings of Machine Learning Research* 134 63–92. PMLR.

* DIEULEVEUT, A., FORT, G., MOULINES, E. and WAI, H.-T. (2023). Stochastic Approximation Beyond Gradient for Signal Processing and Machine Learning. *IEEE Transactions on Signal Processing* 71 3117-3148.

* ELTOFT, T., KIM, T. and LEE, T.-W. (2006). On the multivariate Laplace distribution. *IEEE Signal Processing Letters* 13 300-303.

* GEPPERTH, A. and PFÜLB, B. (2021). Gradient-Based Training of Gaussian Mixture Models for High-Dimensional Streaming Data. *Neural Process. Lett.* 53 4331–4348.

* GILMAN, K., HONG, D., FESSLER, J. A. and BALZANO, L. (2023). Streaming Probabilistic PCA for Missing Data with Heteroscedastic Noise. *arXiv preprint arXiv:2310.06277*.

* HOSSEINI, R. and SRA, S. (2015). Matrix Manifold Optimization for Gaussian Mixtures. In *Advances in Neural Information Processing Systems* (C. CORTES, N. LAWRENCE, D. LEE, M. SUGIYAMA and R. GARNETT, eds.) 28. Curran Associates, Inc.

* JØRGENSEN, B. (1982). Statistical properties of the generalized inverse Gaussian distribution. *Lecture notes in statistics*; v.9. Springer-Verlag, New York.

* KOTZ, S., KOZUBOWSKI, T. and PODGORSKI, K. (2001). *The Laplace Distribution and Generalizations: A Review with Applications to Communications, Economics, Engineering, and Finance*. Progress in Mathematics. Birkhäuser Boston.

* LEE, S. X. and MCLACHLAN, G. J. (2014). Scale mixture distribution. *Wiley StatsRef: Statistics Reference Online* 1–16.

* NGUYEN, H. and FORBES, F. (2022). Global implicit function theorems and the online expectation-maximisation algorithm. *Australian & New Zealand Journal of Statistics* 64.
TAKEKAWA, T. (2022). Fast parallel calculation of modified Bessel function of the second kind and its derivatives. *SoftwareX* **17** 100923.

XU, L. and JORDAN, M. I. (1996). On Convergence Properties of the EM Algorithm for Gaussian Mixtures. *Neural Computation* **8** 129-151.

ZHANG, G., POUPART, P. and TRIMPONIAS, G. (2020). Comparing EM with GD in Mixture Models of Two Components. In *Proceedings of The 35th Uncertainty in Artificial Intelligence Conference 115* 164–174. PMLR.
<figure><img src="image_13.png" /><figcaption>FIG 3. T₁, T₂, δf, B₁ sensitivity, CBV, R maps (columns) with different methods- full matching (top), HD-STM with K = 25 (bottom) for all the slices of the subject presented in Figure 6 of the main paper.</figcaption></figure>