# Multiplicative Diffusion Models:
## Beyond Gaussian Latents

Robert Gruhlke,

Freie Universität at Berlin
Berlin, Germany

r.gruhlke@fu-berlin.de

Valentin Resseguier & Merveille Talla

UR OPAALE, ACTA Team, INRAE
Rennes, France

{valentin.resseguier,merveille.talla}@inrae.fr

## Abstract

We introduce a new class of generative models based on multiplicative score-driven diffusion. In contrast to classical diffusion models that rely on additive Gaussian noise, our construction is driven by skew-symmetric multiplicative noise. It yields conservative forward-backward dynamics inspired by the principles of physics. We prove that the forward process converges exponentially fast to a tractable non-Gaussian latent distribution, and we characterize this limit explicitly. A key property of our diffusion is that it preserves the distribution of data norms, resulting in a latent space that is inherently data-aware. Unlike the standard Gaussian prior, this structure better adapts to heavy-tailed and anisotropic data, providing a closer match between latent and observed distributions. On the algorithmic side, we derive the reverse-time stochastic differential equation and associated probability flow, and show that sliced score matching furnishes a consistent estimator for the backward dynamics. This estimation procedure is equivalent to maximizing an evidence lower bound (ELBO), bridging our framework with established variational principles. Empirically, we demonstrate the advantages of our model in challenging settings, including correlated Cauchy distributions and experimental fluid dynamics images ($d = 1024$). Across these tasks, our approach more accurately captures extreme events and tail behavior than classical diffusion models, particularly in the low-data regime. Our results suggest that multiplicative conservative diffusions open a principled alternative to current score-based generative models, with strong potential for domains where rare but critical events dominate.

## 1 Introduction

Mathematically equivalent (Song et al., 2021), diffusion models and score-based generative models demonstrate impressive skills and are among the current state-of-the-art for the generation of two- and three-dimensional images. Unconditioned sampling scores can be easily modified to conditioned sampling scores to address various inverse problems (Rybchuk et al., 2023; Rozet & Louppe, 2023; Daras et al., 2024; Bao et al., 2025). However, both learning and inference come with significant computational costs. In addition, even with large computational power, the generation of rare and extreme events remains a difficult task (Li et al., 2024; Stamatelopoulos & Sapsis, 2025). Those generative AI challenges may be more easily addressed by introducing physical-based inductive bias in the fully-data-driven approaches. In this paper, we take inspiration from physics and its conservative structure to build a multiplicative score-based generative model. It is inspired by transport noises in fluid dynamics (Kraichnan, 1968; Brzeźniak et al., 1991; Klyatskin, 1994; Piterbarg & Ostrovskii, 1997; Mikulevicius & Rozovskii, 2004; Mémin, 2014; Holm, 2015; Resseguier et al., 2021; Zhen et al., 2023) and, more generally, from slow-fast systems with multiplicative noise (Majda et al.,
1999; Franzke et al., 2005; Gottwald & Melbourne, 2013; Gottwald & Harlim, 2013). Transport noise models may be understood as generative models based on stochastic fluid dynamics rather than fitted neural networks. As other generative models, they suit particularly well to Bayesian inverse problems (Cotter et al., 2020b,a; Rességuier et al., 2022; Dufée et al., 2022).

Here, we might address problems outside the scope of fluid dynamics, though keeping the conservative structure of transport noise. The noising and denoising procedures that we propose maintain a part of the data information: the distribution of norm of the data point. The latent distribution is hence both tractable and close to the data distribution. More specifically, our contributions are the following.

**New generative model paradigm:** We introduce a new type of diffusion model where the noising process is multiplicative. We call it a Multiplicative Score-based Generative Model (MSGM). Involving random rotations around the origin, it greatly differs from previous diffusion models and opens a new research path. The key aspects are summarized in Figure 1.

**Deep theoretical analysis of MSGM:** Assuming a skew-symmetric structure and a rank condition for this noise, we proved several theoretical results, guiding the use of this new generative tool. The first theorem provides the Fokker-Planck equation of forward diffusion and its invariant measures. Then, we separately analyze the norm and direction of the diffusion. The norm is steady, whereas the direction follows a similar multiplicative stochastic differential equation (SDE). Two other theorems show that distributions of the direction and thus of the whole diffusion converge exponentially fast to a white noise in the weak sense. Asymptotically, the norm and direction are independent, and the latter is uniformly distributed over the $d$-sphere.

**General algorithm for MSGM:** We propose to estimate the scaled diffusion score by a neural network using sliced score matching, and our last theorem shows that it is equivalent to maximizing the ELBO. Sampling the non-Gaussian latent vectors reduces to a one-dimensional problem that we address with eCDF. For the denoising process, both SDE and ordinary differential equation (ODE) formulations are proposed.

**Application to extremes in moderate dimension:** We propose a numerical procedure to mimic the heavy-tail distribution with MSGM. We add a first layer to the neural network to perform a spherical decomposition with log-norm, and the latent distribution is characterized by the law of the data log-norm. Compared numerically with a standard diffusion model, MSGM better mimics multidimensional Cauchy distributions and measured fluid vorticity. The proximity between latent and data distributions facilitates the forward and the backward diffusions, and implicitly encompasses the correct distribution tail decays.

**Application in high dimension:** As a first step, we focus on MSGM scalability and design of sparse underlying tensors in the diffusion. While the latter is not covered completely by the theoretical analysis, our numerical experiments show promising image generation results.

<figure><img src="image_1.png" /><figcaption>Figure 1: Illustration of multiplicative score-based generative modeling (ours) compared to additive score-based generative modeling.</figcaption></figure>
# 2 Additive Score-based generative model

## 2.1 Forward and backward SDEs

Diffusion models or score-based generative models (SGM) can be expressed in continuous time with stochastic differential equation (SDE) (Song et al., 2021). The forward SDE is

$$d\vec{x}_s = -\vec{x}_s ds + \sqrt{2}d\vec{B}_s, \quad (2.1)$$

where $\vec{x}_s \in \mathbb{R}^d$ is distributed according to some density $p_s$ for $s > 0$, $s \mapsto \vec{B}_s$ is $d$-dimensional Brownian motion, and $\vec{x}_0$ distributed according to the dataset of interest. It is an Ornstein-Uhlenbeck process: the continuous-time version of a first-order autoregressive (AR) model and the distribution $p_s$ convergences to a standard Gaussian density exponentially for $s \to \infty$, e.g. in total variation or Wasserstein distance. We can then define for $t \in [0, T]$ the backward equation

$$d\vec{x}_t = \vec{x}_t dt + 2\nabla \log p_{T-t}(\vec{x}_t)dt + \sqrt{2}d\vec{B}_t, \quad (2.2)$$

with $t \mapsto \vec{B}_s$ another $d$-dimensional Brownian motion and $\vec{x}_0 \sim p_T$ (identifying the density $p_T$ with its distribution). Then for any $s \in [0, T]$, $\vec{x}_{T-s}$ and $\vec{x}_s$ have the same law $p_s$. In practice, when an approximate score $\nabla \log p_{T-t}$ is available we initialize equation 2.2 with a standard Gaussian distribution $\vec{x}_0 \sim \mathcal{N}(0, I_d)$ and integrate the backward SDE from $t=0$ to $t=T$ (i.e. from $s=T$ to $s=0$), ideally letting $\vec{x}_T$ become statistically similar to our dataset of interest.

## 2.2 A neural network to fit the score

In practice, the score $\nabla \log p_{T-t}(x)$ is approximated by a surrogate model, $s_\theta(x, T-t)$, e.g., a fitted artificial neural network (ANN). Alternatively, one can work on $a_\theta(x, T-t) = \sqrt{2}s_\theta(x, T-t)$ (Huang et al., 2021). For large-dimensional problems, Song et al. (2020) proposes to learn this neural network by *Sliced Score Matching* (SSM). Here, $a_\theta$ is obtained by minimizing the loss function

$$\mathcal{L}_{\text{SSM}}^{\text{SGM}}(\theta) = \int_{0}^{T} \mathbb{E}_{\vec{x}_s} \mathbb{E}_{v \sim \text{Rad}(d)} \left[ \frac{1}{2} \| a_{\theta}(\vec{x}_s, s) \|^2 + (v \cdot \nabla) \left( (\sqrt{2} a_{\theta}(\vec{x}_s, s) - \vec{x}_s) \cdot v \right) \right] ds. \quad (2.3)$$

where $\|\cdot\|$ is the Euclidean norm, Rad(d) denotes the $d$-dimensional Rademacher distribution and $\mathbb{E}_{\vec{x}_s}$ is the expectation along each path realization $\vec{x}_s$. Appendix A details the most common score matching losses and their link to the concept of the Evidence Lower Bound (ELBO).

# 3 Multiplicative Score-based generative model

Rather than relying on additive SDE equation 2.1, we propose a multiplicative SDE and the associated score-based generative model. Taking inspiration from physics, this approach introduces physical-based inductive bias and yields tractable latent distributions closer to the dataset distribution. In this section, we introduce our forward SDE based on skew-symmetric multiplicative noise, its corresponding latents, and backward SDE and analyze the limit properties of the process distribution. To share didactic similarities of the forward and backward processes as in the additive noise case, we will keep the same notation for the forward process $\vec{x}_s$ and the backward process $\vec{x}_s$, respectively.

## 3.1 Forward SDE

Instead of considering a forward SDE with additive noise, we rely on *multiplicative noise model*

$$d\vec{x}_s = \mathbf{G}(\vec{x}_s) \circ d\vec{B}_s, \quad (3.1)$$

where $d \ge 2$, $\mathbf{G}: \mathbb{R}^d \to \mathbb{R}^{d \times d}$ is linear and $\circ$ stands for the *Stratonovich* notation. The readers unfamiliar with this notation may interpret the Stratonovich noise $s \mapsto \circ d\vec{B}_s$ as a process with short correlation time but respecting the usual rules of differential calculus – says the chain rule. The discretized version of equation 3.1 – with an infinitely small time step $ds$ – may also provide insight:

$$\frac{1}{2}(\vec{x}_{s+ds} - \vec{x}_{s-ds}) = \mathbf{G}(\vec{x}_s) \frac{1}{2}(\vec{B}_{s+ds} - \vec{B}_{s-ds}). \quad (3.2)$$
For deeper understanding, Appendix B recalls some important notions of stochastic calculus, including the Stratonovich notation and the relationship to Itô calculus. Let $\mathbf{G}$ be represented by a third-order tensor $[\mathbf{G}_{i,j}^k] \in \mathbb{R}^{d,d,d}$ and define the random matrix $\mathbf{Z}_s = \sum_{k=1}^d \mathbf{G}^k(\vec{B}_s)_k$. Then, equation 3.1 can be written more explicitly as:

$$d\vec{\mathbf{x}}_s = \sum_{k=1}^{d} (\mathbf{G}^k \vec{\mathbf{x}}_s)(\circ d\vec{B}_s)_k = \sum_{k=1}^{d} \mathbf{G}^k (\circ d\vec{B}_s)_k \vec{\mathbf{x}}_s = \circ d\mathbf{Z}_s \vec{\mathbf{x}}_s \approx \frac{1}{2}(\mathbf{Z}_{s+ds} - \mathbf{Z}_{s-ds})\vec{\mathbf{x}}_s, \quad (3.3)$$

where the time increments of the random matrix, $\circ d\mathbf{Z}_s \approx \frac{1}{2}(\mathbf{Z}_{s+ds} - \mathbf{Z}_{s-ds})$, are uncorrelated.

Additionally, we will impose two assumptions valid throughout the paper:

Skew-symmetry: For any $k \in \{1, \dots, d\}$, the matrix $\mathbf{G}^k = (\mathbf{G}_{i,j}^k)_{i,j}$ is skew-symmetric. (A1)

Rank condition: For any $\mathbf{x} \in \mathbb{R}^d \setminus \{0\}$, $\mathrm{rank}(\mathbf{G}(\mathbf{x})) = d-1$. (A2)

In particular equation 3.1 does not describe a geometric Brownian motion, as the noise term $\mathbf{Z}_s$ is not diagonal. It includes zeros on the diagonal due to the skew-symmetry of $\mathbf{G}^k$. A geometric Brownian motion would necessitates $\mathbf{Z}_s = \mathrm{diag}(\vec{B}_s)$. A strategy to obtain a tensor $\mathbf{G}$ that matches assumptions A1 and A2 will be discussed in Section 6. By linearity, the skew-symmetry of all $\mathbf{G}^k$ (assumption A1) implies the skew-symmetry of the whole multiplicative noise matrix $\circ d\mathbf{Z}_s$. This structure is inspired by transport noises in fluid dynamics (Kraichnan, 1968; Piterbarg & Ostrovskii, 1997; Resseguier et al., 2021). In this analogy, $\vec{\mathbf{x}}_s$ would represent an image of temperature, advected by an incompressible fluid flow. Incompressibility leads to the skew-symmetry of the advection operator, and eventually to energy conservation of $\vec{\mathbf{x}}_s$. Here, we might address problems outside the scope of fluid dynamics, though maintaining the noise skew-symmetry (assumption A1) and thus the energy conservation, as discussed in Section 3.2.

With assumption A2, the noise spreads in a large space: $\mathrm{Im}(\mathbf{G}(\vec{\mathbf{x}}_s)) = \vec{\mathbf{x}}_s^\perp$. It ensures sufficient variability in the noising process and, in turn, a tractable distribution for $\vec{\mathbf{x}}_T$ when $T$ becomes large.

**Theorem 3.1.1.** Let the assumptions A1 and A2 hold. Then, the Fokker-Planck equation of equation 3.1 reads

$$\frac{\partial}{\partial s} p_s(\mathbf{x}) = \nabla_{\perp} \cdot \left( \frac{1}{2} \mathbf{\Sigma}(\mathbf{x}) \nabla_{\perp} p_s(\mathbf{x}) \right), \quad \mathbf{x} \in \mathbb{R}^d, \quad (3.4)$$

with conditional noise covariance $\mathbf{\Sigma}(\mathbf{x}) := \mathbf{G}(\mathbf{x})\mathbf{G}(\mathbf{x})^\top$ and $\nabla_{\perp}$ denoting the orthogonal projection of nabla on the tangent plane $\mathbf{x}^{\perp}$, i.e.

$$\nabla_{\perp} := (\mathbf{I}_d - \mathbf{x}^n(\mathbf{x}^n)^\top)\nabla, \quad (3.5)$$

for $\mathbf{x}^n := \mathbf{x}/\|\mathbf{x}\|$ with $\mathbf{x} \in D := \mathbb{R}^d \setminus \{0\}$ and 0 otherwise, the unit vector orthogonal to the d-sphere. Moreover, any stationary density $p_\infty$ of equation 3.4 is rotation-invariant on $\mathbb{R}^d$.

The proof is detailed in Appendix D.2. In order to highlight the connection to diffusion models on Riemannian manifolds, we note that $\nabla_{\perp}$ is the Riemannian gradient on the unit d-sphere: $\mathbb{S}^{d-1} = \{\mathbf{x} \in \mathbb{R}^d | \|\mathbf{x}\| = 1\}$, see Appendix H.1. For a possible extension of the considered diffusion equation 3.1 to the case of non-zero drift, we refer to section D.6.

## 3.2 Dynamics of norm and direction

We now consider for any $s \ge 0$ and $\vec{\mathbf{x}}_s \ne 0$, the spherical decomposition

$$\vec{\mathbf{x}}_s = \|\vec{\mathbf{x}}_s\| \vec{\mathbf{x}}_s^n \quad \text{with} \quad \vec{\mathbf{x}}_s^n := \vec{\mathbf{x}}_s / \|\vec{\mathbf{x}}_s\| \in \mathbb{S}^{d-1}. \quad (3.6)$$

First, we note that the norm, $\|\vec{\mathbf{x}}_s\|$, remains constant throughout the noising process. Indeed, the skew-symmetry of $\circ d\mathbf{Z}_s$ implies that $d\vec{\mathbf{x}}_s = \circ d\mathbf{Z}_s \vec{\mathbf{x}}_s$ is orthogonal to $\vec{\mathbf{x}}_s$ and hence:

$$d\|\vec{\mathbf{x}}_s\|^2 = 2\vec{\mathbf{x}}_s \cdot \circ d\vec{\mathbf{x}}_s = 0, \quad \forall s \ge 0. \quad (3.7)$$
Consequently, $||\vec{\tilde{x}}_s|| = ||\vec{\tilde{x}}_0||$. The vector $\vec{\tilde{x}}_s$ moves randomly on $||\vec{\tilde{x}}_0||\mathbb{S}^{d-1}$, the d-sphere of radius $||\vec{\tilde{x}}_0||$. Therefore, the distribution of the norms of the latent variable is exactly the distribution of the norms of the points of the dataset. We refer to Appendix D.3 for more details. This property will have important consequences for our learning procedure and on the possibility of generating extreme events. Indeed, the property of likely large norm events will be conserved along the diffusion. In particular, the norm distribution is one-dimensional and we can rely on advanced techniques available for this setup, without worrying about the curse of dimensionality. In practice, we will fit the distribution of the log-norm, $F_{\log|\cdot|, \epsilon}$, with eCDF. We refer to Appendix C for details and an overview of sampling from one-dimensional distributions.

We now focus on $p_s^n$, the distribution of the direction, $\vec{\tilde{x}}_s^n$, in particular as $s \to \infty$. For better readability, we postpone the full discussion and the main theorems to Appendix D.4. Lemma D.4.1 introduces the Fokker-Planck equation of the direction on $\mathbb{S}^{d-1}$ and its unique invariant measure, $p_\infty^n$. Then, Theorem D.4.2 shows the exponential convergence of the initial distribution $p_0^n$ to $p_\infty^n$, the uniform distribution on the unit sphere $\mathbb{S}^{d-1}$. Consequently, since $\mathbb{S}^{d-1}$ is compact, this implies convergence in total variation of $p_s^n$ to $p_\infty^n$ and convergence in distribution of $\vec{\tilde{x}}_s^n$ to $\vec{\tilde{x}}_\infty^n \sim U(\mathbb{S}^{d-1})$.

## 3.3 Non-Gaussian latent space

In this section, we characterize the generally non-Gaussian latent distribution. Although this sounds intractable at first glance, it will turn out that we can easily sample it.

In general, the latent space of MSGM is not Gaussian. It becomes Gaussian if and only if the distribution of the squared norms of the dataset points has $\chi^2$ distributions with $d$ degrees of freedom (see Appendix E.2). This property differs from the usual SGM. SGM latent variables are Gaussian, leading to $\chi^2$ distributions for the norms of the latent variables regardless of the data set. According to Lafon et al. (2023), without a heavy tail distribution for the latent variables, it is unlikely that the final samples will be generated with a heavy tail distribution, at least with variational autoencoders (VAE). With our approach, the distribution of the norms of the latent variables has heavy tails if and only if the distribution of the norms of dataset points has heavy tails. Therefore, we expect a significant improvement from our method in generating extreme events. In fact, for heavy-tailed data the KL divergence to the SGM latent distribution is infinite, whereas MSGM yields a finite value, see Appendix E.6. More generally, Appendix E.5 shows that the KL divergence from data to the latent distribution is always smaller under MSGM than SGM. So, only few time steps may be sufficient to integrate the forward and the backward MSGM diffusions. In any case, the MSGM latent vectors are still white noise in the weak sense (see Appendix E.1). Moreover, the norm and direction are independent from each other, which will drastically facilitate the sampling procedure. These results hold for latent vectors $\tilde{x}_\infty \sim p_\infty$. In practice, integrations of forward and backward diffusions are only possible over a finite time $T$. However, the following theorem states that the law of the solution, $\vec{\tilde{x}}_T$, will become close to $p_\infty$ exponentially as fast as $T \to +\infty$. So, we can confidently rely on finite-time integration.

**Theorem 3.3.1.** Let assumptions A1 and A2 hold. Let $\vec{\tilde{x}}_0 \sim p_0 \in C^2(D)$ and let $p_{\|\cdot|}$ be the (radial) density of $||\vec{\tilde{x}}_0||$. Then, the Fokker-Planck equation 3.4 has a unique solution $p_s \in C^2(D) \cap L^2(D)$ for all $s > 0$. Moreover, the Fokker-Planck equation has the stationary distribution

$$p_{\infty}(\tilde{x}) = \frac{p_{\|\cdot|}(\|\tilde{x}\|)}{\|\tilde{x}\|^{d-1} |\mathbb{S}^{d-1}|} \quad (3.8)$$

In particular, $||\vec{\tilde{x}}_s^n||$ and $\vec{\tilde{x}}_s^n$ are asymptotically independent for $s \to +\infty$. Moreover, there exists $\alpha = \alpha(\mathbf{G}, d) > 0$ such that

$$\|p_s - p_\infty\|_{L^2(\mathbb{R}^d)}^2 \le \exp(-\alpha s) \|p_0 - p_\infty\|_{L^2(\mathbb{R}^d)}^2 \quad (3.9)$$

The proof and details on $\alpha$ are given in Appendix D.5 and a specific case is discussed in Appendix J.3. The factor $\|\tilde{x}\|^{1-d}$ in equation 3.8 is expected. Indeed, $\frac{|\mathbb{S}^{d-1}|}{\|\tilde{x}\|^{1-d}}$ is the volume of the scaled d-sphere $\|\tilde{x}\|\mathbb{S}^{d-1}$, i.e. it corresponds to the uniform distribution on the scaled d-sphere $\|\tilde{x}\|\mathbb{S}^{d-1}$.

We will now consider the practical question on how to draw samples from the latent distribution with density $p_\infty$ from equation 3.8. It is of product structure between the radial and the directional
component. So, we can sample the norm $R_\infty$ and the direction $\vec{x}_\infty^n$ separately and multiply them. The norm $R_\infty$ can be sampled from an one-dimensional approximation of the data norm (see Appendix C) and the direction $\vec{x}_\infty^n$ is uniform. So, we can sample a $\vec{x}_\infty^N \sim N(0, I_d)$ and set

$$ \vec{x}_{\infty} = R_{\infty} \vec{x}_{\infty}^{n} \quad \text{with} \quad \vec{x}_{\infty}^{n} = \frac{\vec{x}_{\infty}^{N}}{\|\vec{x}_{\infty}^{N}\|}, \quad R_{\infty} = \hat{f}(\|\vec{x}_{\infty}^{N}\|), \quad (3.10) $$

$$ \text{and} \quad \hat{f}(r) := \exp \left( \hat{F}_{\log|\cdot|,\epsilon}^{-1} (F_{\chi^2(d)}(r^2)) \right) - \epsilon, \quad \forall r > 0. \quad (3.11) $$

If $\vec{x}_\infty^N = 0$, we set $\vec{x}_\infty = 0$. Proposition E.3.2 shows that this procedure leads to samples with the correct distribution, up to the approximation of the log-norm CDF $\hat{F}_{\log|\cdot|,\epsilon} \approx F_{\log|\cdot|,\epsilon}$. Moreover, the direct map $\hat{F}_{\log|\cdot|,\epsilon}$ can transform a latent vector, $\vec{x}_T$, into a Gaussian one (see Appendix E.4). This transformation may be useful for future applications like inverse problems or time evolution fittings.

## 3.4 Reverse ODE/SDE and score matching

From the Itô forward SDE (see Lemma D.1.2), we know that the Stratonovich reverse SDE writes

$$ d\vec{x}_t = \mathbf{G}(\vec{x}_t) \left( \mathbf{G}(\vec{x}_t)^{\top} \nabla \log p_{T-t}(\vec{x}_t) dt + \circ d\vec{B}_t \right), \quad (3.12) $$

and the reverse probability flow ODE is given as

$$ \frac{d\vec{x}_t}{dt} = \frac{1}{2} \mathbf{G}(\vec{x}_t) \left( \mathbf{G}(\vec{x}_t)^{\top} \nabla \log p_{T-t}(\vec{x}_t) \right). \quad (3.13) $$

The corresponding derivations are formulated in Proposition F.1 and Proposition F.2 and are proven in the appendix using Anderson (1982); Song et al. (2021). Following Huang et al. (2021), we directly model $\mathbf{G}(\vec{x}_t)^{\top} \nabla \log p_{T-t}(\vec{x}_t)$ by a neural network $a_\theta(\vec{x}_t, T-t)$. Additionally, we incorporate a spherical input layer, see Appendix L.4.1. We fit the parameters $\theta$ by sliced score matching (SSM) (Song et al., 2020), because in the multiplicative case we do not have an analytic formula for the conditional score $\nabla \log p_s(\vec{x}_s|\vec{x}_0)$ and because of the better scalability of SSM to high-dimensional problems that we would like to address in the future. To this end, we minimize the loss function:

$$ \mathcal{L}_{\text{SSM}}(\boldsymbol{\theta}) = \mathbb{E}_{s \sim \mathcal{U}[0,T]} \mathbb{E}_{\vec{x}_s} \mathbb{E}_{v \sim \text{Rad}(d)} \left[ \frac{1}{2} \|a_{\theta}(\vec{x}_s, s)\|^2 + (v \cdot \nabla)(\mathbf{G}(\vec{x}_s)a_{\theta}(\vec{x}_s, s)) \cdot v \right], \quad (3.14) $$

where Rad(d) denotes the Rademacher distribution in $\mathbb{R}^d$. The following theorem states that even in our multiplicative case, score matching is equivalent to maximize the ELBO, $\mathcal{E}_\infty$. In line with Benton et al. (2024); Ren et al. (2025), this theorem generalizes the result of Huang et al. (2021) and gives a theoretical justification for our score-matching loss equation 3.14. The derivation of this loss from the ELBO below is detailed in Appendix G.7.

**Theorem 3.4.1.** Let assumption A1 holds. Then, there exists a constant C such that

$$ p_0(\boldsymbol{x}) \ge \mathcal{E}_\infty(\boldsymbol{x}) := C - \int_0^T \mathbb{E}_{\vec{x}_s} \left[ \frac{1}{2} \|a_\theta(\vec{x}_s, s)\|^2 + \nabla \cdot (\mathbf{G}(\vec{x}_s) a_\theta(\vec{x}_s, s)) \Big|_{\vec{x}_0 = \boldsymbol{x}} \right] ds. \quad (3.15) $$

We proof this theorem in Appendix G. The first term $C := \mathbb{E}[\log p_T(\vec{x}_T)|\vec{x}_0 = \boldsymbol{x}]$ is a constant w.r.t. $\theta$. So, it has no effect on the optimization procedure. Therefore, even with our multiplicative noise, the minimization of the ELBO corresponds precisely to Implicit Score Matching (ISM), which is itself equivalent to explicit score matching (ESM), denoising score matching (DSM), and SSM (Huang et al., 2021). Note that formally replacing $\mathbf{G}$ by $\sqrt{2}$, we get the SGM SSM loss. For an easier numerical comparison in Section 6, we will also rely on SSM to train our baseline SGM.

# 4 Workflow

Algorithm 1 summarizes the proposed MSGM procedure. Here we make use of color to highlight the differences compared to SGM. For more details, we refer to Appendix L.
**Algorithm 1:** MSGM (Multiplicative Score-Based Generative Model).

**Input:** tensor **G**, one-dimensional distribution model $\hat{f}_{\gamma}$, data $\{\vec{x}_0^m\}_{m=1}^M$, $t_{\epsilon}$, time horizon $T$, time steps $N_T^f$ and $N_T^b$, time scheduler $g$, score model $a_{\theta}$, initial ANN parameter $\theta_0$, iterations $N_{\text{iter}}$
**— Training stage —**

1: $\gamma^* \leftarrow \text{fit_distribution}(\hat{f}_{\gamma}, \{\|\vec{x}_0^m\|\}_{m=1}^M)$ \quad \{Fitting of $\hat{f}_{\gamma}$, see Appendix C\}
2: **for** $n = 0$ to $N_{\text{iter}} - 1$ **do**
3: &nbsp;&nbsp;&nbsp;&nbsp; $\vec{x}_0 \sim \frac{1}{M} \sum_{m=1}^M \delta_{\vec{x}_0^m}$ \quad \{Random mini-batch from dataset\}
4: &nbsp;&nbsp;&nbsp;&nbsp; $s \sim U(I(t_{\epsilon}, T))$ \quad \{Sample uniform gridded time\}
5: &nbsp;&nbsp;&nbsp;&nbsp; $\vec{x}_s \leftarrow \text{SRK4}(s, \lfloor \frac{s}{T} N_T^f \rfloor, \vec{x}_0, 0, g\mathbf{G})$ \quad \{Forward diffusion integration via Algorithm 2\}
6: &nbsp;&nbsp;&nbsp;&nbsp; $\mathbf{v} \sim \text{Rad}(d)$ \quad \{Slicing directions\}
7: &nbsp;&nbsp;&nbsp;&nbsp; $\ell(\theta_n) \leftarrow \mathcal{L}_{\text{SSM}}^{\text{MSGM}}(s, \vec{x}_s, \mathbf{G}, g a_{\theta_n}, \mathbf{v})$ \quad \{Score-matching loss, from equation 3.14\}
8: &nbsp;&nbsp;&nbsp;&nbsp; $\theta_{n+1} \leftarrow \text{optimizer_update}(\theta_n, \ell(\theta_n))$ \quad \{e.g. via ADAM\}
9: **end for**
10: $\theta^* \leftarrow \theta_{N_{\text{iter}}}$ \quad \{Set final ANN parameter\}

**— Generative sampling stage —**

11: $\vec{x}_0 \stackrel{\leftarrow N}{\sim} \mathcal{N}(0, I_d)$, \quad \{Sample strong white noise\}
12: $\vec{x}_0 = \hat{f}_{\gamma^*} \left( \|\vec{x}_0^N\| \right) \frac{\vec{x}_0^N}{\|\vec{x}_0^N\|}$ \quad \{Sample weak white noise, see equation 3.11\}
13: $\vec{x}_T \leftarrow \text{SRK4}(T, N_T^b, \vec{x}_0, g\mathbf{G}a_{\theta^*}, g\mathbf{G})$ \quad \{Reverse diffusion integration via Algorithm 2\}
14: **return** $\gamma^*, \theta^*, \vec{x}_T$

**5 Related works**

Combining machine learning and mechanistic approaches is now a common approach. We may cite physics-informed neural networks (PINNs) (Raissi et al., 2019; Lu & Xu, 2024), physics-based prior covariance (Beauchamp et al., 2025; Clarotto et al., 2024), deep augmentation (Holzschuh et al., 2023; Fan et al., 2025), neural Galerkin (Lee & Carlberg, 2020; Chen et al., 2021; Romor et al., 2023; Finzi et al., 2023; Bruna et al., 2024; Kim et al., 2022), and chaos from energy-based models (Fournier & Pierfrancesco, 2025) among others. Here, we shall focus on score-based generative models. Bastek et al. (2024) add the physical equations inside their score matching loss. Holzschuh et al. (2023) fit a score to correct a backward physical equation but does not propose any generative model. To denoise corrupted images, several authors (e.g. Zhou et al., 2014; Shan et al., 2022; Guha & Acton, 2023) encode the multiplicative structure of speckle noise. Since this noise is not correlated between pixels, this approach strongly differs from ours. Most of these works do not deal with score or generative models. Guha & Acton (2023); Ren et al. (2025); Shetty et al. (2025) do, but their framework simplifies to SGM by considering the pixel-wise logarithm of images. Guth et al. (2022); Lempereur & Mallat (2024) encode a target multiscale structure (e.g. turbulence) by a hierarchy of normalized wavelets conditioned by the larger scales. Chen & Vanden-Eijnden (2025) adapt the noise to that multiscale structure in a stochastic interpolant context. Lobbe et al. (2023, 2025) replaced the Gaussian process involved in the transport-noise equations (Kraichnan, 1968; Piterbarg & Ostrovskii, 1997; Resseguier et al., 2021) by a Shrodinger bridge (De Bortoli et al., 2021). They inserted a SGM inside a transport noise dynamics, whereas we inserted a dynamics similar to transport noise inside a SGM. Following general Bayesian approaches, some of the literature on transport-noise relies on the Girsanov theorem to fit a drift modification or evaluate a likelihood (Cotter et al., 2020a; Singh et al., 2025)(see Appendix G.9). By extending the ELBO of (Huang et al., 2021) to SDEs inspired by transport noise, Theorem 3.4.1 justifies our fit of the backward SDE drift.

Several authors have recently proposed Langevin equations (Arnaudon et al., 2019; Luesink & Street, 2025; Ayala et al., 2025) and SGM (De Bortoli et al., 2022; Huang et al., 2022; Lou et al., 2023; Benton et al., 2024) on Riemannian manifolds in order to generate data lying on a particular manifold. Clearly different, our goal is more classical: generating data in $\mathbb{R}^d$. In our work, neither data nor their noisy versions are restricted to a single manifold. However, each solution path of our forward and backward SDE lies on its particular Riemannian manifold, the scaled $d$-sphere $\|\vec{x}_0\|\mathbb{S}^{d-1}$. De Bortoli et al. (2022) describes diffusions and SGM in the $d$-sphere $\mathbb{S}^{d-1}$. A detailed comparison is given in Appendix H.
Regarding extremes generation, using variational autoencoders (VAE), Lafon et al. (2023) argue that Gaussian latent restricts the generated samples to light-tail distributions. Accordingly, they propose to use fat-tailed latent distributions (see also Jaini et al. (2020); Huster et al. (2021) for normalizing flows (NF) and generative adversarial models (GAN) respectively). Yoon et al. (2023); Shariatian et al.; Pandey et al. (2024); Ren et al. (2025) proposed SGMs with ad hoc heavy-tailed (Lévy and Student-t) latent distributions. Our approach automatically makes the tails of the latent distribution fat when necessary. It learns it from the distribution of the data norm, $p_{|\cdot|}$. Similarly, the diffusion proposed by Dharmakeerthi et al. (2025) adapt to data but through a nonlinear drift and an additive noise. Li et al. (2024); Price et al. (2025); Stamatelopoulos & Sapsis (2025) and references therein show that the usual SGM may correctly represent extremes, especially in "interpolation mode", that is, when extremes lie on the interior of the dataset but have difficulties with extremes lying on the dataset boundaries. Our numerical experiments in Section 6 suggest that our method probably overcomes this limitation. To represent the directionality of extremes, many authors decompose norms and directions of extreme events (Engelke et al., 2019; Palacios-Rodríguez et al., 2020; Lafon et al., 2023; Naveau & Segers, 2024). Large-amplitude criterion (e.g. exceeding a high threshold) or fat-tail model can be applied on the norm. Extreme directions may or may not become asymptotically independent of their magnitude (Engelke et al., 2019; Lafon et al., 2023). Build on random rotations, MSGM naturally suggests such a polar decomposition. The extreme direction of the MSGM latent vector is asymptotically independent of its magnitude. However, the direction of the reverse process does depend on the magnitude (see Appendix H.2).

# 6 Experiments

For our numerical experiments, we choose to define a tensor $G^k$ in a simple way. We sample $d$ random matrices, keep only their skew-symmetric parts, and normalize:

$$G = \frac{\sqrt{d}}{\|G\|_2} \tilde{G} \quad \text{with} \quad \tilde{G}_{i,j}^k = \frac{1}{2}(M_{i,j}^k - M_{j,i}^k) \quad \text{and} \quad M_{i,j}^k \stackrel{iid}{\sim} \mathcal{N}(0, 1). \qquad (6.1)$$

In Appendix J we show that this random tensor $G$ respects conditions A1 and A2 almost surely. Appendix K proposes alternative tensor definitions with sparse structures that allow high-dimensional applications. Following Appendix K.2.2, MSGM can generate images as in Appendix M.7.2. For the test cases below, we also checked in Appendix M.6.2 and Appendix M.7.1 that the MSGM generation skills are equivalent with these sparse and dense tensors. However, these sparse tensors do not match the framework of Section 3.1 so we postpone the associated numerical evaluations to future works.

## 6.1 Multivariate Cauchy distribution

We first illustrate our method with a vector of Cauchy variables, $x_{Ca}$, with scale parameter $\gamma$:

$$(x_{Ca})_i \stackrel{iid}{\sim} p_{Ca} \quad \text{with} \quad p_{Ca}(x) := \frac{\gamma/\pi}{x^2 + \gamma^2}. \qquad (6.2)$$

It is an extreme case of fat-tailed distributions with a power-law tail: $p_{Ca}(x) \propto |x|^{-2}$ for large $|x|$. Real problems often involve both correlation and dimensionality $d > 2$. So, we correlate Cauchy variables, as $x_0 = A x_{Ca}$, with a fixed matrix, $A$, initialized with i.i.d. coefficients $A_{i,j} \sim \mathcal{N}(0, 1)$. Figure 2 confirms that, for $d = 4$, SGM hardly reproduces fat tails and extreme directionality, unlike MSGM. An explanation is the strong dissimilarity between the data distribution and the latent SGM distribution; see Appendix E.5 and Appendix E.6. For a larger number of ADAMS iterations, MSGM becomes more accurate, whereas SGM diverges (see Figure 4a). More plots, numerical comparisons, and experiments with variants of the state-of-the-art SGM can be found in Appendix M.6.1.

## 6.2 Measured vorticity fields

We also test our algorithm on fluid dynamics experimental data: small images of vorticity fields. These fields are two-dimensional curl of fluid velocity measured by Particle Image Velocimetry (PIV) in wind tunnels. Vorticity quantifies the local rotation speed of fluid and is known to have point-wise distributions with tails fatter than Gaussian ones (Wilczek & Friedrich, 2009). We focus on a benchmark fluid flow: a wake flow at Reynolds number 3900 created by a circular cylinder embedded in a mean stream (Parnaureau et al., 2008). Each vorticity sample is evaluated at $d = 16$ spatial
<figure><img src="image_2.png" /><figcaption>Figure 2: Pair plots of generated data (orange lines and dots) compared to ground truth data (blue lines and dots) with the SGM (left) and MSGM (right) for a vector of 4 correlated Cauchy variables. On the diagonal, log-histogram and logarithm of the pdf KDE estimation are superimposed.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 3: Pair plots of generated data (orange lines and dots) compared to ground truth data (blue lines and dots) with the SGM (left) and MSGM (right) trained on 1024 16-dimensional measured vorticity fields. On the diagonal, log-histogram and logarithm of the pdf KDE estimation are superimposed.</figcaption></figure>

points to ensure low dimensionality. We use limited training data (1024 data points) to make rare
events even more rare and learning more challenging. Appendix M.7 provides a deeper description of
this experimental dataset. Figure 3 highlights a larger concentration of points generated by SGM in
the center of the ground truth distribution. Accordingly, the tails of the marginals – i.e. the tails of
the vorticity point-wise distributions – are underestimated: SGM underestimates rare large vorticity
events. MSGM performs better since the MSGM latent distribution – easy to learn – is much closer
to the data distribution than SGM latent distribution, as theoretically suggested by Appendix E.5
and experimentally verified in Appendix M.7. In particular, the MSGM latent distribution seems
to have Laplace tails and to be more accurate in the low-data regime (see Figure 4). Additional we
carried out high-dimensional experiments with *d* = 1024 in Appendix M.7.2 based on sparse tensor
*G* developed in Appendix K. More details on data, preprocessing, illustrations, and other numerical
experiments are given in Appendix M.

# 7 Conclusion and discussion

We have proposed a new type of diffusion model with multiplicative noise. After a theoretical analysis of this ansatz, an algorithm is specified to mimic fat-tailed distributions, surpassing SGM in this task.
<figure><img src="image_4.png" /><figcaption>(a) Single MMD value for the correlated Cauchy distribution ($d = 4$) as a function of the number of effective ADAMS iterations.</figcaption></figure>

<caption>(b) Mean and 80% CI of MMD for the vorticity measurements distribution ($d = 16$) as a function of the number of training samples.</caption>

Figure 4: Convergence behaviors of the Maximum Mean Discrepancy (MMD) for two different test cases. $10^4$ samples are used for each MMD evaluation.

At this point, limitations of the general MSGM framework may be difficult to know. We rather discuss the limitations of the first numerical procedure applied in Section 6. First, the forward SDE has to be integrated numerically since we do not know an analytic solution for large-rank tensors (see Appendix I for the solution with low-rank tensors). It implies either a slower training or a reduced number of iterations compared to SGM. Moreover, we do not know of any analytic solution for the score in finite time. This prevents the use of DSM and force us to use ISM or SSM, a less stable approach. In the next future, we can hope that the active communities of generative models on symmetric Riemannian manifolds and, more generally, of stochastic differential geometry could come up with more efficient sampling algorithms and score evaluation procedures for our diffusions on $d$-spheres. In addition, random matrix theory and free probabilities (Biane, 1997; Delyon & Yao, 2006; Demni, 2008; Lévy, 2008; Delyon, 2010; Demni & Hmidi, 2012; Cébron, 2014) may provide alternative sampling methods and helpful results for large-dimension cases. Indeed, for some choices of $G$, the semigroup of our forward SDE may be expressed as a unitary Brownian matrix, converging for large dimensions to a free multiplicative Brownian matrix. Both theories could facilitate the sampling and the score evaluation of the MSGM forward diffusion. Moreover, a dense third-order tensor $G$ prevents image processing and other large-dimensional applications, related to, say, turbulent fluid dynamics. In fact, dimensions $d$ of such problems are very large – typically $d = O(10^5)$ or more. A dense tensor $G$ as we use in our numerical experiments has $d^3$ coefficients, and the memory and computational costs would become prohibitive in these cases. To address this issue, Appendix K proposes several sparse tensors and alternative to assumptions A1 and A2. Appendix M.7.2 shows first MSGM generated images in dimension $d = 1024$. Furthermore, we are currently developing physics-based sparse tensors $G$. Here, MSGM forward SDE is the spatial discretization of a stochastic partial differential equation involving transport noises (Kraichnan, 1968; Piterbarg & Ostrovskii, 1997; Resseguier et al., 2021). We expect that the physical inductive bias will facilitate both inference and learning, especially in low-data mode. Alternatively, the rank assumption A2 may be expressed more simply with the algebraic properties of $G$, eventually producing simple examples of sparse and efficient tensors.

In addition to the improvements discussed above, many paths remain to be explored. First, our theoretical results could be generalized to other multiplicative diffusions. We have considered dense linear maps $x \mapsto G(x)$ with $\mathrm{Im}(G(x)) = x^{\perp}$ for any $x \neq 0$. We believe that sparse linear maps of Appendix K and non-linear Lipchitz-continuous maps can yield similar theoretical results as long as that image condition is fulfilled for almost every $x \in \mathbb{R}^d$ (see Appendix K.1). The non-linear case would include in particular sphere-wise diffusions of De Bortoli et al. (2022) (see Appendix H.2). Second, we could address dynamical system forecasting. With the Gaussianization of MSGM latent vectors (see Appendix E.4) complex nonlinear dynamics could simplify to uncoupled one-dimensional linear dynamics as in Arbabi & Sapsis (2022). A third path to explore involves Bayesian inverse problems and data assimilation (Rozet & Louppe, 2023; Bao et al., 2025). Finally, our analytic solution issue could be bypassed by a normalizing flow approach: spherical decomposition, stochastic interpolants and flow along scaled $d$-spheres, taking inspiration from normalizing flow along Riemannian manifolds (e.g., Gemici et al., 2016; Mathieu & Nickel, 2020; Wu et al., 2025).
References

Ronald J Adrian and Jerry Westerweel. *Particle image velocimetry*. Number 30. Cambridge university press, 2011.

Brian DO Anderson. Reverse-time diffusion equation models. *Stochastic Processes and their Applications*, 12(3):313–326, 1982.

Hassan Arbabi and Themistoklis Sapsis. Generative stochastic modeling of strongly nonlinear flows with non-gaussian statistics. *SIAM/ASA Journal on Uncertainty Quantification*, 10(2):555–583, 2022.

Alexis Arnaudon, Alessandro Barp, and So Takao. Irreversible langevin mcmc on lie groups. In *Geometric Science of Information: 4th International Conference, GSI 2019, Toulouse, France, August 27–29, 2019, Proceedings 4*, pp. 171–179. Springer, 2019.

Ludwig Arnold. *Stochastic differential equations: theory and applications*. Wiley, 1974.

Mario Ayala, Nicolas Dirr, Grigorios A Pavliotis, and Johannes Zimmer. Reversibility, covariance and coarse-graining for langevin dynamics: On the choice of multiplicative noise. *arXiv preprint arXiv:2511.03347*, 2025.

Feng Bao, Hristo G Chipilski, Siming Liang, Guannan Zhang, and Jeffrey S Whitaker. Nonlinear ensemble filtering with diffusion models: Application to the surface quasigeostrophic dynamics. *Monthly Weather Review*, 153(7):1155–1169, 2025.

Jan-Hendrik Bastek, WaiChing Sun, and Dennis M Kochmann. Physics-informed diffusion models. *arXiv preprint arXiv:2403.14404*, 2024.

Maxime Beauchamp, Ronan Fablet, Simon Benaichouche, Pierre Tandeo, Nicolas Desassis, and Bertrand Chapron. Neural variational data assimilation with uncertainty quantification using spde priors. *Artificial Intelligence for the Earth Systems*, 4(3):240060, 2025.

Joe Benton, Yuyang Shi, Valentin De Bortoli, George Deligiannidis, and Arnaud Doucet. From denoising diffusions to denoising markov models. *Journal of the Royal Statistical Society Series B: Statistical Methodology*, 86(2):286–301, 2024.

Philippe Biane. Free brownian motion, free stochastic calculus and random matrices. *Free probability theory (Waterloo, ON, 1995)*, 12:1–19, 1997.

Joan Bruna, Benjamin Peherstorfer, and Eric Vanden-Eijnden. Neural galerkin schemes with active learning for high-dimensional evolution equations. *Journal of Computational Physics*, 496:112588, 2024.

Zdzislaw Brzeźniak, Marek Capiński, and Franco Flandoli. Stochastic partial differential equations and turbulence. *Mathematical Models and Methods in Applied Sciences*, 1(01):41–59, 1991.

Francesco Paolo Cantelli. Sulla determinazione empirica delle leggi di probabilita. *Giorn. Ist. Ital. Attuari*, 4(421-424), 1933.

Guillaume Cébron. *Processes on the unitary group and free probability*. PhD thesis, Université Pierre et Marie Curie-Paris VI, 2014.

Yanlai Chen, Sigal Gottlieb, Lijie Ji, and Yvon Maday. An eim-degradation free reduced basis method via over collocation and residual hyper reduction-based error estimation. *Journal of Computational Physics*, 444:110545, 2021.

Yifan Chen and Eric Vanden-Eijnden. Scale-adaptive generative flows for multiscale scientific data. *arXiv preprint arXiv:2509.02971*, 2025.

Lucia Clarotto, Denis Allard, Thomas Romary, and Nicolas Desassis. The spde approach for spatio-temporal datasets with advection and diffusion. *Spatial Statistics*, pp. 100847, 2024.

Colin Cotter, Dan Crisan, Darryl Holm, Wei Pan, and Igor Shevchenko. Data assimilation for a quasi-geostrophic model with circulation-preserving stochastic transport noise. *Journal of Statistical Physics*, 179(5):1186–1221, 2020a.
Colin Cotter, Dan Crisan, Darryl D Holm, Wei Pan, and Igor Shevchenko. A particle filter for stochastic advection by lie transport: a case study for the damped and forced incompressible two-dimensional euler equation. *SIAM/ASA Journal on Uncertainty Quantification*, 8(4):1446–1492, 2020b.

Colin J Cotter, Dan Crisan, and Maneesh Kumar Singh. Data assimilation for the stochastic camassa-holm equation using particle filtering: a numerical investigation. In *Stochastic Transport in Upper Ocean Dynamics Annual Workshop*, pp. 137–160. Springer Nature Switzerland Cham, 2023.

Giannis Daras, Hyungjin Chung, Chieh-Hsin Lai, Yuki Mitsufuji, Jong Chul Ye, Peyman Milanfar, Alexandros G Dimakis, and Mauricio Delbracio. A survey on diffusion models for inverse problems. *arXiv preprint arXiv:2410.00083*, 2024.

Valentin De Bortoli, James Thornton, Jeremy Heng, and Arnaud Doucet. Diffusion schrödinger bridge with applications to score-based generative modeling. *Advances in Neural Information Processing Systems*, 34:17695–17709, 2021.

Valentin De Bortoli, Emile Mathieu, Michael Hutchinson, James Thornton, Yee Whye Teh, and Arnaud Doucet. Riemannian score-based generative modelling. *Advances in Neural Information Processing Systems*, 35:2406–2422, 2022.

Bernard Delyon. Concentration inequalities for the spectral measure of random matrices. *Electronic Communications in Probability*, 15:549–562, 2010.

Bernard Delyon and Jian-Feng Yao. On the spectral distribution of gaussian random matrices. *Acta Mathematicae Applicatae Sinica*, 22(2):297–312, 2006.

Nizar Demni. Free jacobi process. *Journal of Theoretical Probability*, 21(1):118–143, 2008.

Nizar Demni and Taoufik Hmidi. Spectral distribution of the free unitary brownian motion: another approach. In *Séminaire de Probabilités XLIV*, pp. 191–206. Springer, 2012.

Prafulla Dhariwal and Alexander Nichol. Diffusion models beat gans on image synthesis. *Advances in neural information processing systems*, 34:8780–8794, 2021.

Kulunu Dharmakeerthi, Yousef El-Laham, Henry H Wong, Vamsi K Potluru, Changhong He, and Taosong He. Beyond linear diffusions: Improved representations for rare conditional generative modeling. *arXiv preprint arXiv:2510.02499*, 2025.

Benjamin Dufée, Etienne Mémin, and Dan Crisan. Stochastic parametrization: An alternative to inflation in ensemble kalman filters. *Quarterly Journal of the Royal Meteorological Society*, 148 (744):1075–1091, 2022.

Sebastian Engelke, Thomas Opitz, and Jennifer Wadsworth. Extremal dependence of random scale constructions. *Extremes*, 22(4):623–666, 2019.

Xiantao Fan, Deepak Akhare, and Jian-Xun Wang. Neural differentiable modeling with diffusion-based super-resolution for two-dimensional spatiotemporal turbulence. *Computer Methods in Applied Mechanics and Engineering*, 433:117478, 2025.

Marc Anton Finzi, Andres Potapczynski, Matthew Choptuik, and Andrew Gordon Wilson. A stable and scalable method for solving initial value pdes with neural networks. In *The Eleventh International Conference on Learning Representations*, 2023.

Samantha J. Fournier and Urbani Pierfrancesco. Generative modeling through internal high-dimensional chaotic activity. *Physical Review E*, 111.4:045304, 2025.

Christian Franzke, Andrew Majda, and Eric Vanden-Eijnden. Low-order stochastic mode reduction for a realistic barotropic model climate. *Journal of the atmospheric sciences*, 62(6):1722–1745, 2005.

Avner. Friedman. *Partial differential equations of parabolic type*. Prentice-Hall, Englewood Cliffs, N.J, 1964.
Mevlana C Gemici, Danilo Rezende, and Shakir Mohamed. Normalizing flows on riemannian manifolds. *arXiv preprint arXiv:1611.02304*, 2016.

Valery Glivenko. Sulla determinazione empirica delle leggi di probabilita. *Gion. Ist. Ital. Attauri.*, 4: 92-99, 1933.

Fabián González, O Deniz Akyildiz, Dan Crisan, and Joaquín Míguez. Nudging state-space models for bayesian filtering under misspecified dynamics: F. gonzález et al. *Statistics and Computing*, 35 (4):112, 2025.

Georg Gottwald and John Harlim. The role of additive and multiplicative noise in filtering complex dynamical systems. *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Science*, 469(2155):20130096, 2013.

Georg Gottwald and Ian Melbourne. Homogenization for deterministic maps and multiplicative noise. *Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences*, 469(2156), 2013.

Georg Gottwald, Daan Crommelin, and Christian Franzke. Stochastic climate theory. In *Nonlinear and Stochastic Climate Dynamics*. Cambridge University Press, 2015.

Arthur Gretton, Karsten M Borgwardt, Malte J Rasch, Bernhard Schölkopf, and Alexander Smola. A kernel two-sample test. *The journal of machine learning research*, 13(1):723-773, 2012.

Soumee Guha and Scott T Acton. Sddpm: Speckle denoising diffusion probabilistic models. *arXiv preprint arXiv:2311.10868*, 2023.

Florentin Guth, Simon Coste, Valentin De Bortoli, and Stephane Mallat. Wavelet score-based generative modeling. *Advances in neural information processing systems*, 35:478-491, 2022.

Darryl Holm. Variational principles for stochastic fluid dynamics. *Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences*, 471(2176), 2015.

Benjamin Holzschuh, Simona Vegetti, and Nils Thuerey. Solving inverse physics problems with score matching. *Advances in Neural Information Processing Systems*, 36, 2023.

Elton P Hsu. *Stochastic analysis on manifolds*. Number 38. American Mathematical Soc., 2002.

Chin-Wei Huang, Jae Hyun Lim, and Aaron C Courville. A variational perspective on diffusion-based generative models and score matching. *Advances in Neural Information Processing Systems*, 34: 22863-22876, 2021.

Chin-Wei Huang, Milad Aghajohari, Joey Bose, Prakash Panangaden, and Aaron C Courville. Riemannian diffusion models. *Advances in Neural Information Processing Systems*, 35:2750-2761, 2022.

Todd Huster, Jeremy Cohen, Zinan Lin, Kevin Chan, Charles Kamhoua, Nandi O Leslie, Cho-Yu Jason Chiang, and Vyas Sekar. Pareto gan: Extending the representational power of gans to heavy-tailed distributions. In *International Conference on Machine Learning*, pp. 4523-4532. PMLR, 2021.

Aapo Hyvärinen and Peter Dayan. Estimation of non-normalized statistical models by score matching. *Journal of Machine Learning Research*, 6(4), 2005.

Priyank Jaini, Ivan Kobyzev, Yaoliang Yu, and Marcus Brubaker. Tails of lipschitz triangular flows. In *International Conference on Machine Learning*, pp. 4673-4681. PMLR, 2020.

Youngkyu Kim, Youngsoo Choi, David Widemann, and Tarek Zohdi. A fast and accurate physics-informed neural network reduced order model with shallow masked autoencoder. *Journal of Computational Physics*, 451:110841, 2022.

Peter E Kloeden, Eckhard Platen, Peter E Kloeden, and Eckhard Platen. *Stochastic differential equations*. Springer, 1992.
Valerii I Klyatskin. Statistical description of the diffusion of a passive tracer in a random velocity field. *Physics-Uspekhi*, 37(5):501, 1994.

Robert Kraichnan. Small-scale structure of a scalar field convected by turbulence. *Physics of Fluids* (1958-1988), 11(5):945-953, 1968.

Hiroshi Kunita. *Stochastic flows and stochastic differential equations*, volume 24. Cambridge university press, 1997.

Nicolas Lafon, Philippe Naveau, and Ronan Fablet. A vae approach to sample multivariate extremes. *arXiv preprint arXiv:2306.10987*, 2023.

John M Lee. *Introduction to Riemannian manifolds*, volume 2. Springer, 2018.

Kookjin Lee and Kevin T Carlberg. Model reduction of dynamical systems on nonlinear manifolds using deep convolutional autoencoders. *Journal of Computational Physics*, 404:108973, 2020.

Etienne Lempereur and Stéphane Mallat. Hierarchic flows to estimate and sample high-dimensional probabilities. *arXiv preprint arXiv:2405.03468*, 2024.

Thierry Lévy. Schur-weyl duality and the heat kernel measure on the unitary group. *Advances in Mathematics*, 218(2):537-575, 2008.

Tianyi Li, Luca Biferale, Fabio Bonaccorso, Martino Andrea Scarpolini, and Michele Buzzicotti. Synthetic lagrangian turbulence by generative diffusion models. *Nature Machine Intelligence*, 6 (4):393-403, 2024.

Alexander Lobbe, Dan Crisan, and Oana Lang. Generative modelling of stochastic rotating shallow water noise. In *Stochastic Transport in Upper Ocean Dynamics Annual Workshop*, pp. 1-23. Springer Nature Switzerland Cham, 2023.

Alexander Lobbe, Dan Crisan, and Oana Lang. Bayesian inference for geophysical fluid dynamics using generative models. *Philosophical Transactions A*, 383(2299):20240321, 2025.

Aaron Lou, Minkai Xu, Adam Farris, and Stefano Ermon. Scaling riemannian diffusion models. *Advances in Neural Information Processing Systems*, 36:80291-80305, 2023.

Yulong Lu and Wuzhe Xu. Generative downscaling of pde solvers with physics-guided diffusion models. *Journal of Scientific Computing*, 101(3):1-23, 2024.

Erwin Luesink and Oliver D Street. Symplectic techniques for stochastic differential equations on reductive lie groups with applications to langevin diffusions. *arXiv preprint arXiv:2504.02707*, 2025.

Calvin Luo. Understanding diffusion models: A unified perspective. *arXiv preprint arXiv:2208.11970*, 2022.

Andrew Majda, Ilya Timofeyev, and Eric Vanden-Eijnden. Models for stochastic climate prediction. *Proceedings of the National Academy of Sciences*, 96(26):14687-14691, 1999.

Emile Mathieu and Maximilian Nickel. Riemannian continuous normalizing flows. *Advances in neural information processing systems*, 33:2503-2515, 2020.

Etienne Mémin. Fluid flow dynamics under location uncertainty. *Geophysical & Astrophysical Fluid Dynamics*, 108(2):119-146, 2014. doi: 10.1080/03091929.2013.836190.

Remigijus Mikulevicius and Boris Rozovskii. Stochastic Navier-Stokes equations for turbulent flows. *SIAM Journal on Mathematical Analysis*, 35(5):1250-1310, 2004.

Robb J Muirhead. *Aspects of multivariate statistical theory*. John Wiley & Sons, 2009.

Philippe Naveau and Johan Segers. Multivariate extreme value theory. *arXiv preprint arXiv:2412.18477*, 2024.

Bernt Oksendal. *Stochastic differential equations*. Spinger-Verlag, 1998.
Fátima Palacios-Rodríguez, Gwladys Toulemonde, Julie Carreau, and Thomas Opitz. Generalized pareto processes for simulating space-time extreme events: an application to precipitation reanalyses. *Stochastic Environmental Research and Risk Assessment*, 34:2033–2052, 2020.

Kushagra Pandey, Jaideep Pathak, Yilun Xu, Stephan Mandt, Michael Pritchard, Arash Vahdat, and Morteza Mardani. Heavy-tailed diffusion models. *arXiv preprint arXiv:2410.14171*, 2024.

Philippe Parnaureau, Johan Carlier, Dominique Heitz, and Eric Lamballais. Experimental and numerical studies of the flow over a circular cylinder at reynolds number 3900. *Physics of fluids*, 20(8), 2008.

Leonid Piterbarg and Alexander Ostrovskii. *Advection and Diffusion in Random Media: Implications for Sea Surface Temperature Anomalies*. Springer Science & Business Media, 1997.

Ilan Price, Alvaro Sanchez-Gonzalez, Ferran Alet, Tom R Andersson, Andrew El-Kadi, Dominic Masters, Timo Ewalds, Jacklynn Stott, Shakir Mohamed, Peter Battaglia, et al. Gencast: Diffusion-based ensemble forecasting for medium-range weather. In *105th Annual AMS Meeting 2025*, volume 105, pp. 449275, 2025.

Maziar Raissi, Paris Perdikaris, and George E Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. *Journal of Computational physics*, 378:686–707, 2019.

Yinuo Ren, Grant M Rotskoff, and Lexing Ying. A unified approach to analysis and design of denoising markov models. *arXiv preprint arXiv:2504.01938*, 2025.

Valentin Resseguier. Maximum likelihood estimation of subgrid flows from tracer image sequences. In *Stochastic Transport in Upper Ocean Dynamics Annual Workshop*, pp. 269–285. Springer Nature Switzerland Cham, 2023.

Valentin Resseguier, Long Li, Gabriel Jouan, Pierre Dérian, Etienne Mémin, and Bertrand Chapron. New trends in ensemble forecast strategy: uncertainty quantification for coarse-grid computational fluid dynamics. *Archives of Computational Methods in Engineering*, 28:215–261, 2021.

Valentin Resseguier, Matheus Ladvig, and Dominique Heitz. Real-time estimation and prediction of unsteady flows using reduced-order models coupled with few measurements. *Journal of Computational Physics*, 471:111631, 2022.

Francesco Romor, Giovanni Stabile, and Gianluigi Rozza. Non-linear manifold reduced-order models with convolutional autoencoders and reduced over-collocation method. *Journal of Scientific Computing*, 94(3):74, 2023.

François Rozet and Gilles Louppe. Score-based data assimilation. *Advances in Neural Information Processing Systems*, 36:40521–40541, 2023.

Alex Rybchuk, Malik Hassanaly, Nicholas Hamilton, Paula Doubrawa, Mitchell J Fulton, and Luis A Martínez-Tossas. Ensemble flow reconstruction in the atmospheric boundary layer from spatially limited measurements through latent diffusion models. *Physics of Fluids*, 35(12), 2023.

Xiujie Shan, Jiebao Sun, Zhichang Guo, Wenjuan Yao, and Zhenyu Zhou. Fractional-order diffusion model for multiplicative noise removal in texture-rich images and its fast explicit diffusion solving. *BIT Numerical Mathematics*, 62(4):1319–1354, 2022.

Dario Shariatian, Umut Simsekli, and Alain Oliviero Durmus. Heavy-tailed diffusion with denoising levy probabilistic models. In *The Thirteenth International Conference on Learning Representations*.

Nishanth Shetty, Madhava Prasath, and Chandra Sekhar Seelamantula. Dale meets langevin: A multiplicative denoising diffusion model. *arXiv preprint arXiv:2510.02730*, 2025.

Maneesh Kumar Singh, Joshua Hope-Collins, Colin J Cotter, and Dan Crisan. Data assimilation using a global girsanov nudged particle filter. *arXiv preprint arXiv:2507.17685*, 2025.

Yang Song, Sahaj Garg, Jiaxin Shi, and Stefano Ermon. Sliced score matching: A scalable approach to density and score estimation. In *Uncertainty in Artificial Intelligence*, pp. 574–584. PMLR, 2020.
Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole. Score-based generative modeling through stochastic differential equations. In *International Conference on Learning Representations*, 2021.

Stamatis Stamatelopoulos and Themistoklis P Sapsis. Can diffusion models capture extreme event statistics? *Computer Methods in Applied Mechanics and Engineering*, 435:117589, 2025.

Michael E. Taylor. *Partial Differential Equations III: Nonlinear Parabolique Equations*, volume 117 of *Applied Mathematical Sciences*. Springer-Verlag, New York, 2nd edition, 2011. ISBN 978-1-4419-7048-0.

Surya T Tokdar, Sheng Jiang, and Erika L Cunningham. Heavy-tailed density estimation. *Journal of the American Statistical Association*, 119(545):163–175, 2024.

Howard G Tucker. A generalization of the glivenko-cantelli theorem. *The Annals of Mathematical Statistics*, 30(3):828–830, 1959.

Michael Wilczek and Rudolf Friedrich. Dynamical origins for non-gaussian vorticity distributions in turbulent flows. *Physical Review E—Statistical, Nonlinear, and Soft Matter Physics*, 80(1):016316, 2009.

Jiawen Wu, Bingguang Chen, Yuyi Zhou, Qi Meng, Rongchan Zhu, and Zhi-Ming Ma. Riemannian neural geodesic interpolant. *arXiv preprint arXiv:2504.15736*, 2025.

Eun Bi Yoon, Keehun Park, Sungwoong Kim, and Sungbin Lim. Score-based generative models with lévy processes. *Advances in Neural Information Processing Systems*, 36:40694–40707, 2023.

Yicun Zhen, Valentin Resseguier, and Bertrand Chapron. Physically constrained covariance inflation from location uncertainty. *Nonlinear Processes in Geophysics*, 30(2):237–251, 2023.

Zhenyu Zhou, Zhichang Guo, Gang Dong, Jiebao Sun, Dazhi Zhang, and Boying Wu. A doubly degenerate diffusion model based on the gray level indicator for multiplicative noise removal. *IEEE Transactions on Image Processing*, 24(1):249–260, 2014.

# Appendix

## Table of Contents

<table><tr><td>A Losses for score matching</td><td>18</td></tr><tr><td>B Stochastic Calculus and Stratonovich Integrals</td><td>18</td></tr><tr><td>B.1 Itô Integrals and SDEs</td><td>18</td></tr><tr><td>B.2 Stratonovich Integrals and Chain rule</td><td>19</td></tr><tr><td>B.3 Conversion between Itô and Stratonovich forms</td><td>19</td></tr><tr><td>B.4 Fokker-Planck equation</td><td>19</td></tr><tr><td>C Sampling from 1D distributions</td><td>20</td></tr><tr><td>D The Fokker-Planck equation and its invariant measures</td><td>20</td></tr><tr><td>D.1 Itô form of the forward SDE</td><td>20</td></tr><tr><td>D.2 Fokker-Planck equation and Theorem 3.1.1</td><td>21</td></tr><tr><td>D.3 Distribution of the norms</td><td>23</td></tr><tr><td>D.4 Fokker-Planck equation of the direction</td><td>25</td></tr><tr><td>D.5 Proof of Theorem 3.3.1: Convergence of Fokker-Planck equation</td><td>27</td></tr></table>
<table>
  <tbody>
    <tr>
      <td>D.6</td>
      <td>Beyond pure Stratonovich noise</td>
      <td>31</td>
    </tr>
    <tr>
      <td>E</td>
      <td>Latent distribution</td>
      <td>31</td>
    </tr>
    <tr>
      <td>E.1</td>
      <td>The invariant measures define white noises in the weak sense</td>
      <td>31</td>
    </tr>
    <tr>
      <td>E.2</td>
      <td>Condition of Gaussianity for the latent vector</td>
      <td>32</td>
    </tr>
    <tr>
      <td>E.3</td>
      <td>A tractable algorithm to sample latent vectors</td>
      <td>32</td>
    </tr>
    <tr>
      <td>E.4</td>
      <td>Gaussianization of the latent vectors</td>
      <td>33</td>
    </tr>
    <tr>
      <td>E.5</td>
      <td>A shorter distance between latent and data distribution</td>
      <td>33</td>
    </tr>
    <tr>
      <td>E.6</td>
      <td>Relevance of MSGM latent space for heavy-tail distributions.</td>
      <td>34</td>
    </tr>
    <tr>
      <td>F</td>
      <td>Backward diffusion</td>
      <td>37</td>
    </tr>
    <tr>
      <td>G</td>
      <td>Proof of Theorem 3.4.1: equivalence between ELBO and score matching</td>
      <td>38</td>
    </tr>
    <tr>
      <td>G.1</td>
      <td>Statement of the theorem</td>
      <td>38</td>
    </tr>
    <tr>
      <td>G.2</td>
      <td>Notations correspondence</td>
      <td>39</td>
    </tr>
    <tr>
      <td>G.3</td>
      <td>Marginal density from Feynman-Kac representation</td>
      <td>39</td>
    </tr>
    <tr>
      <td>G.4</td>
      <td>Change of measure and Jensen's inequality</td>
      <td>40</td>
    </tr>
    <tr>
      <td>G.5</td>
      <td>Girsanov theorem</td>
      <td>40</td>
    </tr>
    <tr>
      <td>G.6</td>
      <td>ELBO evaluation</td>
      <td>41</td>
    </tr>
    <tr>
      <td>G.7</td>
      <td>From ELBO to our SSM loss</td>
      <td>41</td>
    </tr>
    <tr>
      <td>G.8</td>
      <td>Remark on the score parametrization</td>
      <td>42</td>
    </tr>
    <tr>
      <td>G.9</td>
      <td>Girsanov theorem in the transport noise literature</td>
      <td>43</td>
    </tr>
    <tr>
      <td>H</td>
      <td>Comparison with diffusions on Riemannian manifolds</td>
      <td>43</td>
    </tr>
    <tr>
      <td>H.1</td>
      <td>Riemannian Manifolds and Differentiation</td>
      <td>43</td>
    </tr>
    <tr>
      <td>H.2</td>
      <td>Conditional diffusions on scaled d-spheres</td>
      <td>44</td>
    </tr>
    <tr>
      <td>H.3</td>
      <td>Link with neural network architecture</td>
      <td>46</td>
    </tr>
    <tr>
      <td>I</td>
      <td>Analytic illustrations on simplified cases</td>
      <td>46</td>
    </tr>
    <tr>
      <td>I.1</td>
      <td>The two-dimensional case</td>
      <td>46</td>
    </tr>
    <tr>
      <td>I.2</td>
      <td>Tensor built from a single skew-symmetric matrix</td>
      <td>47</td>
    </tr>
    <tr>
      <td>I.3</td>
      <td>Non-commutativity in the general case</td>
      <td>48</td>
    </tr>
    <tr>
      <td>J</td>
      <td>Rank and skew-symmetry conditions for random tensor G</td>
      <td>48</td>
    </tr>
    <tr>
      <td>J.1</td>
      <td>Proof of the rank condition</td>
      <td>49</td>
    </tr>
    <tr>
      <td>J.2</td>
      <td>Tensor renormalization</td>
      <td>52</td>
    </tr>
    <tr>
      <td>J.3</td>
      <td>Mean speed of convergence with renormalized tensor.</td>
      <td>52</td>
    </tr>
    <tr>
      <td>K</td>
      <td>Going beyond the rank condition for MSGM scalability</td>
      <td>53</td>
    </tr>
    <tr>
      <td>K.1</td>
      <td>Weaker assumptions</td>
      <td>53</td>
    </tr>
    <tr>
      <td>K.2</td>
      <td>Sparse tensors</td>
      <td>55</td>
    </tr>
    <tr>
      <td>K.3</td>
      <td>Discussion about local and non local structure</td>
      <td>58</td>
    </tr>
    <tr>
      <td>L</td>
      <td>Numerical scheme</td>
      <td>59</td>
    </tr>
    <tr>
      <td>L.1</td>
      <td>Numerical integration of SDEs</td>
      <td>59</td>
    </tr>
    <tr>
      <td>L.2</td>
      <td>Scheduling</td>
      <td>59</td>
    </tr>
    <tr>
      <td>L.3</td>
      <td>Loss evaluation</td>
      <td>61</td>
    </tr>
    <tr>
      <td>L.4</td>
      <td>Neural network architecture</td>
      <td>62</td>
    </tr>
    <tr>
      <td>M</td>
      <td>Details about our numerical experiments</td>
      <td>64</td>
    </tr>
    <tr>
      <td>M.1</td>
      <td>Test cases</td>
      <td>64</td>
    </tr>
    <tr>
      <td>M.2</td>
      <td>Data preprocessing</td>
      <td>64</td>
    </tr>
    <tr>
      <td>M.3</td>
      <td>Comparison strategy</td>
      <td>64</td>
    </tr>
    <tr>
      <td>M.4</td>
      <td>Swiss roll</td>
      <td>65</td>
    </tr>
    <tr>
      <td>M.5</td>
      <td>Anisotropic Gaussian distribution</td>
      <td>66</td>
    </tr>
    <tr>
      <td>M.6</td>
      <td>Multivariate Cauchy distribution</td>
      <td>66</td>
    </tr>
  </tbody>
</table>
$\mathbb{R}^{d \times m}$ satisfying appropriate integrability conditions, the *Itô integral* of $X$ with respect to $B$ is defined as the mean-square limit:

$$ \int_{0}^{T} X_s dB_s := \lim_{|\Pi| \to 0} \sum_{[t_i, t_{i+1}] \in \Pi} X_{t_i} (B_{t_{i+1}} - B_{t_i}), \quad (B.1) $$

where the sum is taken over a partition $\Pi$ of $[0, T]$.

An SDE interpreted in the *Itô* sense reads:

$$ dX_t = f(X_t, t) dt + G(X_t, t) dB_t, \quad (B.2) $$

where $f: \mathbb{R}^d \times \mathbb{R}_+ \to \mathbb{R}^d$ is the drift, and $G: \mathbb{R}^d \times \mathbb{R}_+ \to \mathbb{R}^{d \times m}$ is the diffusion coefficient.

## B.2 Stratonovich Integrals and Chain rule

Unlike the *Itô* integral, the *Stratonovich integral* is defined using a symmetric discretization:

$$ \int_{0}^{T} X_s \circ dB_s := \lim_{|\Pi| \to 0} \sum_{[t_i, t_{i+1}] \in \Pi} \frac{X_{t_i} + X_{t_{i+1}}}{2} (B_{t_{i+1}} - B_{t_i}). \quad (B.3) $$

A Stratonovich SDE is written as:

$$ dX_t = f_S(X_t, t) dt + G(X_t, t) \circ dB_t. \quad (B.4) $$

A key advantage of the Stratonovich formulation is that it satisfies the classical chain rule. For any smooth function $\phi: \mathbb{R}^d \to \mathbb{R}$, we have:

$$ d\phi(X_t) = \nabla\phi(X_t)^\top f_S(X_t, t) dt + \nabla\phi(X_t)^\top G(X_t, t) \circ dB_t. \quad (B.5) $$

Moreover, in multiscale deterministic or stochastic equations, if the fast component is a continuous process with infinitesimal correlation time, the slow component generally converges to the solution of another SDE. In this other SDE, the fast component is often replaced by a Stratonovich integral (Arnold, 1974). Note that it is not always true for nonlinear dynamics (Gottwald & Melbourne, 2013; Gottwald et al., 2015). Accordingly, the readers may interpret the Stratonovich noise $s \mapsto \circ dB_s$ as a formal representation of a process with a short correlation time that nevertheless respects the classical rules of differential calculus, in particular, the chain rule.

## B.3 Conversion between Itô and Stratonovich forms

Given a Stratonovich SDE, it is always possible to convert it to the equivalent *Itô* form:

$$ dX_t = \left( f_S(X_t, t) + \frac{1}{2} \sum_{j=1}^{m} G_j(X_t, t) \cdot \nabla G_j(X_t, t) \right) dt + G(X_t, t) dB_t, \quad (B.6) $$

where $G_j(x, t)$ is the j-th column of the diffusion matrix $G(x, t)$. The additional drift term arises from the correction due to the non-zero quadratic variation of the noise.

## B.4 Fokker-Planck equation

An *Itô* SDE of the form

$$ dX_t = f(X_t, t) dt + G(X_t, t) dB_t, \quad (B.7) $$

induces a time-evolution equation for the probability density $p(x, t)$ of $X_t$. This is known as the *Fokker-Planck equation*, given by:

$$ \partial_t p(x, t) = -\nabla \cdot (f(x, t) p(x, t)) + \frac{1}{2} \nabla \cdot (\nabla \cdot (\Sigma(x, t) p(x, t)))^\top, \quad (B.8) $$

where $\Sigma(x, t) := G(x, t)G(x, t)^\top$ is the diffusion tensor. The Fokker-Planck equation describes the deterministic evolution of the probability density associated with the stochastic process.
C Sampling from 1D distributions

Let us denote by $p_{|\cdot|}$ the distribution of the norms, $\|\vec{\mathbf{x}}_T\|$. In MSGM, it is also the distribution of $\|\vec{\mathbf{x}}_0\|$ (see Proposition D.3.2). This distribution is arbitrary, but is a one-dimensional distribution. So, it is straightforward to learn and sample from, e.g., using an empirical cumulative distribution function (eCDF) (Cantelli, 1933; Glivenko, 1933; Tucker, 1959). Norms are positive and might be close to zero, so in practice we work with a regularized log-norm: $\log \|\mathbf{x}\|_{\epsilon} := \log(\|\mathbf{x}\| + \epsilon)$ with $\epsilon$ small, typically $\epsilon = 10^{-6}$. From a data set of the log-norms of $M$ training samples, $(\log \|\vec{\mathbf{x}}_T^{(i)}\|_{\epsilon})_i = (\log \|\vec{\mathbf{x}}_0^{(i)}\|_{\epsilon})_i$, we define eCDF $\hat{F}_{\log|\cdot|_{\epsilon}}$ as

$$ \hat{F}_{\log|\cdot|_{\epsilon}} : \mathbb{R} \to [0, 1], \qquad (C.1) $$

$$ R \mapsto \frac{1}{M} \sum_{i=1}^{M} \mathbf{1}_{\{R \ge \log \|\vec{\mathbf{x}}_T^{(i)}\|_{\epsilon}\}}. \qquad (C.2) $$

Then, we approximate the distribution of the norms, $p_{\log|\cdot|_{\epsilon}}(R)dR$, by the empirical one, $\hat{p}_{\log|\cdot|_{\epsilon}}(R)dR := \hat{F}_{\log|\cdot|_{\epsilon}}(dR)$. In particular, we can sample a new norm of latent variables, $\|\vec{\mathbf{x}}_T\|$, from a uniform one-dimensional variable $u \sim \mathcal{U}(0, 1)$ as follows

$$ \|\vec{\mathbf{x}}_T\| = \|\vec{\mathbf{x}}_0\| = \exp(\log \|\vec{\mathbf{x}}_0\|_{\epsilon}) - \epsilon = \exp(\hat{F}_{\log|\cdot|_{\epsilon}}^{-1}(u)) - \epsilon. \qquad (C.3) $$

eCDF is an efficient tool, but it cannot generalize the distribution $\hat{p}_{|\cdot|}$ outside the training set $(\|\vec{\mathbf{x}}_T^{(i)}\|)_i$. For better generalization, instead, one could use a one-dimensional kernel density estimation or fitting of parametric distributions. In the case of one-dimensional distributions with fat tails, classical kernel density estimation (KDE) suffers from bias in the tail estimation or peaks due to sparse data in the tails. In that case, one could consider more robust approaches that, in general, do not require the existence of moments of the target distribution Tokdar et al. (2024).

In this paper, we rely on the eCDF.

D The Fokker-Planck equation and its invariant measures

D.1 Itô form of the forward SDE

Define the conditional noise covariance $\Sigma(\mathbf{x})$ as

$$ \Sigma(\mathbf{x}) := \mathbf{G}(\mathbf{x})\mathbf{G}(\mathbf{x})^{\top} = \mathbb{E}[(d\vec{\mathbf{x}}_s)(d\vec{\mathbf{x}}_s)^{\top}|\vec{\mathbf{x}}_s = \mathbf{x}]. \qquad (D.1) $$

We begin with a lemma.

**Lemma D.1.1.** Let the skew-symmetry assumption A1 hold. Then,

$$ \frac{1}{2} \nabla \cdot (\Sigma(\mathbf{x})) = \frac{1}{2} \sum_{k=1}^{d} (G^k)^2 \mathbf{x}. \qquad (D.2) $$

*Proof.* Let us explicitly state the matrix divergence. For $k = 1, \dots, d$ define $\Sigma^k(\mathbf{x}) = [\Sigma_{ij}^k(\mathbf{x})] := G^k \mathbf{x} \mathbf{x}^\top (G^k)^\top$, then we decompose $\Sigma$ as

$$ \Sigma(\mathbf{x}) := \mathbf{G}(\mathbf{x})\mathbf{G}^{\top}(\mathbf{x}) = \sum_{k=1}^{d} \Sigma^k(\mathbf{x}). \qquad (D.3) $$
Then, taking the divergence of the j-th column of $\Sigma^k(\mathbf{x})$, we obtain

$$ (\nabla \cdot (\Sigma^k(\mathbf{x})))_j = \nabla \cdot [\Sigma^k(\mathbf{x})]_{:,j} = \sum_{i=1}^{d} \frac{\partial}{\partial x_i} \Sigma_{ij}^k(\mathbf{x}), \qquad (D.4) $$

$$ = \sum_{i=1}^{d} \frac{\partial}{\partial x_i} \sum_{p,q=1}^{d} G_{ip}^{k} x_{p} x_{q} G_{jq}^{k}, \qquad (D.5) $$

$$ = \sum_{i,p,q=1}^{d} G_{ip}^{k}(\delta_{ip}x_{q} + x_{p}\delta_{iq})G_{jq}^{k}, \qquad (D.6) $$

$$ = \sum_{p,q=1}^{d} G_{pp}^{k} x_{q} G_{jq}^{k} + G_{qp}^{k} x_{p} G_{jq}^{k}, \qquad (D.7) $$

$$ = [\mathrm{tr}(\mathbf{G}^k)\mathbf{G}^k x + \mathbf{G}^k(\mathbf{G}^k x)]_j. \qquad (D.8) $$

By skew-symmetry trace($\mathbf{G}^k$) = 0 and consequently

$$ \nabla \cdot (\Sigma(\mathbf{x})) = \sum_{k=1}^{d} (\nabla \cdot (\Sigma^k(\mathbf{x}))) = \sum_{k=1}^{d} \mathbf{G}^k(\mathbf{G}^k x). $$

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

**Lemma D.1.2. (Forward SDE - Itô)** Let the skew-symmetry assumption equation A1 hold. Then, the Itô form of the forward SDE equation 3.1 of $\vec{\mathbf{x}}_s$ is given by

$$ d\vec{\mathbf{x}}_s = \frac{1}{2}(\nabla \cdot \Sigma)(\vec{\mathbf{x}}_s)ds + \mathbf{G}(\vec{\mathbf{x}}_s)d\vec{\mathbf{B}}_s. \qquad (D.9) $$

*Proof.* Using the standard Stratonovich-to-Itô formula (e.g. Kunita, 1997), it holds

$$ d\vec{\mathbf{x}}_s = \frac{1}{2} d\langle \mathbf{G}(\vec{\mathbf{x}}_s), \vec{\mathbf{B}}_s \rangle_s + \mathbf{G}(\vec{\mathbf{x}}_s) d\vec{\mathbf{B}}_s, \qquad (D.10) $$

$$ = \frac{1}{2} \sum_{k=1}^{d} d\langle \mathbf{G}^k \vec{\mathbf{x}}_s, (\vec{\mathbf{B}}_s)_k \rangle_s + \mathbf{G}(\vec{\mathbf{x}}_s) d\vec{\mathbf{B}}_s, \qquad (D.11) $$

$$ = \frac{1}{2} \sum_{k=1}^{d} (\mathbf{G}^k)^2 \vec{\mathbf{x}}_s ds + \mathbf{G}(\vec{\mathbf{x}}_s) d\vec{\mathbf{B}}_s, \qquad (D.12) $$

$$ = \frac{1}{2}(\nabla \cdot \Sigma)(\vec{\mathbf{x}}_s)ds + \mathbf{G}(\vec{\mathbf{x}}_s)d\vec{\mathbf{B}}_s, \qquad (D.13) $$

where the last equality comes from Lemma D.1.1.

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

## D.2 Fokker-Planck equation and Theorem 3.1.1

Let $\vec{x}_0 \sim p_{\vec{x}_0}$ with $p_{\vec{x}_0} \in C^2(\mathbb{R}^d)$ be twice continuously differentiable. Let $p_s$ denote the density of the distribution of $\vec{x}_s$. For each $\mathbf{x} \in \mathbb{R}^d$ we define the normalized vector $\mathbf{x}^n := \frac{\mathbf{x}}{\|\mathbf{x}\|}$ for $\mathbf{x} \neq 0$ and 0 otherwise, which is orthogonal to the $d$-sphere $\mathbb{S}^{d-1}$. Furthermore, let $\nabla_{\perp}$ be the orthogonal projection of the gradient $\nabla$ in the tangent plane $\mathbf{x}^{\perp} = T_{\mathbf{x}}\mathbb{S}^{d-1}$, the tangent space on the Riemannian manifold $\mathbb{S}^{d-1}$ at the point $\mathbf{x}$, defined for $f \in C^2(\mathbb{R}^d)$ as

$$ \nabla_{\perp} f := \nabla f - (\mathbf{x}^n \cdot \nabla f) \mathbf{x}^n. \qquad (D.14) $$

**Lemma D.2.1.** It holds for any smooth vector field $f$ that

$$ \nabla \cdot f = (\mathbf{x}^n \cdot \nabla)(\mathbf{x}^n \cdot f) + \nabla_{\perp} \cdot f. \qquad (D.15) $$

*Proof.* Let $f$ be a smooth vector field, then

$$ (\mathbf{x}^n(\mathbf{x}^n)^{\top}\nabla) \cdot f = \sum_{i=1}^{d} \sum_{j=1}^{d} x_i^n x_j^n \frac{\partial}{\partial x_j} f_i = \sum_{i=1}^{d} (\mathbf{x}^n)_i \langle \mathbf{x}^n, \nabla f \rangle f_i = \langle \langle \mathbf{x}^n, \nabla f \rangle, \mathbf{x}^n \rangle. \qquad (D.16) $$
It holds for each $j = 1, \dots, d$ that

$$[(\mathbf{x}^n \cdot \nabla)\mathbf{x}^n]_j = \left[ \left( \sum_{i=1}^d (\mathbf{x}^n)_i \frac{\partial}{\partial x_i} \right) \frac{\mathbf{x}}{\|\mathbf{x}\|} \right]_j = (\mathbf{x}^n)_j \frac{1}{\|\mathbf{x}\|} - \sum_{i=1}^d (\mathbf{x}^n)_i \frac{x_i}{\|\mathbf{x}\|^3} x_j, \quad (\text{D.17})$$

$$= x_j / \|\mathbf{x}\|^2 - \left( \sum_{i=1}^{d} x_i^2 / \|\mathbf{x}\|^4 \right) x_j, \quad (\text{D.18})$$

$$= 0. \quad (\text{D.19})$$

Consequently,

$$(\mathbf{x}^n \cdot \nabla)(\mathbf{x}^n \cdot f) = (\mathbf{x}^n \cdot \nabla)\mathbf{x}^n \cdot f + (\mathbf{x}^n \cdot \nabla)f \cdot \mathbf{x}^n = (\mathbf{x}^n \cdot \nabla)f \cdot \mathbf{x}^n. \quad (\text{D.20})$$

Using the decomposition $\nabla = \mathbf{x}^n(\mathbf{x}^n)^\top\nabla + (I - \mathbf{x}^n(\mathbf{x}^n)^\top)\nabla = \mathbf{x}^n(\mathbf{x}^n)^\top\nabla + \nabla_\perp$, equation D.16 and equation D.20 we conclude

$$\nabla \cdot f = (\mathbf{x}^n(\mathbf{x}^n)^\top\nabla + \nabla_\perp) \cdot f = (\mathbf{x}^n \cdot \nabla)f \cdot \mathbf{x}^n + \nabla_\perp \cdot f = (\mathbf{x}^n \cdot \nabla)(\mathbf{x}^n \cdot f) + \nabla_\perp \cdot f.$$

□

Define the conditional noise covariance $\Sigma(\mathbf{x})$ as

$$\Sigma(\mathbf{x}) := \mathbf{G}(\mathbf{x})\mathbf{G}(\mathbf{x})^\top = \mathbb{E}[(d\vec{\mathbf{x}}_s)(d\vec{\mathbf{x}}_s)^\top|\vec{\mathbf{x}}_s = \mathbf{x}]. \quad (\text{D.21})$$

We can now state and proof Theorem 3.1.1.

**Theorem D.2.1.** Let the assumptions A1 and A2 hold. Then, the Fokker-Planck equation of equation 3.1 reads

$$\left\{ \begin{array}{l} \frac{\partial}{\partial s} p_s(\mathbf{x}) = \nabla_{\perp} \cdot \left( \frac{1}{2} \Sigma(\mathbf{x}) \nabla_{\perp} p_s(\mathbf{x}) \right), \quad \mathbf{x} \in \mathbb{R}^d, \\ p_0 = p_{\vec{\mathbf{x}}_0}. \end{array} \right. \quad (\text{D.22})$$

Moreover, any stationary density $p_\infty$ of equation D.22 is rotation-invariant on $\mathbb{R}^d$.

*Proof.* From the Itô SDE equation D.9, the Fokker-Planck equation describing the evolution of $(p_s)_{s \ge 0}$ is given as

$$\frac{\partial}{\partial s} p_s = \nabla \cdot \left( -\frac{1}{2} \nabla \cdot (\Sigma(\mathbf{x})) p_s(\mathbf{x}) + \frac{1}{2} \nabla \cdot (\Sigma(\mathbf{x}) p_s(\mathbf{x})) \right), \quad (\text{D.23})$$

$$= \nabla \cdot \left( -\frac{1}{2} \nabla \cdot (\Sigma(\mathbf{x})) p_s(\mathbf{x}) + \frac{1}{2} \nabla \cdot (\Sigma(\mathbf{x})) p_s(\mathbf{x}) + \frac{1}{2} \Sigma(\mathbf{x}) \nabla p_s(\mathbf{x}) \right), \quad (\text{D.24})$$

$$= \nabla \cdot \left( \frac{1}{2} \Sigma(\mathbf{x}) \nabla p_s(\mathbf{x}) \right). \quad (\text{D.25})$$

The skew-symmetry condition in assumption A1 implies for any $\mathbf{x} \in \mathbb{R}^d$ that

$$((\mathbf{x}^n)^\top \Sigma(\mathbf{x}))^\top = \Sigma(\mathbf{x}) \mathbf{x}^n = \sum_{k=1}^d G_k \mathbf{x} \|\mathbf{x}\| \underbrace{(\mathbf{x}^n)^\top G_k^\top \mathbf{x}^n}_{=0} = 0. \quad (\text{D.26})$$

Combining this with the result of equation D.16, it holds that

$$\nabla \cdot (\Sigma(\mathbf{x}) \nabla p_s(\mathbf{x})) = (\mathbf{x}^n \cdot \nabla)(\mathbf{x}^n \cdot \Sigma(\mathbf{x}) \nabla p_s(\mathbf{x})) + \nabla_\perp \cdot (\Sigma(\mathbf{x}) \nabla p_s(\mathbf{x})). \quad (\text{D.27})$$

The decomposition $\nabla = \mathbf{x}^n(\mathbf{x}^n)^\top\nabla + (I - \mathbf{x}^n(\mathbf{x}^n)^\top)\nabla = \mathbf{x}^n(\mathbf{x}^n)^\top\nabla + \nabla_\perp$ and equation D.26 yields

$$\nabla_\perp \cdot (\Sigma(\mathbf{x}) \nabla p_s(\mathbf{x})) = \nabla_\perp \cdot (\Sigma(\mathbf{x}) \nabla_\perp p_s(\mathbf{x})) + \nabla_\perp \cdot (\Sigma(\mathbf{x}) \mathbf{x}^n (\mathbf{x}^n)^\top \nabla p_s(\mathbf{x})), \quad (\text{D.28})$$

$$= \nabla_\perp \cdot (\Sigma(\mathbf{x}) \nabla_\perp p_s(\mathbf{x})). \quad (\text{D.29})$$

Hence, by linearity

$$\frac{\partial}{\partial s} p_s = \nabla_\perp \cdot \left( \frac{1}{2} \Sigma(\mathbf{x}) \nabla_\perp p_s(\mathbf{x}) \right). \quad (\text{D.30})$$
We shall now explore the set of possible invariant densities $\rho_{\infty}$ of the Fokker-Planck equation D.22. We will show that $\rho_{\infty}$ is stationary if and only if it is rotation-invariant.

Let $p_{\infty}$ be rotation-invariant, i.e. $\nabla_{\perp} p_{\infty} = 0$ then it is a stationary solution of the Fokker-Planck equation. The set of rotation-invariant measures is not empty, e.g. containing the isotropic normal distributions $N(0, I_d)$.

Conversely, let $p_{\infty}$ be a stationary solution of the Fokker-Planck equation, in particular

$$ \nabla_{\perp} \cdot \left( \frac{1}{2} \Sigma(\mathbf{x}) \nabla_{\perp} p_{\infty}(\mathbf{x}) \right) = 0. \qquad (D.31) $$

Integrating over the test function $\phi = p_{\infty}$ gives a necessary condition for $p_{\infty}$ to be an invariant measure:

$$ 0 = - \int_{\mathbb{R}^d} p_{\infty}(\mathbf{x}) \nabla_{\perp} \cdot (\Sigma(\mathbf{x}) \nabla_{\perp} p_{\infty}(\mathbf{x})) d\mathbf{x}, \qquad (D.32) $$

$$ = \int_{\mathbb{R}^d} \nabla_{\perp} p_{\infty}(\mathbf{x})^{\top} \Sigma(\mathbf{x}) \nabla_{\perp} p_{\infty}(\mathbf{x}) d\mathbf{x}, \qquad (D.33) $$

$$ = \int_{\mathbb{R}^d} \underbrace{\| \mathbf{G}^{\top}(\mathbf{x}) \nabla_{\perp} p_{\infty}(\mathbf{x}) \|^2}_{\ge 0} d\mathbf{x}. \qquad (D.34) $$

Hence, for a.e. $\mathbf{x} \in \mathbb{R}^d$ it holds that $\nabla_{\perp} p_{\infty}(\mathbf{x}) \in \ker(\mathbf{G}^{\top}(\mathbf{x}))$. By assumption A2, this kernel has dimension 1. Moreover, $\mathbf{G}^{\top}(\mathbf{x})\mathbf{x} = 0$ and by definition of $\nabla_{\perp}$ we have that $\nabla_{\perp} p_{\infty}(\mathbf{x}) \perp \mathbf{x}$. That means

$$ \nabla_{\perp} p_{\infty}(\mathbf{x}) \in \ker(\mathbf{G}^{\top}(\mathbf{x})) \cap \mathbf{x}^{\perp} = \operatorname{span}\{\mathbf{x}\} \cap \mathbf{x}^{\perp} = \{0\}. \qquad (D.35) $$

We conclude that $\nabla_{\perp} p_{\infty}(\mathbf{x}) = 0$ almost everywhere on $\mathbb{R}^d$, i.e. the measure is rotation-invariant. $\square$

## D.3 Distribution of the norms

In this section, we see more precisely that the norm of the MSGM SDE solution remains constant along the noising process whereas in SGM the norm dynamics is random with a mean going to $\sqrt{d}$ asymptotically.

### D.3.1 Norm dynamics in SGM

The dynamics of the SGM diffusion norm is stochastic. The following proposition states that the norm of the SGM latent concentrates around its mean, $\sqrt{d}$, for large dimension $d$.

**Proposition D.3.1.** If $\vec{\mathbf{x}}_s$ is an Ornstein Uhlenbeck process then

$$ \mathbb{E} \left[ \| \vec{\mathbf{x}}_s \|^2 | \vec{\mathbf{x}}_0 \right] = e^{-2s} \| \vec{\mathbf{x}}_0 \|^2 + (1 - e^{-2s}) d \underset{s \to \infty}{\longrightarrow} d. \qquad (D.36) $$

and

$$ \| \vec{\mathbf{x}}_s \|^2 = \mathbb{E} \left[ \| \vec{\mathbf{x}}_s \|^2 | \vec{\mathbf{x}}_0 \right] + \sqrt{d} I_s + e^{-s} K_s = d \left( 1 + \frac{1}{\sqrt{d}} \mathcal{O} \left( \frac{1}{s} \right) \right), \qquad (D.37) $$

with both $\mathbb{E}K_s^2$ and $\mathbb{E}I_s^2$ bounded for large time $s$ and $\mathbb{E}I_s^2$ independent of the dimension $d$.

*Proof.* To get the dynamics of the squared norm mean in SGM, we can take the expectation of the following Itô equation

$$ d \| \vec{\mathbf{x}}_s \|^2 = 2 \vec{\mathbf{x}}_s \cdot d \vec{\mathbf{x}}_s + d \langle \vec{\mathbf{x}}^{\top}, \vec{\mathbf{x}} \rangle_s = 2 (\| \vec{\mathbf{x}}_s \|^2 - d) ds + 2 \sqrt{2} \vec{\mathbf{x}}_s \cdot d \vec{\mathbf{B}}_s, \quad \forall s \ge 0. \qquad (D.38) $$

Thus,

$$ \mathbb{E} \left[ \| \vec{\mathbf{x}}_s \|^2 | \vec{\mathbf{x}}_0 \right] = e^{-2s} \| \vec{\mathbf{x}}_0 \|^2 + (1 - e^{-2s}) d \underset{s \to \infty}{\longrightarrow} d. \qquad (D.39) $$

To obtain the full nom dynamics from equation D.38, we note that $t \to e^{-2s}$ has finite variations. Accordingly,

$$ d(e^{-2s} (\| \vec{\mathbf{x}}_s \|^2 - d)) = 2\sqrt{2} e^{-2s} \vec{\mathbf{x}}_s \cdot d \vec{\mathbf{B}}_s, \qquad (D.40) $$
and a temporal integration and the analytic expression of the Ornstein Uhlenbeck process yields:

$$ \|\vec{\mathbf{x}}_s\|^2 = \mathbb{E} \left[ \|\vec{\mathbf{x}}_s\|^2 | \vec{\mathbf{x}}_0 \right] + 2\sqrt{2} \int_0^s e^{-2(s-s')} \vec{\mathbf{x}}_{s'} \cdot d\vec{\mathbf{B}}_{s'}, \quad (D.41) $$

$$ = \mathbb{E} \left[ \|\vec{\mathbf{x}}_s\|^2 |\vec{\mathbf{x}}_0\right] + 2\sqrt{2} e^{-s} \vec{\mathbf{x}}_0 \cdot \int_0^s e^{-(s-s')} d\vec{\mathbf{B}}_{s'} \\ + 2\sqrt{2} \int_0^s \int_0^{s'} e^{-(2s-s'-s'')} d\vec{\mathbf{B}}_{s'} \cdot d\vec{\mathbf{B}}_{s''}, \quad (D.42) $$

$$ = \mathbb{E} \left[ \|\vec{\mathbf{x}}_s\|^2 |\vec{\mathbf{x}}_0\right] + \sqrt{d} I_s + e^{-s} K_s, \quad (D.43) $$

with the martingales

$$ I_s = \sqrt{\frac{8}{d}} \int_0^s \int_0^{s'} e^{-(2s-s'-s'')} d\vec{\mathbf{B}}_{s'} \cdot d\vec{\mathbf{B}}_{s''}, \quad (D.44) $$

$$ K_s = \sqrt{8} \vec{\mathbf{x}}_0 \cdot \int_0^s e^{-(s-s')} d\vec{\mathbf{B}}_{s'}. \quad (D.45) $$

$K_s$ corresponds to the martingale part of the Ornstein Uhlenbeck solution projected on $\vec{\mathbf{x}}_0$. It is well known that $\mathbb{E}K_s^2$ is bounded for large time $s$. $\mathbb{E}I_s^2$ may be less known and we shall evaluate it below:

$$ \mathbb{E}I_s^2 = \frac{8}{d} e^{-4s} \mathbb{E} \left( \sum_{p=1}^{d} \int_0^s \int_0^{s'} e^{s'+s''} d(\vec{B}_{s'})_p \cdot d(\vec{B}_{s''})_p \right)^2, \quad (D.46) $$

$$ = \frac{8}{d} e^{-4s} \mathbb{E} \sum_{p_1, p_2=1}^{d} \int_0^s \int_0^{s_1'} \int_0^s \int_0^{s_2'} e^{s_1' + s_1'' s_2' + s_2''} d(\vec{B}_{s_1'})_{p_1} \cdot d(\vec{B}_{s_1''})_{p_1} d(\vec{B}_{s_2''})_{p_2} \cdot d(\vec{B}_{s_2'})_{p_2}, \quad (D.47) $$

$$ = \frac{8}{d} e^{-4s} \mathbb{E} \sum_{p_1, p_2=1}^{d} \int_0^s \int_0^{s_1'} \int_0^s \int_0^{s_2'} e^{s_1' + s_1'' + s_2' + s_2''} \delta_{p_1, p_2} \delta(s_1' - s_2') \delta_{p_1, p_2} \delta(s_1'' - s_2'') ds_1' ds_1'' ds_2' ds_2'' \quad (D.48) $$

$$ = 8e^{-4s} \int_0^s \int_0^{s'} e^{2(s'+s'')} ds' ds'', \quad (D.49) $$

$$ = 4e^{-4s} \int_0^s e^{2s'} (e^{2s'} - 1) ds', \quad (D.50) $$

$$ = e^{-4s} ((e^{4s} - 1) + 2(e^{2s} - 1)), \quad (D.51) $$

$$ \underset{s \to \infty}{\longrightarrow} 1. \quad (D.52) $$

□

### D.3.2 Norm dynamics in MSGM

For MSGM, the norm follows totally different dynamics. We recall that the skew-symmetry of $\circ d\mathbb{Z}_s$ implies that $d\vec{\mathbf{x}}_s = \circ d\mathbb{Z}_s \vec{\mathbf{x}}_s$ is orthogonal to $\vec{\mathbf{x}}_s$ and hence:

$$ d\|\vec{\mathbf{x}}_s\|^2 = 2\vec{\mathbf{x}}_s \cdot \circ d\vec{\mathbf{x}}_s = 0, \quad \forall s \ge 0. \quad (D.53) $$

Consequently, $\vec{\mathbf{x}}_s$ moves randomly on $\|\vec{\mathbf{x}}_0\|\mathbb{S}^{d-1}$, the $d$-sphere of radius $\|\vec{\mathbf{x}}_0\|$, and the increments $d\vec{\mathbf{x}}_s$ are tangent to the $d$-sphere. In particular, we obtain the following result.

**Proposition D.3.2.** Let the skew-symmetry assumption A1 hold. Let $\vec{\mathbf{x}}_0$ be a random variable. Then, for all $s \ge 0$ the distribution of $\|\vec{\mathbf{x}}_s\|$ equals the distribution of $\|\vec{\mathbf{x}}_0\|$.

Therefore, the distribution of the norms of the latent variable is exactly the distribution of the norms of the points of the dataset. Moreover, $\vec{\mathbf{x}}_s \equiv 0$ if and only if $\vec{\mathbf{x}}_0 = 0$. As a consequence, we can exclude all points exactly equal to zero from a dataset, treat them aside, and hence consider, without loss of generality, that $\vec{\mathbf{x}}_T \neq 0$ almost surely.
## D.4 Fokker-Planck equation of the direction

This subsection is devoted to the analysis of the Fokker-Planck equation on the unit sphere $S^{d-1}$, i.e. the distribution of $\vec{x}_s^n$, in particular as $s \to \infty$.

### D.4.1 Main results on the distribution of directions

We saw in Appendix D.3 that $\vec{x}_s$ moves randomly on the $d$-sphere of radius $\|\vec{x}_0\|$ and that the increments, $d\vec{x}_s = \mathbf{G}(\vec{x}_s) \circ d\vec{B}_s$, are tangent to the $d$-sphere. If the rank condition, assumption A2 is verified, then the support of the noise distribution $\mathbf{G}(\vec{x}_s) \circ d\vec{B}_s$ coincides with the $d-1$-dimensional tangent space, i.e. it will likely explore all local directions around $\vec{x}_s$. With time, the support of the solution distribution will gradually cover the whole $d$-sphere, i.e. every direction $\vec{x}_s^n$ will become equiprobable. Lemma D.4.1 illustrates and precises this claim.

**Lemma D.4.1.** Let assumptions A1 and A2 hold. Let an initial density $p_0^n \in C^2(S^{d-1})$ and $\Sigma(\mathbf{x}^n) := \mathbf{G}(\mathbf{x}^n)\mathbf{G}(\mathbf{x}^n)^T$. Then, the Fokker-Planck equation

$$ \frac{\partial}{\partial s} p_s^n(\mathbf{x}^n) = \nabla_{\perp} \cdot (\frac{1}{2} \Sigma(\mathbf{x}^n) \nabla_{\perp} p_s^n(\mathbf{x}^n)), \quad \mathbf{x}^n \in S^{d-1}, \qquad (D.54) $$

has a unique density solution $p_s^n \in C^2(S^{d-1})$ for all $s > 0$. Moreover, there is a unique invariant measure $p_\infty^n$ of that Fokker-Planck equation, i.e. the uniform distribution on the $d$-sphere $U(S^{d-1})$, with density

$$ p_{\infty}^{n}(\mathbf{x}^{n}) := \frac{1}{|S^{d-1}|}, \quad \forall \mathbf{x}^{n} \in S^{d-1}, \qquad (D.55) $$

with $|S^{d-1}| = 2\pi^{d/2}/\Gamma(d/2)$ the volume of the $d$-sphere $S^{d-1}$ and $\Gamma$ the gamma function.

Lemma D.4.1 is a consequence of Theorem 3.1.1 as shown in Appendix D.4.2. Note that in this case $\nabla_{\perp} = \nabla_{S^{d-1}}$ is the Riemannian gradient on $S^{d-1}$, see Appendix H.1.

Given the unique invariant measure of Fokker-Planck equation formulated on $S^{d-1}$, we can also show that we have exponential convergence of the initial distribution $p_0^n$ to $p_\infty^n$, the uniform distribution on the unit sphere $S^{d-1}$.

**Theorem D.4.1.** Let assumptions A1 and A2 hold. Then, there exists $\alpha = \alpha(\mathbf{G}, d) > 0$ with

$$ \|p_s^n - p_\infty^n\|_{L^2(S^{d-1})}^2 \le \exp(-\alpha s) \|p_0^n - p_\infty^n\|_{L^2(S^{d-1})}^2. \qquad (D.56) $$

The convergence rate $\alpha$ is given as

$$ \alpha(\mathbf{G}, d) = (d-1) \min_{(\mathbf{x}, \mathbf{y}) \in S} \|\mathbf{G}^T(\mathbf{x})\mathbf{y}\|^2, \quad S = \{(\mathbf{x}, \mathbf{y}) \in S^{d-1} \times S^{d-1} | \mathbf{x} \perp \mathbf{y}\}. \qquad (D.57) $$

Consequently, since $S^{d-1}$ is compact this implies convergence in total variation of $p_s^n$ to $p_\infty^n$ and convergence in distribution of $\vec{x}_s^n$ with $\vec{x}_\infty^n \sim U(S^{d-1})$. The full proof is detailed in Appendix D.4.3.

### D.4.2 Proof of Lemma D.4.1

*Proof. Existence and Uniqueness:*

Consider $L(p_s^n) = \nabla_{\perp} \cdot (\frac{1}{2} \Sigma(\mathbf{x})\nabla_{\perp} p_s^n(\mathbf{x})) - \frac{\partial}{\partial s} p_s^n(\mathbf{x})$. $L$ is a parabolic type operator according to Friedman (1964) since $\mathbf{x} \mapsto \Sigma(\mathbf{x})$ is positive definite by assumption A2 on $S^{d-1}$. Indeed, for any $\mathbf{y} \in T_x S^{d-1}$ the tangential (linear) space of $S^{d-1}$ at $\mathbf{x}$,

$$ \mathbf{y}^T \Sigma(\mathbf{x}) \mathbf{y} = \| \mathbf{G}^T(\mathbf{x}) \mathbf{y} \|^2 \ge 0. \qquad (D.58) $$

with equality if and only if $\mathbf{G}^T(\mathbf{x})\mathbf{y} = 0$. Then, the rank condition A2 implies $\mathbf{y} = 0$ as previously in equation D.35. Consequently, the associated spatial operator $L_0$ defined by

$$ L_0 p_s^n = \nabla_{\perp} \cdot \left( \frac{1}{2} \Sigma(\mathbf{x}) \nabla_{\perp} p_s^n(\mathbf{x}) \right) \qquad (D.59) $$

is an elliptic operator on $S^{d-1}$, a compact manifold without boundary such that the semi-group $e^{sL_0}$ is strongly continuous on $C^2(S^{d-1})$, $s \ge 0$. As $p_0^n \in C^2(S^{d-1})$, according to chapter 1, proposition
1.1 in Taylor (2011), there exists a unique solution $p_s^n \in C^2(\mathbb{S}^{d-1})$, for $s \in [0, T]$ of equation D.54. As the semigroup is well-defined for all $s > 0$, this extends the uniqueness of the solution to all $s > 0$.

**Invariant measure:** Repeating the lines in the proof of Theorem 3.1.1 given in Appendix D.2 it follows that $p_\infty^n$ is rotation-invariant. The only rotation-invariant distribution on the $d$-sphere is the uniform distribution. $\square$

### D.4.3 Proof of Theorem D.4.1: Limit behavior of Fokker-Planck equation of the direction

**Theorem D.4.2.** Let assumptions A1 and A2 hold. Then, there exists $\alpha = \alpha(\mathbf{G}, d) > 0$ with

$$ \|p_s^n - p_\infty^n\|_{L^2(\mathbb{S}^{d-1})}^2 \le \exp(-\alpha s) \|p_0^n - p_\infty^n\|_{L^2(\mathbb{S}^{d-1})}^2. \quad (D.60) $$

The convergence rate $\alpha$ is given as

$$ \alpha(\mathbf{G}, d) = (d-1) \min_{(x,y) \in S} \|\mathbf{G}^\top(x)y\|^2, \quad S = \{(x, y) \in \mathbb{S}^{d-1} \times \mathbb{S}^{d-1} | x \perp y\}. \quad (D.61) $$

**Proof.** Let $p_s^n$ denoting the density of $x_s^n$. Define $e_s^n = p_s^n - p_\infty^n$ with $p_\infty^n = |\mathbb{S}^{d-1}|^{-1}$ being the uniform distribution on $\mathbb{S}^{d-1}$. Then, by linearity of the Fokker-Planck equation, $e_t^n$ satisfies

$$ \partial_t e_t^n = \nabla_\perp \cdot \left( \frac{1}{2} \mathbf{\Sigma}(x) \nabla_\perp e_t^n(x) \right). \quad (D.62) $$

Since $p_s^n$ and $p_\infty^n$ are densities on $\mathbb{S}^{d-1}$, we have $\int_{\mathbb{S}^{d-1}} e_s^n d\mathbf{x} = 0$ for all $s \ge 0$. Consequently, since $\mathbb{S}^{d-1}$ is a compact manifold without boundary, Poincaré inequality holds, i.e.

$$ \|e_t^n\|_{L^2(\mathbb{S}^{d-1})}^2 \le \frac{1}{d-1} \|\nabla_{\mathbb{S}^{d-1}} e_t^n\|_{L^2(\mathbb{S}^{d-1})}^2, \quad (D.63) $$

with

$$ \nabla_{\mathbb{S}^{d-1}} e_t^n(y)|_{y=x} = \operatorname{Proj}_{\mathcal{T}_{x,\mathbb{S}^{d-1}}} \nabla e_t^n(y)|_{y=x} = \nabla_{\perp} e_t^n(x). \quad (D.64) $$

Consequently, integration by part on $\mathbb{S}^{d-1}$ leads to

$$ \frac{1}{2} \frac{d}{dt} \|e_t^n\|_{L^2(\mathbb{S}^{d-1})}^2 = \int_{\mathbb{S}^{d-1}} e_t^n(x) \nabla_\perp \cdot \left( \frac{1}{2} \mathbf{\Sigma}(x) \nabla_\perp e_t^n(x) \right) dx, \quad (D.65) $$

$$ = - \int_{\mathbb{S}^{d-1}} \nabla_\perp e_t^n(x)^\top \mathbf{\Sigma}(x) \nabla_\perp e_t^n(x) dx. \quad (D.66) $$

We will now bound $y^\top \mathbf{\Sigma}(x)y^\top$ from below for any $y \in x^\perp$ and $x \in \mathbb{S}^{d-1}$, in particular with a bound independent of $x$. Since $\mathbf{\Sigma}(x)$ is symmetric, it is real diagonalizable with eigen-basis denoted as $v_1(x), \dots, v_d(x) \in \mathbb{R}^d$ and eigenvalues $\lambda_1(x), \dots, \lambda_d(x)$. By construction $\mathbf{\Sigma}(x)x = 0$, hence we can set $v_d(x) := x/\|x\|$ and $\lambda_d(x) \equiv 0$. Moreover, by the rank condition A2, $\lambda_i \ne 0$ for $i \ne d$. By orthonormality of the eigenvectors $v_1(x), \dots, v_{d-1}(x)$ then span the tangent plane $x^\perp$ at $x$ on $\mathbb{S}^{d-1}$. For any $i = 1, \dots, d-1$, we have that

$$ \lambda_i(x) = v_i(x)^\top \mathbf{\Sigma}(x) v_i(x) = \|\mathbf{G}^\top(x) v_j(x)\|^2 \ge \min_{y \in x^\perp} \frac{\|\mathbf{G}^\top(x)y\|^2}{\|y\|^2}. \quad (D.67) $$

The polynomial $(x, y) \mapsto P(x, y) = \|\mathbf{G}^\top(x)y\|^2$ on the compact $S = \{(x, y) \in \mathbb{S}^{d-1} \times \mathbb{S}^{d-1} | x \perp y\}$ attains its minimum $P^*$, which from the rank condition satisfies $P(x, y) \ge P^* > 0$ for all $(x, y) \in S$. As a consequence $\lambda_i(x) \ge P^*$ for $i = 1, \dots, d-1$ and $x \in \mathbb{S}^{d-1}$, which implies for all $y \in x^\perp$ that

$$ y^\top \mathbf{\Sigma}(x)y = \|\mathbf{G}^\top(x)y\|^2 \ge P^* \|y\|^2. \quad (D.68) $$

Therefore, combining equation D.63, equation D.64 and equation D.66, we obtain

$$ \frac{1}{2} \frac{d}{dt} \|e_t^n\|_{L^2(\mathbb{S}^{d-1})}^2 \le -P^* \|\nabla_\perp e_s^n\|_{L^2(\mathbb{S}^{d-1})}^2 \le -P^*(d-1) \|e_s^n\|_{L^2(\mathbb{S}^{d-1})}^2. \quad (D.69) $$

Then, by Gronwall for $\alpha = P^*(d-1) > 0$, we conclude that

$$ \|p_s^n - p_\infty^n\|_{L^2(\mathbb{S}^{d-1})}^2 = \|e_s^n\|_{L^2(\mathbb{S}^{d-1})}^2 \le \|e_0^n\|_{L^2(\mathbb{S}^{d-1})}^2 \exp(-\alpha s). \quad (D.70) $$

$\square$
### D.5 Proof of Theorem 3.3.1: Convergence of Fokker-Planck equation

This section is devoted to the analysis of the Fokker-Planck equation in the whole domain $\mathbb{R}^d$. Due to the fact that the norm of a point does not change in the SDE process as shown in equation 3.7 and the fact that $\Sigma(\mathbf{x}) = 0$ for $\mathbf{x} = 0$, we exclude the origin in the analysis.

**Theorem D.5.1.** Let $D = \mathbb{R}^d \setminus \{0\}$ for $d > 1$. Let assumptions A1 and A2 hold. Let $\vec{\mathbf{x}}_0 \sim p_0 \in C^2(D)$ and let $p_{|\cdot|}$ be the (radial) density of $\|\vec{\mathbf{x}}_0\|$. Then, the Fokker-Planck equation

$$ \frac{\partial}{\partial s} p_s(\mathbf{x}) = \nabla_{\perp} \cdot \left( \frac{1}{2} \Sigma(\mathbf{x}) \nabla_{\perp} p_s(\mathbf{x}) \right), \quad \mathbf{x} \in D, \qquad (D.71) $$

has a unique solution $p_s \in C^2(D) \cap L^2(D)$ for all $s > 0$. Moreover, the Fokker-Planck equation has the stationary distribution

$$ p_{\infty}(\mathbf{x}) = \frac{p_{|\cdot|}(\|\mathbf{x}\|)}{\|\mathbf{x}\|^{d-1} |\mathbb{S}^{d-1}|} \qquad (D.72) $$

In particular, $\|\vec{\mathbf{x}}_s\|$ and $\vec{\mathbf{x}}_s^n$ are asymptotically independent for $s \to +\infty$. Moreover, there exists $\alpha = \alpha(\mathbf{G}, d) > 0$ such that

$$ \|p_s - p_{\infty}\|_{L^2(\mathbb{R}^d)}^2 \leq \exp(-\alpha s) \|p_0 - p_{\infty}\|_{L^2(\mathbb{R}^d)}^2. $$

The convergence rate $\alpha$ is given as

$$ \alpha(\mathbf{G}, d) = (d-1) \min_{(\mathbf{x}, \mathbf{y}) \in S} \|\mathbf{G}^{\top}(\mathbf{x})\mathbf{y}\|^2, \quad S = \{(\mathbf{x}, \mathbf{y}) \in \mathbb{S}^{d-1} \times \mathbb{S}^{d-1} | \mathbf{x} \perp \mathbf{y}\}. \qquad (D.73) $$

*Proof.* We will proof existence, uniqueness, regularity, invariant property and convergence separately.

*Existence:* Let $p_0(\mathbf{x}^n \| \vec{\mathbf{x}}_0 \| = r)$ be the start value of the FP equation D.54 on $\mathbb{S}^{d-1}$ of Lemma D.4.1. This gives rise to a smooth unique density solution $p_s^n(\mathbf{x}^n \| \vec{\mathbf{x}}_0 \| = r)$ for $s > 0$ and any $r > 0$. Moreover, $p_s^n(\mathbf{x}^n \| \vec{\mathbf{x}}_0 \| = r) = p_s^n(\mathbf{x}^n \| \vec{\mathbf{x}}_s \| = r)$ since $d \| \vec{\mathbf{x}}_s \| = 0$. Now define

$$ \rho_s(\mathbf{x}) = p_s^n(\mathbf{x}^n \| \vec{\mathbf{x}}_s \| = \|\mathbf{x}\|) p_{|\cdot|}(\|\mathbf{x}\|) \|\mathbf{x}\|^{1-d}, $$

where we denote $\mathbf{x}^n = \frac{\mathbf{x}}{\|\mathbf{x}\|}$. $\rho_s$ is a density since

$$ \int_{\mathbb{R}^d} \rho_s(\mathbf{x}) \, d\mathbf{x} = \int_{\mathbb{R}^d} p_s^n(\mathbf{x}^n \| \vec{\mathbf{x}}_s \| = \|\mathbf{x}\|) p_{|\cdot|}(\|\mathbf{x}\|) \|\mathbf{x}\|^{1-d} \, d\mathbf{x}, \qquad (D.74) $$

$$ = \int_{\mathbb{R}_+} \int_{\mathbb{S}^{d-1}} p_s^n(\mathbf{x}^n \| \vec{\mathbf{x}}_s \| = r) p_{|\cdot|}(r) r^{1-d} r^{d-1} \, dr \, d\mathbf{x}^n, \qquad (D.75) $$

$$ = \int_{\mathbb{R}_+} \int_{\mathbb{S}^{d-1}} p_s^n(\mathbf{x}^n \| \vec{\mathbf{x}}_s \| = r) p_{|\cdot|}(r) \, dr \, d\mathbf{x}^n, \qquad (D.76) $$

$$ = \int_{\mathbb{R}_+} p_{|\cdot|}(r) \left( \int_{\mathbb{S}^{d-1}} p_s^n(\mathbf{x}^n \| \vec{\mathbf{x}}_s \| = r) \, d\mathbf{x}^n \right) \, dr, \qquad (D.77) $$

$$ = 1. \qquad (D.78) $$

We have $\nabla_{\perp} = \frac{1}{\|\mathbf{x}\|} \nabla_{\mathbb{S}^{d-1}}$ and $\nabla_{\perp}$ does not act on radial functions. Besides, $\Sigma(\mathbf{x}) := \mathbf{G}(\mathbf{x})\mathbf{G}(\mathbf{x})^{\top}$ with $\mathbf{G}$ linear so $\Sigma(\mathbf{x}) = \Sigma(\|\mathbf{x}\| \frac{\mathbf{x}}{\|\mathbf{x}\|}) = \|\mathbf{x}\|^2 \Sigma(\mathbf{x}^n)$. Hence

$$ \nabla_{\perp} \rho_s(\mathbf{x}) = \nabla_{\perp} p_s^n(\mathbf{x}^n \| \vec{\mathbf{x}}_s \| = \|\mathbf{x}\|) p_{|\cdot|}(\|\mathbf{x}\|) \|\mathbf{x}\|^{1-d}, \qquad (D.79) $$

$$ = p_{|\cdot|}(\|\mathbf{x}\|) \|\mathbf{x}\|^{1-d} \nabla_{\perp} p_s^n(\mathbf{x}^n \| \vec{\mathbf{x}}_s \| = \|\mathbf{x}\|), \qquad (D.80) $$

$$ = p_{|\cdot|}(\|\mathbf{x}\|) \|\mathbf{x}\|^{1-d} \left( \frac{1}{\|\mathbf{x}\|} \nabla_{\mathbb{S}^{d-1}} \right) p_s^n(\mathbf{x}^n \| \vec{\mathbf{x}}_s \| = \|\mathbf{x}\|), \qquad (D.81) $$
and

$$ \Sigma(\mathbf{x})\nabla_{\perp}\rho_s(\mathbf{x}) = \Sigma(\mathbf{x})p_{\cdot|\cdot}(\|\mathbf{x}\|)\|\mathbf{x}\|^{1-d} \frac{1}{\|\mathbf{x}\|} \nabla_{\mathbb{S}^{d-1}} p_s^n(\mathbf{x}^n)\|\vec{\mathbf{x}}_s\| = \|\mathbf{x}\|), \quad (D.82) $$

$$ = \|\mathbf{x}\|^2 \Sigma(\mathbf{x}^n) p_{\cdot|\cdot}(\|\mathbf{x}\|) \|\mathbf{x}\|^{1-d} \frac{1}{\|\mathbf{x}\|} \nabla_{\mathbb{S}^{d-1}} p_s^n(\mathbf{x}^n) \|\vec{\mathbf{x}}_s\| = \|\mathbf{x}\|), \quad (D.83) $$

$$ = \|\mathbf{x}\|^2 p_{\cdot|\cdot}(\|\mathbf{x}\|) \|\mathbf{x}\|^{1-d} \frac{1}{\|\mathbf{x}\|} \Sigma(\mathbf{x}^n) \nabla_{\mathbb{S}^{d-1}} p_s^n(\mathbf{x}^n) \|\vec{\mathbf{x}}_s\| = \|\mathbf{x}\|) \quad (D.84) $$

$$ \nabla_{\perp} \cdot (\Sigma(\mathbf{x}) \nabla_{\perp} \rho_s(\mathbf{x})) = \frac{1}{\|\mathbf{x}\|} \nabla_{\mathbb{S}^{d-1}} \cdot \left( \|\mathbf{x}\| p_{\cdot|\cdot}(\|\mathbf{x}\|) \|\mathbf{x}\|^{1-d} \Sigma(\mathbf{x}^n) \nabla_{\mathbb{S}^{d-1}} p_s^n(\mathbf{x}^n) \|\vec{\mathbf{x}}_s\| = \|\mathbf{x}\| \right), \quad (D.85) $$

$$ = p_{\cdot|\cdot}(\|\mathbf{x}\|) \|\mathbf{x}\|^{1-d} \nabla_{\mathbb{S}^{d-1}} \cdot \left( \Sigma(\mathbf{x}^n) \nabla_{\mathbb{S}^{d-1}} p_s^n(\mathbf{x}^n) \|\vec{\mathbf{x}}_s\| = \|\mathbf{x}\| \right), \quad (D.86) $$

$$ = p_{\cdot|\cdot}(\|\mathbf{x}\|) \|\mathbf{x}\|^{1-d} 2 \frac{\partial}{\partial s} p_s^n(\mathbf{x}^n) \|\vec{\mathbf{x}}_s\| = \|\mathbf{x}\|), \quad (D.87) $$

$$ = 2 \frac{\partial}{\partial s} \left( p_{\cdot|\cdot}(\|\mathbf{x}\|) \|\mathbf{x}\|^{1-d} p_s^n(\mathbf{x}^n) \|\vec{\mathbf{x}}_s\| = \|\mathbf{x}\| \right), \quad (D.88) $$

$$ = 2 \frac{\partial}{\partial s} \rho_s(\mathbf{x}), \quad (D.89) $$

i.e. $\frac{\partial}{\partial s} \rho_s(\mathbf{x}) = \frac{1}{2} \nabla_{\perp} \cdot (\Sigma(\mathbf{x}) \nabla_{\perp} \rho_s(\mathbf{x}))$. Then, $\rho_s$ solves the Fokker-Planck equation on $\mathbb{R}^d$.

*Uniqueness:* Assume there exists another $\tilde{\rho}$ solving the FP on $\mathbb{R}^d$ verify

$$ \tilde{\rho}_s(\mathbf{x}) = \rho_{1,s}(\mathbf{x}^n \|\mathbf{x}\|) \rho_{2,s}(\|\mathbf{x}\|). $$

Since $d\|\vec{\mathbf{x}}_s\| = 0$, by marginalizing $\tilde{\rho}_s$ (integrating on $\mathbb{S}^{d-1}$), we have the uniqueness of the radial density $\rho_{2,s}(r) = \int_{\mathbb{S}^{d-1}} \tilde{\rho}_s(r\mathbf{x}^n)d\mathbf{x}^n = p_{\|\vec{\mathbf{x}}_s\|}(r)r^{1-d} = p_{\cdot|\cdot}(r)r^{1-d}$.

Since $d\|\vec{\mathbf{x}}_s\| = 0$, we have $\frac{\partial}{\partial s} \rho_{2,s}(\|\mathbf{x}\|) = 0$. Therefore,

$$ \frac{\partial}{\partial s} \tilde{\rho}_s(\mathbf{x}) = \frac{\partial}{\partial s} (\rho_{1,s}(\mathbf{x}^n \|\mathbf{x}\|) \rho_{2,s}(\|\mathbf{x}\|)), \quad (D.90) $$

$$ = \rho_{2,s}(\|\mathbf{x}\|) \frac{\partial}{\partial s} \rho_{1,s}(\mathbf{x}^n \|\mathbf{x}\|) + \rho_{1,s}(\mathbf{x}^n \|\mathbf{x}\|) \frac{\partial}{\partial s} \rho_{2,s}(\|\mathbf{x}\|), \quad (D.91) $$

$$ = \rho_{2,s}(\|\mathbf{x}\|) \frac{\partial}{\partial s} \rho_{1,s}(\mathbf{x}^n \|\mathbf{x}\|). \quad (D.92) $$

In addition,

$$ 2 \frac{\partial}{\partial s} \tilde{\rho}_s(\mathbf{x}) = \nabla_{\perp} \cdot (\Sigma(\mathbf{x}) \nabla_{\perp} \tilde{\rho}_s(\mathbf{x})), \quad (D.93) $$

$$ = \nabla_{\perp} \cdot (\Sigma(\mathbf{x}) \nabla_{\perp} \rho_{1,s}(\mathbf{x}^n \|\mathbf{x}\|) \rho_{2,s}(\|\mathbf{x}\|)), \quad (D.94) $$

$$ = \rho_{2,s}(\|\mathbf{x}\|) \|\mathbf{x}\|^2 \nabla_{\perp} \cdot (\Sigma(\mathbf{x}^n) \nabla_{\perp} \rho_{1,s}(\mathbf{x}^n \|\mathbf{x}\|)), \quad (D.95) $$

$$ = \rho_{2,s}(\|\mathbf{x}\|) \nabla_{\mathbb{S}^{d-1}} \cdot (\Sigma(\mathbf{x}^n) \nabla_{\mathbb{S}^{d-1}} \rho_{1,s}(\mathbf{x}^n \|\mathbf{x}\|)), \quad (D.96) $$

and finally,

$$ \rho_{2,s}(\|\mathbf{x}\|) \left( 2 \frac{\partial}{\partial s} \rho_{1,s}(\mathbf{x}^n \|\mathbf{x}\|) - \nabla_{\mathbb{S}^{d-1}} \cdot (\Sigma(\mathbf{x}^n) \nabla_{\mathbb{S}^{d-1}} \rho_{1,s}(\mathbf{x}^n \|\mathbf{x}\|)) \right) = 0. \quad (D.97) $$

Then, $\rho_{1,s}(\mathbf{x}^n \|\mathbf{x}\|)$ is a solution of the de Fokker-Planck equation on the sphere, for any $\mathbf{x}$ such that $\|\mathbf{x}\| \in A := \{r \in \mathbb{R}^+ | \rho_{2,s}(r) > 0\}$. If $\|\mathbf{x}\| \notin A$, then $\rho_{2,s}(\|\mathbf{x}\|) = 0$ and

$$ \rho_s(\mathbf{x}) = p_s^n(\mathbf{x}^n) \|\vec{\mathbf{x}}_s\| = \|\mathbf{x}\| p_{\cdot|\cdot}(\|\mathbf{x}\|) \|\mathbf{x}\|^{1-d}. $$

For $\|\mathbf{x}\| \in A$, then $\rho_{2,s}(\|\mathbf{x}\|) \neq 0$ and $\rho_{1,s}(\cdot \|\mathbf{x}\|)$ is solution of the Fokker Planck equation D.54. According to lemma D.4.1, the Fokker Planck equation D.54 has a unique density solution $p_s^n(\mathbf{x}^n \|\vec{\mathbf{x}}_0\| = \|\mathbf{x}\|)$. Hence, for any $\mathbf{x}$ such that $p_{\cdot|\cdot}(\|\mathbf{x}\|) > 0$, we have

$$ \rho_s(\mathbf{x}) = p_s^n(\mathbf{x}^n) \|\vec{\mathbf{x}}_s\| = \|\mathbf{x}\| p_{\cdot|\cdot}(\|\mathbf{x}\|) \|\mathbf{x}\|^{1-d}. $$
It is true for any $\boldsymbol{x}$. So, $\rho_s(\boldsymbol{x})$ is the unique solution of the Fokker Planck in $\mathbb{R}^d$.

*Regularity:* By definition of the marginal density, we have

$$p_{|\cdot|}(r) := \int_{\mathbb{S}^{d-1}} p_0(\Phi(r, \boldsymbol{\theta})) r^{d-1} d\boldsymbol{\theta} = r^{d-1} \int_{\mathbb{S}^{d-1}} p_0(\Phi(r, \boldsymbol{\theta})) d\boldsymbol{\theta}.$$

with $\Phi(r, \boldsymbol{\theta}) = r\boldsymbol{\theta}$. According to the assumption $\vec{\boldsymbol{x}}_0 \sim p_0 \in C^2(\mathbb{R}^d \setminus \{0\})$ and compactness of $\mathbb{S}^{d-1}$, one can conclude that $p_{|\cdot|} \in C^2([0, \infty[)$.

Since $p_0(\boldsymbol{x}^n \| \|\vec{\boldsymbol{x}}_0\| = r)$ is $C^2$, we have that $p_s^n(\boldsymbol{x}^n \mid |\vec{\boldsymbol{x}}_s| = r)$ is smooth by Lemma D.4.1 for any $s > 0$. Consequently, $\rho_s(\boldsymbol{x})$ is smooth on $D$ for any $s > 0$.

*Invariant distribution:* The distribution

$$p_{\infty}(\boldsymbol{x}) = \frac{p_{|\cdot|}(\|\boldsymbol{x}\|)}{\|\boldsymbol{x}\|^{d-1} |\mathbb{S}^{d-1}|}.$$

is radial function in $\|\boldsymbol{x}\|$. The operator $\nabla_{\perp}$ does not act on radial functions and $\frac{1}{|\mathbb{S}^{d-1}|}$ is in the kernel of $\nabla_{\perp}$ such that $\nabla_{\perp}(\frac{1}{|\mathbb{S}^{d-1}|}) = 0$. Hence

$$\frac{\partial}{\partial s} p_{\infty}(\boldsymbol{x}) = \nabla_{\perp} \cdot \left( \frac{1}{2} \boldsymbol{\Sigma}(\boldsymbol{x}) \nabla_{\perp} p_{\infty}(\boldsymbol{x}) \right), \quad \boldsymbol{x} \in D, \qquad (D.98)$$

$$= \nabla_{\perp} \cdot \left( \frac{1}{2} \boldsymbol{\Sigma}(\boldsymbol{x}) \nabla_{\perp} \left( \frac{p_{|\cdot|}(\|\boldsymbol{x}\|)}{\|\boldsymbol{x}\|^{d-1} |\mathbb{S}^{d-1}|} \right) \right), \qquad (D.99)$$

$$= \nabla_{\perp} \cdot \left( \frac{1}{2} p_{|\cdot|} (\|\boldsymbol{x}\|) \boldsymbol{\Sigma}(\boldsymbol{x}) \nabla_{\perp} \left( \frac{1}{|\mathbb{S}^{d-1}|} \right) \right), \qquad (D.100)$$

$$= 0. \qquad (D.101)$$

Therefore, the Fokker-Planck distribution $p_{\infty}$ is the stationary. In addition, $p_{\infty}$ is a density since

$$\int_{\mathbb{R}^d} p_{\infty}(\boldsymbol{x}) d\boldsymbol{x} = \int_{\mathbb{R}^d} \frac{p_{|\cdot|}(\|\boldsymbol{x}\|)}{\|\boldsymbol{x}\|^{d-1} |\mathbb{S}^{d-1}|} d\boldsymbol{x}, \qquad (D.102)$$

$$= \int_{\mathbb{R}_+} \int_{\mathbb{S}^{d-1}} \frac{p_{|\cdot|}(r)}{r^{d-1} |\mathbb{S}^{d-1}|} r^{d-1} dr d\boldsymbol{x}^n, \qquad (D.103)$$

$$= \int_{\mathbb{R}_+} \int_{\mathbb{S}^{d-1}} p_{|\cdot|}(r) \frac{1}{|\mathbb{S}^{d-1}|} dr d\boldsymbol{x}^n, \qquad (D.104)$$

$$= \int_{\mathbb{R}_+} p_{|\cdot|}(r) dr \int_{\mathbb{S}^{d-1}} \frac{1}{|\mathbb{S}^{d-1}|} d\boldsymbol{x}^n, \qquad (D.105)$$

$$= 1. \qquad (D.106)$$
Convergence: Hence, we obtain for $p_s = \rho_s$ that, we can bound the speed of convergence

$$ \|p_s - p_\infty\|_{L^2(\mathbb{R}^d)}^2 = \int_{\mathbb{R}^d} \left| p_s^n(\mathbf{x}^n \mid \vec{\mathbf{x}}_s) - \|\mathbf{x}\| p_{|\cdot|}(\|\mathbf{x}\|) \|\mathbf{x}\|^{1-d} - \frac{p_{|\cdot|}(\|\mathbf{x}\|)}{\|\mathbf{x}\|^{d-1}} \frac{1}{|\mathbb{S}^{d-1}|} \right|^2 d\mathbf{x}, \quad (D.107) $$

$$ = \int_{\mathbb{R}^d} \left| p_s^n(\mathbf{x}^n \mid \vec{\mathbf{x}}_s) - \|\mathbf{x}\| \right) - \frac{1}{|\mathbb{S}^{d-1}|} \left| \frac{p_{|\cdot|}(\|\mathbf{x}\|)^2}{\|\mathbf{x}\|^{2d-2}} \right| d\mathbf{x}, \quad (D.108) $$

$$ = \int_{\mathbb{R}_+} \left( \int_{\mathbb{S}^{d-1}} \left| p_s(\boldsymbol{\theta} \mid \vec{\boldsymbol{x}}_s) - \frac{1}{|\mathbb{S}^{d-1}|} \right|^2 d\boldsymbol{\theta} \right) \frac{p_{|\cdot|}(r)^2}{r^{d-1}} dr, \quad (D.109) $$

$$ = \int_{\mathbb{R}_+} \left( \left\| p_s(\cdot \mid \vec{\mathbf{x}}_s) - \frac{1}{|\mathbb{S}^{d-1}|} \right\|_{L^2(\mathbb{S}^{d-1})}^2 \right) \frac{p_{|\cdot|}(r)^2}{r^{d-1}} dr, \quad (D.110) $$

$$ \leq \exp(-\alpha s) \int_{\mathbb{R}_+} \left( \left\| p_0(\boldsymbol{\theta} \mid \vec{\mathbf{x}}_0) - \frac{1}{|\mathbb{S}^{d-1}|} \right\|_{L^2(\mathbb{S}^{d-1})}^2 \right) \frac{p_{|\cdot|}(r)^2}{r^{d-1}} dr, \quad (D.111) $$

$$ = \exp(-\alpha s) \int_{\mathbb{R}^d} \left| p_0(\mathbf{x}^n \| \vec{\mathbf{x}}_0 \| = \|\mathbf{x}\|) - \frac{1}{|\mathbb{S}^{d-1}|} \right|^2 \frac{p_{|\cdot|}(\|\mathbf{x}\|)^2}{\|\mathbf{x}\|^{2d-2}} d\mathbf{x}, \quad (D.112) $$

$$ = \exp(-\alpha s) \int_{\mathbb{R}^d} \left| p_0(\mathbf{x}^n \| \vec{\mathbf{x}}_0 \| = \|\mathbf{x}\|) p_{|\cdot|}(\|\mathbf{x}\|) \|\mathbf{x}\|^{1-d} - \frac{p_{|\cdot|}(\|\mathbf{x}\|)}{\|\mathbf{x}\|^{d-1}} \frac{1}{|\mathbb{S}^{d-1}|} \right|^2 d\mathbf{x}, \quad (D.113) $$

$$ = \exp(-\alpha s) \|p_0 - p_\infty\|_{L^2(\mathbb{R}^d)}^2, \quad (D.114) $$

where in the inequality we used Theorem D.4.2. The upper bound is finite since $p_\infty \in L^2(\mathbb{R}^d)$. In order to see this, we will show that the function $p_{|\cdot|}(r) r^{\frac{1-d}{2}}$ is in $L^2(0, \infty)$. Recall that

$$ p_{|\cdot|}(r) := \int_{\mathbb{S}^{d-1}} p_0(\Phi(r, \boldsymbol{\theta})) r^{d-1} d\boldsymbol{\theta} = r^{d-1} \int_{\mathbb{S}^{d-1}} p_0(\Phi(r, \boldsymbol{\theta})) d\boldsymbol{\theta}. $$

Then, application of Jensen inequality leads

$$ \int_{\mathbb{R}_+} p_{|\cdot|}(r)^2 r^{1-d} dr = \int_{\mathbb{R}_+} \left( r^{d-1} \int_{\mathbb{S}^{d-1}} p_0(\Phi(r, \boldsymbol{\theta})) d\boldsymbol{\theta} \right)^2 r^{1-d} dr, \quad (D.115) $$

$$ = \int_{\mathbb{R}_+} \left( \int_{\mathbb{S}^{d-1}} p_0(\Phi(r, \boldsymbol{\theta})) d\boldsymbol{\theta} \right)^2 r^{d-1} dr, \quad (D.116) $$

$$ \leq |\mathbb{S}^{d-1}| \int_{\mathbb{R}_+} \int_{\mathbb{S}^{d-1}} p_0(\Phi(r, \boldsymbol{\theta}))^2 r^{d-1} d\boldsymbol{\theta} dr, \quad (D.117) $$

$$ = |\mathbb{S}^{d-1}| \|p_0\|_{L^2(\mathbb{R}^d)}^2. \quad (D.118) $$

Consequently,

$$ \|p_\infty\|_{L^2(\mathbb{R}^d)}^2 = \left\| \frac{p_{|\cdot|}(\|\mathbf{x}\|)}{\|\mathbf{x}\|^{d-1}} \frac{1}{|\mathbb{S}^{d-1}|} \right\|_{L^2(\mathbb{R}^d)}^2, \quad (D.119) $$

$$ = \frac{1}{|\mathbb{S}^{d-1}|^2} \int_{\mathbb{S}^{d-1}} \int_{\mathbb{R}_+} \frac{p_{|\cdot|}(r)^2}{r^{d-1}} dr d\boldsymbol{\theta} \le \frac{1}{|\mathbb{S}^{d-1}|} \|p_0\|_{L^2(\mathbb{R}^d)}^2 < \infty. \quad (D.120) $$

□
## D.6 Beyond pure Stratonovich noise

A possible extension to the described diffusion in equation 3.1 would be to add a drift term, i.e. considering

$$d\vec{x}_s = A\vec{x}_s ds + G(\vec{x}_s) \circ d\vec{B}_s,$$

with a skew-symmetric matrix $A \in \mathbb{R}^{d,d}$. Then, the associated Fokker-Planck equations will additional involve an advection term $(Ax) \cdot \nabla_{\perp} p_s$, which can be used to improve the speed of convergence of the dynamics.

## E Latent distribution

The latent vectors $\vec{x}_{\infty} \sim p_{\infty}$ of additive SGM are Gaussian white noises. This is not the case for MSGM in general. This appendix will elaborate on this point. First, we will show that MSGM latent vectors are white noise in the weak sense. Then, we will discuss the conditions for these latent vectors to be Gaussian, how to sample them, and how to transform map them to another latent space which is Gaussian. We also show that the MSGM latent distribution is always closer than the SGM latent distribution to the data distribution. Finally, we focus on the case of Cauchy data distribution, where SGM leads to singularity, unlike MSGM.

### E.1 The invariant measures define white noises in the weak sense

In additive SGM, latent vectors $\vec{x}_{\infty} \sim p_{\infty}$ are Gaussian white noise in the strong sense, i.e. for any $i \neq j$, the coordinates $(\vec{x}_{\infty})_i$ and $(\vec{x}_{\infty})_j$ are centered, independent, and identically distributed. In contrast, the latent vectors of MSGM are white noises in the weak sense, as stated by the following proposition. For any $i \neq j$, the coordinates $(\vec{x}_{\infty})_i$ and $(\vec{x}_{\infty})_j$ are uncorrelated but neither Gaussian nor independent, in general.

**Proposition E.1.1.** Let the assumptions A1 and A2 hold, $E\|\vec{x}_{\infty}\|^2 < +\infty$, and $p_{\infty}$ be a stationary density of the Fokker-Planck equation D.22. Then, $\vec{x}_{\infty} \sim p_{\infty}$ is a white noise in the weak sense, i.e. $E\vec{x}_{\infty} = 0$, $E(\vec{x}_{\infty})_i^2 < +\infty$ independent of $i$, and $E(\vec{x}_{\infty})_i(\vec{x}_{\infty})_j = 0, \forall i, j \in \{1, \dots, d\}$ with $i \neq j$.

*Proof.* From Theorem D.2.1, $p_{\infty}$ is rotation-variant. So there exist a function $h: \mathbb{R}^{+} \to \mathbb{R}^{+}$ such that for any $x \in \mathbb{R}^d, p_{\infty}(x) = h(\|x\|)$. Then,

$$E(\vec{x}_{\infty})_i = \int_{\mathbb{R}^d} x_i h(\|x\|) dx = \int_{\mathbb{R}^{d-1}} \left( \int_{\mathbb{R}} x_i h(\|x\|) dx_i \right) \Pi_{k \neq i} dx_k = 0, \quad (E.1)$$

since the function $x_i \to x_i h(\|x\|)$ is even.

Similarly, for $i \neq j$ in $\{1, \dots, d\}$

$$E(\vec{x}_{\infty})_i(\vec{x}_{\infty})_j = \int_{\mathbb{R}^d} x_i x_j h(\|x\|) dx, \quad (E.2)$$

$$= \int_{\mathbb{R}^{d-1}} x_j \underbrace{\left( \int_{\mathbb{R}} x_i h(\|x\|) dx_i \right)}_{=0} \Pi_{k \neq i} dx_k, \quad (E.3)$$

$$= 0. \quad (E.4)$$

Moreover, for $i$ in $\{1, \dots, d\}$, we have

$$+\infty > E\|\vec{x}_{\infty}\|^2, \quad (E.5)$$

$$= E \sum_{i=1}^{d} (\vec{x}_{\infty})_i^2, \quad (E.6)$$

$$\geq E(\vec{x}_{\infty})_i^2, \quad (E.7)$$

$$= \int_{\mathbb{R}^{d-1}} \left( \int_{\mathbb{R}} x_i^2 h(\|x\|) dx_i \right) \Pi_{k \neq i} dx_k, \quad (E.8)$$
which does not depends of i.

**Remark 1.** Since $\mathbb{E}(\vec{\mathbf{x}}_\infty)^2_i$ does not depend on i, we can easily evaluate it from Theorem 3.3.1 and Proposition D.3.2

$$ \mathbb{E}(\vec{\mathbf{x}}_\infty)^2_i = \frac{1}{d} \mathbb{E}\|\vec{\mathbf{x}}_\infty\|^2 = \lim_{s \to +\infty} \frac{1}{d} \mathbb{E}\|\vec{\mathbf{x}}_s\|^2 = \frac{1}{d} \mathbb{E}\|\vec{\mathbf{x}}_0\|^2, \quad (E.9) $$

and thus

$$ \mathbb{E}(\vec{\mathbf{x}}_\infty \vec{\mathbf{x}}_\infty^\top) = \frac{1}{d} \mathbb{E}\|\vec{\mathbf{x}}_0\|^2 \mathbf{I}_d. \quad (E.10) $$

Therefore, monitoring the covariance of $\vec{\mathbf{x}}_T$ and its distance to $\frac{1}{d}\mathbb{E}\|\vec{\mathbf{x}}_0\|^2 \mathbf{I}_d$ are convenient proxies of the forward SDE convergence.

## E.2 Condition of Gaussianity for the latent vector

**Proposition E.2.1.** Let assumptions A1 and A2 hold and $p_{\parallel 1,2}$ be the density of $\|\vec{\mathbf{x}}_0\|^2$. Then, the latent distribution $p_\infty$ is Gaussian if and only if $p_{\parallel 1,2}$ is a scaled $\chi^2$ distribution with d degrees of freedom, denoted $\alpha^2\chi_d^2$ with $\alpha \ge 0$.

*Proof.* From Theorem D.2.1, we know that $p_\infty$ is rotation invariant, i.e. it is a function of $\|\mathbf{x}\|$. If this distribution is Gaussian, it has to be of the form $\mathcal{N}(0, \alpha^2 \mathbf{I}_d)$ with $\alpha \ge 0$. Then, $\|\vec{\mathbf{x}}_0\|^2 = \|x_\infty\|^2 \sim \alpha^2 \chi_d^2$. Reciprocally, if there exists $\alpha \ge 0$ such that $p_{\parallel 1,2} = \alpha^2 \chi_d^2$ then $p_{\parallel 1,2} = \alpha \chi_d$, where we denote by $\alpha \chi_d$ the distribution of a positive random variable $X = \sqrt{\alpha^2 R}$ such that $R \sim \chi_d^2$. From Theorem 3.3.1, we know that

$$ p_{\infty}(\mathbf{x}) = p_{\parallel 1,2}(\|\mathbf{x}\|) \frac{\|\mathbf{x}\|^{1-d}}{|\mathbb{S}^{d-1}|} = p_{\alpha \chi_d}(\|\mathbf{x}\|) \frac{\|\mathbf{x}\|^{1-d}}{|\mathbb{S}^{d-1}|}. \quad (E.11) $$

It is the distribution $\mathcal{N}(0, \alpha^2 \mathbf{I}_d)$ written in spherical form. So, the latent distribution $p_\infty$ is Gaussian.

**Remark 2.** Isotropic Gaussian data $\vec{\mathbf{x}}_0 \sim \mathcal{N}(0, \alpha^2 \mathbf{I}_d)$ will hence leads to Gaussian latent space. But the contrapositive is not true. To see this, let us consider a general spherical decomposition of the data distribution $p_0$:

$$ p_0(\mathbf{x}) = p^{\otimes} \left( \|\mathbf{x}\|, \frac{\mathbf{x}}{\|\mathbf{x}\|} \right) \|\mathbf{x}\|^{1-d} = p_{\parallel 1,2}(\|\mathbf{x}\|) p_0^n \left( \frac{\mathbf{x}}{\|\mathbf{x}\|} \middle| \|\mathbf{x}\| \right) \|\mathbf{x}\|^{1-d}. \quad (E.12) $$

The latent distribution would be Gaussian as long as the distribution of norm is $p_{\parallel 1,2} = p_{\alpha \chi_d}$. But the conditional distribution of direction can be any valid conditional distribution on the d-sphere. For instance,

$$ p_0^n \left( \frac{\mathbf{x}}{\|\mathbf{x}\|} \middle| \|\mathbf{x}\| \right) = \delta \left( \frac{\mathbf{x}}{\|\mathbf{x}\|} - e^{(1)} \right), \quad \text{with} \quad e^{(1)} = (1, 0, \dots, 0), \quad (E.13) $$

is a valid candidate even though $p_0$ is not Gaussian (since its support is $\mathbb{R}^+ \times \{0\}^{d-1}$).

## E.3 A tractable algorithm to sample latent vectors

With the following proposition, if we know the distribution of norms, $p_{\parallel 1,2}$, we can sample latent vectors from $p_\infty$.

**Proposition E.3.1.** Let $\vec{\mathbf{x}}_0 \sim \mathcal{N}(0, \mathbf{I}_d)$ and $\vec{\mathbf{x}}_0 = F^{-1}(\vec{\mathbf{x}}_0^N)$ with $F^{-1}(\mathbf{x}) := f(\|\mathbf{x}\|) \frac{\mathbf{x}}{\|\mathbf{x}\|}$ if $\mathbf{x} \ne 0$ and $0$ otherwise,

$$ f(r) := F_{\parallel 1,2}^{-1}(F_{\chi^2(d)}(r^2)), \quad \forall r > 0, \quad (E.14) $$

for the generalized inverse CDF $F_{\parallel 1,2}^{-1}$ of $p_{\parallel 1,2}$ and $F_{\chi^2}$ is the CDF of the $\chi^2$ distribution with d degrees of freedom. Then $\vec{\mathbf{x}}_0 \sim p_\infty$.
**Proof.** Since $\vec{\hat{x}}_0^N \sim \mathcal{N}(0, I_d)$, we know that $\|\vec{\hat{x}}_0^N\|^2 \sim \chi_{d-1}^2$, i.e. $F_{\chi^2}(\|\vec{\hat{x}}_0^N\|^2) \sim \mathcal{U}(0, 1)$ and then $R := f(\|\vec{\hat{x}}_0^N\|) = F_{|\cdot|}^{-1}(F_{\chi^2}(\|\vec{\hat{x}}_0^N\|^2)) \sim p_{|\cdot|}$. In addition, the normalized vector is $\frac{\vec{\hat{x}}_0^N}{\|\vec{\hat{x}}_0^N\|} \sim \mathcal{U}(\mathbb{S}^{d-1})$. The norm $\|\vec{\hat{x}}_0^N\|$ and the direction $\frac{\vec{\hat{x}}_0^N}{\|\vec{\hat{x}}_0^N\|}$ are independent. Therefore, $R$ and $\frac{\vec{\hat{x}}_0^N}{\|\vec{\hat{x}}_0^N\|}$ are also independent. We can conclude that $\vec{\hat{x}}_0 = R \frac{\vec{\hat{x}}_0^N}{\|\vec{\hat{x}}_0^N\|}$ follows the correct distribution. $\square$

In practice, we do not exactly know the distribution of the data norm $p_{|\cdot|}$. So, we do not have access to $F_{|\cdot|}$ or $f$. Instead, we approximate the distribution of $\log \|\vec{\hat{x}}_0\|_{\epsilon}$ with $\|\vec{x}\|_{\epsilon} := \|\vec{x}\| + \epsilon$, denoted $p_{\log|\cdot|_{\epsilon}}$, by a model $\hat{p}_{\log|\cdot|_{\epsilon}}$, or equivalently $F_{\log|\cdot|_{\epsilon}}$ by a model $\hat{F}_{\log|\cdot|_{\epsilon}}$ (see Appendix C). We perform a similar sampling procedure for the latent vectors, replacing $F$ by our approximation. We obtain samples of an approximate latent distribution $\hat{p}_{\infty}$, as stated by Proposition E.3.2.

**Proposition E.3.2.** Let $\vec{\hat{x}}_0^N \sim \mathcal{N}(0, I_d)$ and $\vec{\hat{x}}_0 = \hat{F}^{-1}(\vec{\hat{x}}_0^N)$ with $\hat{F}^{-1}(\vec{x}) := \hat{f}(\|\vec{x}\|) \frac{\vec{x}}{\|\vec{x}\|}$ if $\vec{x} \neq 0$ and 0 otherwise,

$$ \hat{f}(r) := \exp \left( \hat{F}_{\log|\cdot|_{\epsilon}}^{-1} (F_{\chi^2(d)}(r^2)) \right) - \epsilon, \quad \forall r > 0, \quad (E.15) $$

for the generalized inverse of the approximated CDF $\hat{F}_{\log|\cdot|_{\epsilon}}^{-1}$ associated to the approximated PDF $\hat{p}_{\log|\cdot|_{\epsilon}}$, and $F_{\chi^2}$ is the CDF of the $\chi^2$ distribution with $d$ degrees of freedom. Then $\vec{\hat{x}}_0 \sim \hat{p}_{\infty}$, where $\hat{p}_{\infty}$ is the empirical approximation of $p_{\infty}$, that is $\hat{p}_{\infty}(\vec{x}) := \hat{p}_{\log|\cdot|_{\epsilon}}(\log \|\vec{x}\|_{\epsilon}) \frac{\|\vec{x}\|^{1-d}}{\|\mathbb{S}^{d-1}\|}$, $\forall \vec{x} \in \mathbb{R}^d$.

**Proof.** We can follow the same proof that for Proposition E.3.1 replacing $F_{|\cdot|}$, $f$, $p_{|\cdot|}$, and $p_{\infty}$ by $\hat{F}_{\log|\cdot|_{\epsilon}}$, $\hat{f}$, $\hat{p}_{\log|\cdot|_{\epsilon}}$, and $\hat{p}_{\infty}$ respectively. $\square$

## E.4 Gaussianization of the latent vectors

If needed, we can easily build a second latent space with standard Gaussian vectors. As stated by the following proposition, for any (non-Gaussian) latent vector $\vec{\hat{x}}_T$, we can create a Gaussian vector $\vec{\hat{x}}_T^N$

$$ \vec{\hat{x}}_T^N = R_T \vec{\hat{x}}_T^n, \quad \text{with} \quad \vec{\hat{x}}_T^n = \vec{\hat{x}}_T / \|\vec{\hat{x}}_T\|, \quad \text{and} \quad R_T = \hat{f}^{-1}(\|\vec{\hat{x}}_T\|). \quad (E.16) $$

If $\vec{\hat{x}}_T$ is zero, we just set $\vec{\hat{x}}_T^N$ to zero.

**Proposition E.4.1.** Let $\vec{\hat{x}}_T \sim \hat{p}_{\infty}$, where $\hat{p}_{\infty}$ is the empirical approximation of $p_{\infty}$, that is $\hat{p}_{\infty}(\vec{x}) := \hat{p}_{\log|\cdot|_{\epsilon}}(\log \|\vec{x}\|_{\epsilon}) \frac{\|\vec{x}\|^{1-d}}{\|\mathbb{S}^{d-1}\|}$, $\forall \vec{x} \in \mathbb{R}^d$, and $\vec{\hat{x}}_T^N = \hat{F}(\vec{\hat{x}}_T)$ with $\hat{F}(\vec{x}) := \hat{f}^{-1}(\|\vec{x}\|) \frac{\vec{x}}{\|\vec{x}\|}$ if $\vec{x} \neq 0$ and 0 otherwise,

$$ \hat{f}^{-1}(r) = \sqrt{(F_{\chi^2}^{-1}(\hat{F}_{\log|\cdot|_{\epsilon}}(r)))}, \quad \forall r > 0, \quad (E.17) $$

for the approximated CDF $\hat{F}_{\log|\cdot|_{\epsilon}}$ associated to the approximated PDF $\hat{p}_{\log|\cdot|_{\epsilon}}$, and $F_{\chi^2}$ is the CDF of the $\chi^2$ distribution with $d$ degrees of freedom. Then $\vec{\hat{x}}_T^N \sim \mathcal{N}(0, I_d)$.

**Proof.** We can follow the same proof that for Proposition E.3.1 replacing $F_{|\cdot|}^{-1}$, $f$, $\chi^2$, $\mathcal{N}(0, I_d)$, $p_{|\cdot|}$, and $p_{\infty}$ by $\hat{F}_{\log|\cdot|_{\epsilon}}$, $\hat{f}^{-1}$, $\hat{p}_{\log|\cdot|_{\epsilon}}$, $\hat{p}_{\infty}$, $\chi^2$, and $\mathcal{N}(0, I_d)$ respectively. $\square$

## E.5 A shorter distance between latent and data distribution

The following result states, that the latent space of MSGM is closer to the data distribution compared to the SGM latent distribution in KL-divergence.

**Proposition E.5.1.** Let the assumptions A1 and A2 hold, $p_{|\cdot|}^2$ be the density of $\|\vec{\hat{x}}_0\|^2$, $p_{\infty}$ and $p_{\infty}^N = \mathcal{N}(0, I_d)$ be the MSGM and the SGM latent distributions respectively, then

$$ D_{KL}(p_{\infty}\|p_0) \le D_{KL}(p_{\infty}^N\|p_0), \quad (E.18) $$

with equality if and only if $p_{|\cdot|}^2$ is a $\chi^2$ distribution with $d$ degrees of freedom.
Proof. We recall that the MSGM latent pdf is

$$p_{\infty}(\boldsymbol{x}) = \frac{p_{\parallel,\parallel}(\|\boldsymbol{x}\|)}{\|\boldsymbol{x}\|^{d-1}} \frac{1}{\mathbb{S}^{d-1}}. \qquad (E.19)$$

and the data distribution reads

$$p_0(\boldsymbol{x}) = \frac{p_{\parallel,\parallel}(\|\boldsymbol{x}\|)}{\|\boldsymbol{x}\|^{d-1}} p_0(\boldsymbol{x}^n \|\vec{\boldsymbol{x}}_0\| = \|\boldsymbol{x}\|). \qquad (E.20)$$

Let denotes $p_{\chi_d^2}$ the $\chi^2$ distribution with $d$ degrees of freedom

$$p_0^{\mathcal{LN}}(\boldsymbol{x}) = \frac{p_{\chi_d^2}(\|\boldsymbol{x}\|)}{\|\boldsymbol{x}\|^{d-1}} p_0(\boldsymbol{x}^n \|\vec{\boldsymbol{x}}_0\| = \|\boldsymbol{x}\|). \qquad (E.21)$$

It is the distribution of $\vec{\boldsymbol{x}}_0^{\mathcal{LN}} = F(\vec{\boldsymbol{x}}_0)$ with $F(\boldsymbol{x}) := f^{-1}(\|\boldsymbol{x}\|) \frac{\boldsymbol{x}}{\|\boldsymbol{x}\|}$ if $\boldsymbol{x} \neq 0$ and 0 otherwise, and

$$f^{-1}(r) = \sqrt{(F_{\chi^2}^{-1}(F_{\parallel,\parallel}(r)))}, \quad \forall r > 0, \qquad (E.22)$$

and $F_{\parallel,\parallel}(R) = \int_0^R p_{\parallel,\parallel}(r)dr$ the cdf associated to $p_{\parallel,\parallel}$.

We have

$$0 \le D_{KL}(p_0 \| p_0^{\mathcal{LN}}), \qquad (E.23)$$

$$= \int p_0(\boldsymbol{x}) \log \frac{p_0}{p_0^{\mathcal{LN}}} d\boldsymbol{x}, \qquad (E.24)$$

$$= \int p_0(\boldsymbol{x}) \log \frac{p_{\parallel,\parallel}(\|\boldsymbol{x}\|)}{p_{\chi_d^2}(\|\boldsymbol{x}\|)} d\boldsymbol{x}, \qquad (E.25)$$

$$= \int p_0(\boldsymbol{x}) \log \frac{p_{\parallel,\parallel}(\|\boldsymbol{x}\|)}{\|\boldsymbol{x}\|^{d-1} \|\mathbb{S}^{d-1}\|} \frac{\|\boldsymbol{x}\|^{d-1} \|\mathbb{S}^{d-1}\|}{p_{\chi_d^2}(\|\boldsymbol{x}\|)} d\boldsymbol{x}, \qquad (E.26)$$

$$= \int p_0(\boldsymbol{x}) \log \frac{p_{\infty}(\boldsymbol{x})}{p_{\infty}^{\mathcal{N}}(\boldsymbol{x})} d\boldsymbol{x}, \qquad (E.27)$$

$$= \int p_0(\boldsymbol{x}) \log \frac{p_0(\boldsymbol{x})}{p_{\infty}^{\mathcal{N}}(\boldsymbol{x})} \frac{p_{\infty}(\boldsymbol{x})}{p_0(\boldsymbol{x})} d\boldsymbol{x}, \qquad (E.28)$$

$$= \int p_0(\boldsymbol{x}) \log \frac{p_0(\boldsymbol{x})}{p_{\infty}^{\mathcal{N}}(\boldsymbol{x})} d\boldsymbol{x} - \int p_0(\boldsymbol{x}) \log \frac{p_0(\boldsymbol{x})}{p_{\infty}(\boldsymbol{x})} d\boldsymbol{x}, \qquad (E.29)$$

$$= D_{KL}(p_0 \| p_{\infty}^{\mathcal{N}}) - D_{KL}(p_0 \| p_{\infty}), \qquad (E.30)$$

with equality if and only if $p_0 = p_0^{\mathcal{LN}}$ i.e. $p_{\parallel,\parallel} = p_{\chi_d^2}$.

□

## E.6 Relevance of MSGM latent space for heavy-tail distributions.

This appendix provides an analysis of why the latent space of MSGM is better suited to heavy-tailed data distribution as compared to the latent space of SGM. This subsection can be viewed as an extension of Proposition E.5.1. In particular the derived inequality in Proposition E.5.1 becomes meaning less if both sides are not finite. However, as we will see for example of heavy tail distribution such as the (product) Cauchy distribution, this is not the case. To this end we will show in Appendix E.6.1 that we KL divergence of data distribution and SGM latent space is not finite and that it is finite for the data distribution and the MSGM latent space in Appendix E.6.2.

We note that the analysis can be extended to a broader class of heavy tailed distributions and more general SGM latent spaces such as general Gaussian distributions.

### E.6.1 Infinite KL Divergence between Cauchy distribution and standard Gaussian

Let $\phi(\boldsymbol{x}) = \frac{1}{\sqrt{2\pi}} e^{-\boldsymbol{x}^2/2}$ be the density of the standard Gaussian $N(0, I)$, and let $p_0(\boldsymbol{x}) = \prod_{i=1}^d \frac{1}{\pi(1+x_i^2)}$ be the product density of univariate Cauchy distributions. Then, the following holds.
**Lemma E.6.1.**

$$D_{\text{KL}}(p_0 \| \phi) = \infty.$$

*Proof.* Let $L > 1$ and define the set

$$M = \{ \boldsymbol{x} \in \mathbb{R}^d \mid x_1 \ge L, |x_j| \le 1, \quad j = 2, \dots, d \}.$$

Then for $\boldsymbol{x} \in M$ and $C := \frac{1}{\pi^{d/2} 2^{d-1}}$

$$p_0(\boldsymbol{x}) = \prod_{i=1}^{d} \frac{1}{\pi(1+x_i^2)} \ge \frac{1}{\pi^d} \cdot \frac{1}{1+x_1^2} \cdot \prod_{i=2}^{d} \frac{1}{1+1} = \frac{1}{\pi^d 2^{d-1}} \cdot \frac{1}{1+x_1^2} = C \frac{1}{1+x_1^2}.$$

Moreover, for any $\boldsymbol{x} \in \mathbb{R}^d$, it holds that

$$\phi(\boldsymbol{x}) = (2\pi)^{-d/2} e^{-(x_1^2 + \sum_{i=2}^{d} x_i^2)/2} \le (2\pi)^{-d/2} e^{-x_1^2/2}.$$

Consequently, for $L$ large enough,

$$p_0(\boldsymbol{x}) \log \frac{p_0(\boldsymbol{x})}{\phi(\boldsymbol{x})} \ge \frac{x_1^2}{2} + \mathcal{O}(\log x_1),$$

where $\mathcal{O}$ refers to Landau-symbol of big-O notation. Together, for $\boldsymbol{x} \in M$ and $L$ large enough, there exists $\underline{C} > 0$ such that

$$p_0(\boldsymbol{x}) \log \frac{p_0(\boldsymbol{x})}{\phi(\boldsymbol{x})} \ge \frac{C}{1+x_1^2} \frac{x_1^2}{4} \ge \underline{C} > 0.$$

Consequently,

$$\int_M p_0(\boldsymbol{x}) \log \frac{p_0(\boldsymbol{x})}{\phi(\boldsymbol{x})} d\boldsymbol{x} \ge \int_T^\infty \left( \int_{|x_j| \le 1, j \ge 2} \int_{\underline{C}} d x_2 \cdots d x_d \right) d x_1 = \infty.$$

By Lebesgue decomposition, we conclude $D_{\text{KL}}(p_0 \| \phi) = \infty$. $\square$

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

**E.6.2 Finite KL divergence between Cauchy distribution and its related $\rho_\infty$**

Let $d \ge 2$ and again consider the product of Cauchy densities

$$p_0(\boldsymbol{x}) = \prod_{i=1}^{d} \frac{1}{\pi(1+x_i^2)}.$$

Let $\boldsymbol{x}_0 \sim p_0$ and let $p_R$ be the density of $R = \|\vec{\boldsymbol{x}}_0\|$. Then, motivated by our latent space distribution equation 3.8, consider the density

$$p_{\infty}(\boldsymbol{x}) = \frac{p_R(\|\boldsymbol{x}\|)}{\|\boldsymbol{x}\|^{d-1} |\mathbb{S}^{d-1}|}, \quad \boldsymbol{x} \in \mathbb{R}^d \setminus \{0\}, \qquad (\text{E.31})$$

Then, the following holds

**Lemma E.6.2.**

$$D_{\text{KL}}(p_0 \| p_{\infty}) < \infty.$$

*Proof.* It holds that

$$\log \frac{p_0(\boldsymbol{x})}{p_{\infty}(\boldsymbol{x})} = \log p_0(\boldsymbol{x}) - \log p_R(\|\boldsymbol{x}\|) + (d-1) \log \|\boldsymbol{x}\| + \log |\mathbb{S}^{d-1}|.$$

Hence,

$$D_{\text{KL}}(p_0 \| p_{\infty}) = \mathbb{E}_p[\log p_0(\vec{\boldsymbol{x}}_0)] - \mathbb{E}_p[\log p_R(\|\vec{\boldsymbol{x}}_0\|)] + (d-1)\mathbb{E}_p[\log \|\vec{\boldsymbol{x}}_0\|] + \log |\mathbb{S}^{d-1}|,$$

where $\mathbb{E}_p$ denotes the expectation with respect to the probability measure $p_0 d\boldsymbol{x}$. We will show, that each term separately is finite. We start with the first term, followed by the third. The finiteness of the second term turns out to be a consequence of the finiteness of the second term.
• First term: It holds that

$$ \log p_0(\vec{x}_0) = -d \log \pi - \sum_{i=1}^{d} \log(1 + (\vec{x}_0)_i^2) $$

for $\vec{x}_0 \sim p_0$. Since coordinates of $(\vec{x}_0)_i$ are iid it is enough to check the marginal integrals for finiteness. In particular it holds that

$$ \int_{-\infty}^{\infty} \frac{\log(1 + x^2)}{\pi(1 + x^2)} dx = \log(4) < \infty. $$

Consequently

$$ |\mathbb{E}_p[\log p_0(\vec{x}_0)]| < \infty. \qquad (E.32) $$

• Third term: For the second term, let $R = \|\vec{x}_0\|$ and $M = \max_{i=1,...,d} |(\vec{x}_0)_i|$. Then, almost surely

$$ M \le R \le \sqrt{dM} \Rightarrow \log M \le \log R \le \log M + \frac{1}{2} \log d. $$

Consequently,

$$ |\mathbb{E}_p[\log R] - \mathbb{E}_p[\log M]| \le \frac{1}{2} \log d. $$

Thus if we show $\mathbb{E}_p[\log M] < \infty$, then $\mathbb{E}_p[\log R] < \infty$ as well, since both expectation only differ up to a finite factor. Using the CDF $F$ of $(\vec{x}_0)_1$ e.g. for $(\vec{x}_0)_1$ it holds that

$$ \mathbb{P}(\| (\vec{x}_0)_1 \| \le t) = F(t) - F(-t) = \frac{2}{\pi} \arctan t, \quad t \ge 0. $$

Consequently, since $(\vec{x}_0)_1, \dots, (\vec{x}_0)_d$ are iid, the CDF $F_M$ of $M$ satisfies

$$ F_M(t) = \mathbb{P}(M \le t) = \left(\frac{2}{\pi} \arctan t\right)^d, \quad t \ge 0. $$

Hence, the density $f_M$ of $M$ is given (for $d \ge 2$) as

$$ f_M(t) = \frac{\partial}{\partial t} F_M(t) = d \left(\frac{2}{\pi}\right)^d (\arctan t)^{d-1} \frac{1}{1+t^2}. $$

Now, by a integral splitting we find that

$$ \mathbb{E}_p[\log M] = \int_0^\infty \log t f_M(t) dt = \int_0^1 \log t f_M(t) dt + \int_1^\infty \log(t) f_M(t) dt. \qquad (E.33) $$

By noting that for $0 \le t$, $f_M(t) \le C t^{d-1}$ for some $C > 0$ and

$$ \int_0^1 f_M(t)(-\log(t))dt \le \int_0^1 t^{a-1}(-\log(t))dt = \frac{1}{a^2}, \quad a > 0, $$

the first integrant of equation E.33 is finite using $a = d$. For $t \ge 1$ $\arctan(t) \le \pi/2$ and hence $f_M(t) \le C' \frac{1}{1+t^2}$ for some $C' > 0$ and the second integrant of equation E.33 is finite since

$$ \int_1^\infty \frac{\log(t)}{1+t^2} dt = 1 < \infty. $$

It follows that $\mathbb{E}_p[\log M]$ is finite.
• Second term: Recall that

$$p_R(r) = \int_{\mathbb{S}^{d-1}} p_0(r\boldsymbol{\theta})r^{d-1}d\sigma(\boldsymbol{\theta}) = r^{d-1} \int_{\mathbb{S}^{d-1}} p_0(r\boldsymbol{\theta})d\sigma(\boldsymbol{\theta}).$$

Since for $\boldsymbol{x} = r\boldsymbol{\theta}$ with $\boldsymbol{\theta} = (\theta_1, \dots, \theta_d)$, with $\theta_i^2 \le 1$ and using the fact that

$$p_0(r\boldsymbol{\theta}) = \prod_{i=1}^{d} \frac{1}{\pi(1 + r^2 \theta_i^2)}$$

we conclude

$$p_0(r\boldsymbol{\theta}) \ge \prod_{i=1}^{d} \frac{1}{\pi(1 + r^2)} = \frac{1}{\pi^d(1 + r^2)^d}.$$

Therefore,

$$p_R(r) \ge r^{d-1} \frac{1}{\pi^d(1+r^2)^d} \cdot |\mathbb{S}^{d-1}| := C_d \frac{r^{d-1}}{(1+r^2)^d}.$$

Hence,

$$\log p_R(r) \ge \log C_d + (d-1) \log r - d \log(1+r^2),$$

which yields

$$\mathbb{E}_p[\log p_R(\|\vec{\boldsymbol{x}}_0\|)] \ge \log C_d + (d-1)\mathbb{E}_p[\log R] - d\mathbb{E}_p[\log(1+R^2)]. \quad (E.34)$$

For the third term in equation E.34 it holds that $R^2 = \sum_{i=1}^d (\vec{x}_0)_i^2$. Now for $M \le 1$ we have since $R \le \sqrt{d}M \le \sqrt{d}$ that $\log(1 + R^2) \le \log(1 + \sqrt{d})$ is independent of $R$. For $M \ge 1$, $\log(1 + R^2) \le \log(1 + dM^2) \le \log(dM^2 + dM^2) = \log(2d) + 2\log(M)$. Since we already showed that $\mathbb{E}_p[\log R]$ is finite, we conclude that the lower bound in equation E.34 is finite. For the upper bound, note that

$$p_0(r\boldsymbol{\theta}) \le \frac{1}{\pi^d}, \quad \forall r > 0.$$

Thus,

$$p_R(r) \le r^{d-1} \frac{1}{\pi^d} |\mathbb{S}^{d-1}| =: C_d r^{d-1}$$

for some $C_d > 0$. So

$$\log p_R(r) \le \log C_d + (d-1) \log r.$$

And finally,

$$\mathbb{E}_p[\log_R(R)] \le \log C_d + (d-1)\mathbb{E}_p[\log R] < \infty$$

$$\text{since } \mathbb{E}_p[\log R] < \infty$$

• Fourth term: Finite since volume of the finite dimensional hypersphere.

□

## F Backward diffusion

This section is devoted to the derivation of the reverse SDE and ODE of our proposed MSGM in Itô and Stratonovich form.

**Proposition F.1. (Backward SDE)** Let the skew-symmetry assumption A1 hold. Then, the Itô form of the reverse SDE associated to the forward SDE 3.1 is given by the SDE

$$d\vec{\boldsymbol{x}}_t = \frac{1}{2}(\nabla \cdot \boldsymbol{\Sigma})(\vec{\boldsymbol{x}}_t)dt + \mathbf{G}(\vec{\boldsymbol{x}}_t)\mathbf{G}(\vec{\boldsymbol{x}}_t)^{\top}\nabla \log p_{T-t}(\vec{\boldsymbol{x}}_t)dt + \mathbf{G}(\vec{\boldsymbol{x}}_t)d\vec{\boldsymbol{B}}_t. \quad (F.1)$$

In the Stratonovich form, it reads:

$$d\vec{\boldsymbol{x}}_t = \mathbf{G}(\vec{\boldsymbol{x}}_t) \left( \mathbf{G}(\vec{\boldsymbol{x}}_t)^{\top} \nabla \log p_{T-t}(\vec{\boldsymbol{x}}_t) dt + \circ d\vec{\boldsymbol{B}}_t \right). \quad (F.2)$$
*Proof.* From Anderson (1982); Song et al. (2021) and the Itô forward SDE (see Lemma D.1.2), we know that the Itô reverse SDE with negative ds writes

$$d\vec{\tilde{x}}_s = \frac{1}{2}(\nabla \cdot \boldsymbol{\Sigma})(\vec{\tilde{x}}_s)ds - (\nabla \cdot \boldsymbol{\Sigma})(\vec{\tilde{x}}_s)ds - \mathbf{G}(\vec{\tilde{x}}_s)\mathbf{G}(\vec{\tilde{x}}_s)^{\mathrm{T}}\nabla \log p_s(\vec{\tilde{x}}_s)ds \\ + \mathbf{G}(\vec{\tilde{x}}_s)d\vec{\tilde{B}}_s, \qquad (\mathrm{F.3})$$

$$= -\frac{1}{2}(\nabla \cdot \boldsymbol{\Sigma})(\vec{\tilde{x}}_s)ds - \mathbf{G}(\vec{\tilde{x}}_s)\mathbf{G}(\vec{\tilde{x}}_s)^{\mathrm{T}}\nabla \log p_s(\vec{\tilde{x}}_s)ds + \mathbf{G}(\vec{\tilde{x}}_s)d\vec{\tilde{B}}_s. \qquad (\mathrm{F.4})$$

Replacing the decreasing $s \in [0, T]$ by $s = T - t$ with increasing $t \in [0, T]$ and using another Brownian motion $\vec{\tilde{B}}$, we obtain the Itô backward SDE with positive $dt$

$$d\vec{\tilde{x}}_t = \frac{1}{2}(\nabla \cdot \boldsymbol{\Sigma})(\vec{\tilde{x}}_t)dt + \mathbf{G}(\vec{\tilde{x}}_t)\mathbf{G}(\vec{\tilde{x}}_t)^{\mathrm{T}}\nabla \log p_{T-t}(\vec{\tilde{x}}_t)dt + \mathbf{G}(\vec{\tilde{x}}_t)d\vec{\tilde{B}}_t. \qquad (\mathrm{F.5})$$

Then, Lemma D.1.1 and the standard Stratonovich-to-Itô formula (e.g. Kunita, 1997) yields the Stratonovich form of the backward SDE:

$$d\vec{\tilde{x}}_t = -\frac{1}{2}d\langle \mathbf{G}(\vec{\tilde{x}}_t), \vec{\tilde{B}}_t \rangle_t + \frac{1}{2}(\nabla \cdot \boldsymbol{\Sigma})(\vec{\tilde{x}}_t)dt + \mathbf{G}(\vec{\tilde{x}}_t)\mathbf{G}(\vec{\tilde{x}}_t)^{\mathrm{T}}\nabla \log p_{T-t}(\vec{\tilde{x}}_t)dt \\ + \mathbf{G}(\vec{\tilde{x}}_t) \circ d\vec{\tilde{B}}_t, \qquad (\mathrm{F.6})$$

$$= -\frac{1}{2}(\nabla \cdot \boldsymbol{\Sigma})(\vec{\tilde{x}}_t)dt + \frac{1}{2}(\nabla \cdot \boldsymbol{\Sigma})(\vec{\tilde{x}}_t)dt$$

$$+ \mathbf{G}(\vec{\tilde{x}}_t) \left( \mathbf{G}(\vec{\tilde{x}}_t)^{\mathrm{T}}\nabla \log p_{T-t}(\vec{\tilde{x}}_t)dt + \circ d\vec{\tilde{B}}_t \right), \qquad (\mathrm{F.7})$$

$$= \mathbf{G}(\vec{\tilde{x}}_t) \left( \mathbf{G}(\vec{\tilde{x}}_t)^{\mathrm{T}}\nabla \log p_{T-t}(\vec{\tilde{x}}_t)dt + \circ d\vec{\tilde{B}}_t \right). \qquad (\mathrm{F.8})$$

□

**Proposition F.2. (Backward probability flow ODE)** Let the skew-symmetry assumption A1 hold. Then, the reverse probability flow associated to the forward SDE 3.1 is given by the ODE

$$\frac{d\vec{\tilde{x}}_t}{dt} = \frac{1}{2}\mathbf{G}(\vec{\tilde{x}}_t)\mathbf{G}(\vec{\tilde{x}}_t)^{\mathrm{T}}\nabla \log p_{T-t}(\vec{\tilde{x}}_t). \qquad (\mathrm{F.9})$$

*Proof.* From Song et al. (2021) and the Itô forward SDE (see Lemma D.1.2), we know that the reverse probability flow writes with negative ds

$$d\vec{\tilde{x}}_s = \frac{1}{2}(\nabla \cdot \boldsymbol{\Sigma})(\vec{\tilde{x}}_s)ds - \frac{1}{2}(\nabla \cdot \boldsymbol{\Sigma})(\vec{\tilde{x}}_s)ds - \frac{1}{2}\mathbf{G}(\vec{\tilde{x}}_s)\mathbf{G}(\vec{\tilde{x}}_s)^{\mathrm{T}}\nabla \log p_s(\vec{\tilde{x}}_s)dt, \qquad (\mathrm{F.10})$$

$$= -\frac{1}{2}\mathbf{G}(\vec{\tilde{x}}_s)\mathbf{G}(\vec{\tilde{x}}_s)^{\mathrm{T}}\nabla \log p_s(\vec{\tilde{x}}_s)dt. \qquad (\mathrm{F.11})$$

Replacing the decreasing $s \in [0, T]$ by $s = T - t$ with increasing $t \in [0, T]$ and using another Brownian motion $\vec{\tilde{B}}$, we obtain the Itô backward SDE with positive $dt$

$$\frac{d\vec{\tilde{x}}_t}{dt} = \frac{1}{2}\mathbf{G}(\vec{\tilde{x}}_t)\mathbf{G}(\vec{\tilde{x}}_t)^{\mathrm{T}}\nabla \log p_{T-t}(\vec{\tilde{x}}_t). \qquad (\mathrm{F.12})$$

□

# G Proof of Theorem 3.4.1: equivalence between ELBO and score matching

This appendix derives a score-matching-based ELBO for MSGM training. In this work, we focus on the simple forward multiplicative SDE equation 3.1. Nevertheless, we here derive a slightly more general theorem, where we include a possibly non-zero Stratonovich drift $f_S$.

## G.1 Statement of the theorem

**Theorem G.1.1.** Let us consider the forward SDE

$$d\vec{\tilde{x}}_s = f_S(\vec{\tilde{x}}_s)ds + \mathbf{G}(\vec{\tilde{x}}_s) \circ d\vec{\tilde{B}}_s, \qquad (\mathrm{G.1})$$
where assumption A1 holds. Then, we have

$$p_0(\mathbf{x}|\boldsymbol{\theta}) \geq \mathcal{E}_{\infty}(\mathbf{x}|\boldsymbol{\theta}), \qquad (G.2)$$

with the following ELBO

$$\mathcal{E}_{\infty}(\mathbf{x}|\boldsymbol{\theta}) := \mathbb{E} \left[ \log p_T(\vec{\mathbf{x}}_T) |\vec{\mathbf{x}}_0 = \mathbf{x} \right] \qquad (G.3) \\ - \int_0^T \mathbb{E}_{\vec{\mathbf{x}}_s} \left[ \frac{1}{2} \|\mathbf{a}_{\theta}(\vec{\mathbf{x}}_s, s)\|^2 + \nabla \cdot (\mathbf{G}(\vec{\mathbf{x}}_s) \mathbf{a}_{\theta}(\vec{\mathbf{x}}_s, s)) - f_S(\vec{\mathbf{x}}_s) \Big|_{\vec{\mathbf{x}}_0 = \mathbf{x}} \right] ds.$$

*Proof.* Here, we review the work of Huang et al. (2021) on SGM and generalize some of their results to derive an ELBO and justify score matching for MSGM. Note that Benton et al. (2024); Ren et al. (2025) proposes a very general SGM framework with associated ELBO and score matching losses. The MSGM ELBO and thus the above theorem can be understood as a particular case of their work. The explicit dependence in $\boldsymbol{\theta}$ is omitted for readability.

## G.2 Notations correspondence

The forward and backward processes are denoted $Y_s$ and $X_t$ in Huang et al. (2021) and $\vec{\mathbf{x}}_s$ and $\vec{\mathbf{x}}_t$ in this paper. The forward Itō equation of Huang et al. (2021) is denoted:

$$dY_s = f(Y_s, s)ds + g(Y_s, s)d\hat{B}_s, \qquad (G.4)$$

Lemma D.1.1 gives the forward Itō equation of MSGM. It yields the following notation correspondence:

$$g(\mathbf{x}, s) = \mathbf{G}(\mathbf{x}), \qquad (G.5)$$

$$D(\mathbf{x}) = \frac{1}{2}g(\mathbf{x})g(\mathbf{x})^{\top} = \frac{1}{2}\mathbf{\Sigma}(\mathbf{x}), \qquad (G.6)$$

$$f(\mathbf{x}) = \frac{1}{2}(\nabla \cdot \mathbf{\Sigma})^{\top}(\mathbf{x}) + \mathbf{f}_S(\mathbf{x}). \qquad (G.7)$$

And the backward equation is:

$$d\vec{\mathbf{x}}_t = \mu(\vec{\mathbf{x}}_t, t)dt + \mathbf{G}(\vec{\mathbf{x}}_t, t)d\vec{\mathbf{B}}_t, \qquad (G.8)$$

with a drift

$$\mu(\mathbf{x}, t) = -f(\mathbf{x}) + 2(\nabla \cdot D)^{\top}(\mathbf{x}) + 2D(\mathbf{x})\nabla \log p_{T-t}(\mathbf{x}), \qquad (G.9)$$

$$= -\mathbf{f}_S(\mathbf{x}) + \frac{1}{2}(\nabla \cdot \mathbf{\Sigma})^{\top}(\mathbf{x}) + \mathbf{\Sigma}(\mathbf{x})\nabla \log p_{T-t}(\mathbf{x}), \qquad (G.10)$$

where we would arrive at the approximate backward SDE of Figure 1 if we replace $\nabla \log p_{T-t}(\vec{\mathbf{x}}_t)$ by $s_{\theta}(\vec{\mathbf{x}}_t, T-t)$ also parametrized as $\mathbf{a}_{\theta} = \mathbf{G}^{\top}s_{\theta}$. We note that in our case, $\mathbf{f}_S = 0$, the drift reads $\boldsymbol{\mu} = \frac{1}{2}(\nabla \cdot \mathbf{\Sigma})^{\top} + \mathbf{\Sigma}\nabla \log p_{T-t}$, and the SDE simplifies with Stratonovich notations equation 3.12.

## G.3 Marginal density from Feynman-Kac representation

The Appendix D of Huang et al. (2021) treats the general case of multiplicative noise. It states that

$$p_0(\mathbf{x}) = \mathbb{E} \left[ p_T(\vec{\mathbf{x}}_T) \exp \left( \int_0^T \left( -\nabla \cdot \boldsymbol{\mu}(\vec{\mathbf{x}}_s, T-s) + \nabla \cdot \frac{1}{2}(\nabla \cdot \mathbf{\Sigma})^{\top}(\vec{\mathbf{x}}_s, T-s) \right) ds \right) \Big|_{\vec{\mathbf{x}}_0 = \mathbf{x}} \right], \qquad (G.11)$$

$$= \mathbb{E} \left[ p_T(\vec{\mathbf{x}}_T) \exp \left( - \int_0^T \nabla \cdot \left( \boldsymbol{\mu} - \frac{1}{2}(\nabla \cdot \mathbf{\Sigma})^{\top} \right) (\vec{\mathbf{x}}_s, T-s) ds \right) \Big|_{\vec{\mathbf{x}}_0 = \mathbf{x}} \right], \qquad (G.12)$$

where

$$d\vec{\mathbf{x}}_s = -\tilde{\mu}(\vec{\mathbf{x}}_s, T-s)ds + \mathbf{G}(\vec{\mathbf{x}}_s, T-s)d\mathbf{B}'_s, \qquad (G.13)$$

$$\tilde{\mu}(\mathbf{x}, t) = \mathbf{\mu}(\mathbf{x}, t) - (\nabla \cdot \mathbf{\Sigma})^{\top}(\mathbf{x}), \qquad (G.14)$$

and $\mathbf{B}'_s$ is a Brownian motion.
**Remark 3.** In our case, $\tilde{\mu}(x, t) = \mu(x, t) - (\nabla \cdot \Sigma)^{\mathrm{T}}(x) = \Sigma(x)\nabla \log p_{T-t}(x) - \frac{1}{2}(\nabla \cdot \Sigma)^{\mathrm{T}}(x)$.

**Remark 4.** Note that $\tilde{\mu} - \mu = -(\nabla \cdot \Sigma)^{\mathrm{T}} = -(\nabla \cdot \Sigma)^{\mathrm{T}}$ is twice the Itō to Stratonovich correction of the backward SDE equation G.8 (see Lemma D.1.2). It is expected since this SDE can be reversed in time once written with Stratonovich notations equation 3.12 (Kunita, 1997). Then, changing back from Stratonovich to Itō notations but with a different sign in front of the drift, we obtain the forward SDE equation G.14 verified by $\vec{x}_s$ including twice the Itō to Stratonovich correction.

## G.4 Change of measure and Jensen's inequality

From the Feynman-Kac representation equation G.12 and Jensen's inequality, we obtain an ELBO as in Huang et al. (2021).

Let $(\Omega, \mathcal{F}, \mathbb{P})$ be the underlying probability space for which $\mathcal{B}'$ is a Brownian motion. Suppose $\mathbb{Q}$ is another probability measure on $(\Omega, \mathcal{F})$ equivalent to $\mathbb{P}$ (i.e., they have the same measure zero sets). We can hence apply the change-of-measure

$$p_0(\mathbf{x}) = \mathbb{E} \left[ \frac{d\mathbb{P}}{d\mathbb{Q}} p_T(\vec{\mathbf{x}}_T) \exp \left( -\int_0^T \nabla \cdot (\mu - \frac{1}{2}(\nabla \cdot \mathbf{\Sigma})^{\mathrm{T}})(\vec{\mathbf{x}}_s, T-s)ds \right) \middle| \vec{\mathbf{x}}_0 = \mathbf{x} \right] \quad (G.15)$$

Then, we apply Jensen's inequality:

$$\log p_0(\mathbf{x}) \ge \underbrace{\mathbb{E} \left[ \log \frac{d\mathbb{P}}{d\mathbb{Q}} + \log p_T(\vec{\mathbf{x}}_T) - \int_0^T \nabla \cdot (\mu - \frac{1}{2}(\nabla \cdot \mathbf{\Sigma})^{\mathrm{T}})(\vec{\mathbf{x}}_s, T-s)ds \middle| \vec{\mathbf{x}}_0 = \mathbf{x} \right]}_{=\mathcal{E}^{\infty}} \quad (G.16)$$

Compared to Huang et al. (2021), we have the additional term $-\frac{1}{2}(\nabla \cdot \mathbf{\Sigma})^{\mathrm{T}}$, that is, $-\frac{1}{2}(\nabla \cdot \mathbf{\Sigma})^{\mathrm{T}}$.

## G.5 Girsanov theorem

Huang et al. (2021) apply the Girsanov theorem to the following forward SDE equation (17) of Huang et al. (2021):

$$d\vec{x}_s = (-\mu + \mathbf{G}\mathbf{a})ds + \mathbf{G}d\hat{\mathbf{B}}_s, \quad (G.17)$$

since the Itō to Stratonovich correction $\frac{1}{2}(\mu - \tilde{\mu}) = \frac{1}{2}(\nabla \cdot \mathbf{\Sigma})^{\mathrm{T}}$ is zero in Huang et al. (2021). However, it is not the case in MSGM and here we use the Girsanov theorem to this forward SDE instead:

$$d\vec{x}_s = (-\tilde{\mu} + \mathbf{G}\mathbf{a})ds + \mathbf{G}d\hat{\mathbf{B}}_s. \quad (G.18)$$

The Girsanov theorem (Oksendal, 1998, Theorem 8.6.3) states the following. Let $\hat{\mathbf{B}}$ be an Itō process solving

$$d\hat{\mathbf{B}}_s = -\mathbf{a}(\omega, s)ds + d\mathcal{B}'_s, \quad (G.19)$$

for $\omega \in \Omega$ and $\hat{\mathbf{B}}_0 = 0$ where $\mathbf{a}$ satisfies the Novikov's condition. Then $\hat{\mathbf{B}}$ is a Brownian motion with respect to $\mathbb{Q}$ and:

$$\mathbb{E} \left[ \log \frac{d\mathbb{P}}{d\mathbb{Q}} \middle| \vec{\mathbf{x}}_0 = \mathbf{x} \right] = \mathbb{E} \left[ \int_0^T \mathbf{a}(\omega, s) \cdot d\mathcal{B}'_s - \frac{1}{2} \int_0^T \|\mathbf{a}(\omega, s)\|_2^2 ds \middle| \vec{\mathbf{x}}_0 = \mathbf{x} \right], \quad (G.20)$$

$$= -\frac{1}{2} \int_0^T \mathbb{E}_{\vec{\mathbf{x}}_s} \left[ \|\mathbf{a}(\omega, s)\|_2^2 \middle| \vec{\mathbf{x}}_0 = \mathbf{x} \right] ds, \quad (G.21)$$

since $T \mapsto \int_0^T \mathbf{a}(\omega, s) \cdot d\mathcal{B}'_s$ is a martingale and thus $\mathbb{E} \left[ \int_0^T \mathbf{a}(\omega, s) \cdot d\mathcal{B}'_s \right] = 0$ (Oksendal, 1998, Theorem 3.2.1).
## G.6 ELBO evaluation

Equation G.21 enable us to evaluate the ELBO $\mathcal{E}^{\infty}$ given by equation G.16. To evaluate the divergence term, we note that:

$$ (\mu - \frac{1}{2}(\nabla \cdot \mathbf{\Sigma})^{\mathrm{T}})(\mathbf{x}, T - s) = -f_{S}(\mathbf{x}) + \frac{1}{2}(\nabla \cdot \mathbf{\Sigma})^{\mathrm{T}}(\mathbf{x}) + \mathbf{\Sigma}(\mathbf{x})s_{\theta}(\mathbf{x}, s) - \frac{1}{2}(\nabla \cdot \mathbf{\Sigma})^{\mathrm{T}}(\mathbf{x}), \qquad (G.22) $$

$$ = -f_S(\mathbf{x}) + \mathbf{G}(\mathbf{x})a_{\theta}(\mathbf{x}, s). \qquad (G.23) $$

Then, the ELBO simplifies to:

$$ \mathcal{E}^{\infty}(\mathbf{x}) = \mathbb{E} \left[ \log \frac{d\mathbb{P}}{d\mathbb{Q}} \middle| \vec{\mathbf{x}}_0 = \mathbf{x} \right] + \mathbb{E} \left[ \log p_T(\vec{\mathbf{x}}_T) \middle| \vec{\mathbf{x}}_0 = \mathbf{x} \right] \qquad (G.24) $$

$$ + \int_{0}^{T} \mathbb{E}_{\vec{\mathbf{x}}_s} \left[ -\nabla \cdot \left( \mu - \frac{1}{2}(\nabla \cdot \mathbf{\Sigma})^{\mathrm{T}} \right) \middle| \vec{\mathbf{x}}_0 = \mathbf{x} \right] ds, $$

$$ = \mathbb{E} \left[ \log p_T(\vec{\mathbf{x}}_T) \middle| \vec{\mathbf{x}}_0 = \mathbf{x} \right] \qquad (G.25) $$

$$ - \int_{0}^{T} \mathbb{E}_{\vec{\mathbf{x}}_s} \left[ \frac{1}{2} \| a_{\theta}(\vec{\mathbf{x}}_s, s) \|_{2}^{2} + \nabla \cdot (\mathbf{G}(\vec{\mathbf{x}}_s) a_{\theta}(\vec{\mathbf{x}}_s, s) - f_S(\vec{\mathbf{x}}_s)) \middle| \vec{\mathbf{x}}_0 = \mathbf{x} \right] ds. $$

□

We recall that in our case, $f_S$ cancels out. The first term $\mathbb{E}[\log p_T(\vec{\mathbf{x}}_T)\middle|\vec{\mathbf{x}}_0 = \mathbf{x}]$ is a constant w.r.t. $\theta$. So, if when maximizing the ELBDO, this term has no effect on the optimization procedure. Therefore, even with our multiplicative noise, the minimization of the ELBO corresponds precisely to Implicit Score Matching (ISM), which is itself equivalent to Explicit Score Matching (ESM), Sliced Score Matching (SSM) and Denoising Score Matching (DSM) (Huang et al., 2021).

## G.7 From ELBO to our SSM loss

Here we show how to derive our practical SSM loss equation 3.14 from Theorem 3.4.1. We assume the skew-symmetry condition A1 and zero Stratonovich drift, i.e. $f_s = 0$. The theorem states the

$$ p_0(x_0|\theta) \geq \mathcal{E}_{\infty}(x_0|\theta) := C(x_0) - \mathcal{L}_{\infty}(x_0|\theta) \qquad (G.26) $$

with $C$ being a constant with respect to the parameters $\theta$ to be learned. More precisely,

$$ C(x_0) = \mathbb{E} \left[ \log p_T(\vec{\mathbf{x}}_T) \middle| \vec{\mathbf{x}}_0 = x_0 \right], \qquad (G.27) $$

$$ \mathcal{L}_{\infty}(x_0|\theta) = \int_{0}^{T} \mathbb{E}_{\vec{\mathbf{x}}_s} \left[ \frac{1}{2} \| a_{\theta}(\vec{\mathbf{x}}_s, s) \|^2 + \nabla \cdot (\mathbf{G}(\vec{\mathbf{x}}_s) a_{\theta}(\vec{\mathbf{x}}_s, s) - f_S(\vec{\mathbf{x}}_s)) \middle| \vec{\mathbf{x}}_0 = x_0 \right] ds. \qquad (G.28) $$

Then, we average over the data $x_0$ to obtain the following lower bound for the likelihood of the dataset:

$$ \mathbb{E}_{x_0} p_0(x_0|\theta) \geq \mathbb{E}_{x_0} \mathcal{E}_{\infty}(x_0|\theta) = \mathbb{E}_{x_0} C(x_0) - \mathbb{E}_{x_0} \mathcal{L}_{\infty}(x_0|\theta). \qquad (G.29) $$

Our objective is to find the neural network parameters $\theta$, that try to maximize the likelihood of the data set, $\mathbb{E}_{x_0} p_0(x_0|\theta)$. Since $\mathbb{E}_{x_0} C(x_0)$ is a constant with respect to $\theta$, we maximize $-\mathbb{E}_{x_0} \mathcal{L}_{\infty}(x_0|\theta)$. Let us explicit the two terms above with the Hutchinson trick, $\mathbb{E}_{\mathbf{v} \sim \mathrm{Rad}(d)}[\mathbf{v}\mathbf{v}^{\mathrm{T}}] = \mathbf{I}_d$ (Song et al.,
2020)

$$ \mathbb{E}_{x_0} C(\mathbf{x}_0) = \mathbb{E}_{x_T} [\log p_T(\vec{\mathbf{x}}_T)], \qquad (G.30) $$

$$ \mathbb{E}_{x_0} \mathcal{L}_{\infty}(x_0 | \theta) = \int_{0}^{T} \mathbb{E}_{\vec{\mathbf{x}}_s} \left[ \frac{1}{2} \| \mathbf{a}_{\theta}(\vec{\mathbf{x}}_s, s) \|^2 + \nabla \cdot (\mathbf{G}(\vec{\mathbf{x}}_s) \mathbf{a}_{\theta}(\vec{\mathbf{x}}_s, s) - f_S(\vec{\mathbf{x}}_s)) \right] ds, \qquad (G.31) $$

$$ = T \int_{0}^{T} \mathbb{E}_{\vec{\mathbf{x}}_s} \left[ \frac{1}{2} \| \mathbf{a}_{\theta}(\vec{\mathbf{x}}_s, s) \|^2 + \nabla \cdot (\mathbf{G}(\vec{\mathbf{x}}_s) \mathbf{a}_{\theta}(\vec{\mathbf{x}}_s, s)) \right] \frac{1}{T} ds, \qquad (G.32) $$

$$ = T \mathbb{E}_{s \sim \mathcal{U}[0,T]} \mathbb{E}_{\vec{\mathbf{x}}_s} \left[ \frac{1}{2} \| \mathbf{a}_{\theta}(\vec{\mathbf{x}}_s, s) \|^2 + \nabla \cdot (\mathbf{G}(\vec{\mathbf{x}}_s) \mathbf{a}_{\theta}(\vec{\mathbf{x}}_s, s)) \right], \qquad (G.33) $$

$$ = T \mathbb{E}_{s \sim \mathcal{U}[0,T]} \mathbb{E}_{\vec{\mathbf{x}}_s} \left[ \frac{1}{2} \| \mathbf{a}_{\theta}(\vec{\mathbf{x}}_s, s) \|^2 + \nabla \cdot (\mathbb{E}_{\mathbf{v} \sim \mathrm{Rad}(d)} [\mathbf{v} \mathbf{v}^{\top}] \mathbf{G}(\vec{\mathbf{x}}_s) \mathbf{a}_{\theta}(\vec{\mathbf{x}}_s, s)) \right], \qquad (G.34) $$

$$ = T \mathbb{E}_{s \sim \mathcal{U}[0,T]} \mathbb{E}_{\vec{\mathbf{x}}_s} \mathbb{E}_{\mathbf{v} \sim \mathrm{Rad}(d)} \left[ \frac{1}{2} \| \mathbf{a}_{\theta}(\vec{\mathbf{x}}_s, s) \|^2 + (\mathbf{v} \cdot \nabla)(\mathbf{G}(\vec{\mathbf{x}}_s) \mathbf{a}_{\theta}(\vec{\mathbf{x}}_s, s)) \cdot \mathbf{v} \right]. \qquad (G.35) $$

$$ = T \mathcal{L}_{\mathrm{SSM}}(\theta). \qquad (G.36) $$

Therefore, maximizing the ELBO, $\mathbb{E}_{x_0} \mathcal{E}_{\infty}(x_0 | \theta) = \mathbb{E}_{x_0} C(x_0) - \mathbb{E}_{x_0} \mathcal{L}_{\infty}(x_0 | \theta)$, is equivalent to minimizing our practical score-matching loss, $\mathcal{L}_{\mathrm{SSM}}(\theta)$.

## G.8 Remark on the score parametrization

Following Huang et al. (2021), we directly model $\mathbf{G}(\vec{\mathbf{x}}_t)^{\top} \nabla \log p_s(\mathbf{x})$ by a neural network $\mathbf{a}_{\theta}(\mathbf{x}, s)$. If needed, the projected score, $\nabla_{\perp} \log p_s$, can be retrieved directly from $\mathbf{a}_{\theta}$ as shown below. Note that the full score,

$$ \nabla \log p_s = \nabla_{\perp} \log p_s + (\mathbf{x}^n \cdot \nabla) \log p_s, \qquad (G.37) $$

involves a radial term, $(\mathbf{x}^n \cdot \nabla) \log p_s$ that cannot be directly estimated in MSGM.

**Proposition G.1.** We assume that assumptions A1 and A2 hold, and that we have an approximation, $\mathbf{a}_{\theta}$, of the scaled score and an orthonormal basis $\mathbf{u}_2(\mathbf{x}), \dots, \mathbf{u}_d(\mathbf{x})$ of $\mathbf{x}^{\perp}$, that we concatenate in $\mathcal{U}(\mathbf{x}) = [\mathbf{u}_2(\mathbf{x}), \dots, \mathbf{u}_d(\mathbf{x})] \in \mathbb{R}^{d \times (d-1)}$. Then,

$$ [\mathcal{U}^{\top}(\mathbf{x}) \mathbf{\Sigma}(\mathbf{x}) \mathcal{U}(\mathbf{x})]^{-1} \mathcal{U}^{\top}(\mathbf{x}) \mathbf{G}(\mathbf{x}) \mathbf{a}_{\theta}(\mathbf{x}, s). \qquad (G.38) $$

approximates the projected score

$$ \mathcal{U}^{\top}(\mathbf{x}) \nabla_{\perp} \log p_s(\mathbf{x}). \qquad (G.39) $$

*Proof.* Since $\mathbb{R}^d = \mathbb{R}\mathbf{x}^n \perp \mathbf{x}^{\perp}$, we have $\mathcal{I}_d = \mathbf{x}^n(\mathbf{x}^n)^{\top} + \mathcal{U}(\mathbf{x})\mathcal{U}^{\top}(\mathbf{x})$. Using $\mathbf{\Sigma}(\mathbf{x})\mathbf{x}^n = 0$, we obtain

$$ \mathcal{U}^{\top}(\mathbf{x}) \mathbf{G}(\mathbf{x}) \mathbf{a}_{\theta}(\mathbf{x}, s) \approx \mathcal{U}^{\top}(\mathbf{x}) \mathbf{G}(\mathbf{x}) \mathbf{G}(\mathbf{x})^{\top} \nabla \log p_s(\mathbf{x}), \qquad (G.40) $$

$$ = \mathcal{U}^{\top}(\mathbf{x}) \mathbf{\Sigma}(\mathbf{x}) [\mathbf{x}^n(\mathbf{x}^n)^{\top} + \mathcal{U}(\mathbf{x}) \mathcal{U}^{\top}(\mathbf{x})] \nabla \log p_s(\mathbf{x}), \qquad (G.41) $$

$$ = \mathcal{U}^{\top}(\mathbf{x}) \mathbf{\Sigma}(\mathbf{x}) \mathcal{U}(\mathbf{x}) \mathcal{U}^{\top}(\mathbf{x}) \nabla \log p_s(\mathbf{x}), \qquad (G.42) $$

$$ = \mathcal{U}^{\top}(\mathbf{x}) \mathbf{\Sigma}(\mathbf{x}) \mathcal{U}(\mathbf{x}) \mathcal{U}^{\top}(\mathbf{x}) \nabla_{\perp} \log p_s(\mathbf{x}). \qquad (G.43) $$

$\mathcal{U}^{\top}(\mathbf{x}) \mathbf{\Sigma}(\mathbf{x}) \mathcal{U}(\mathbf{x}) \in \mathbb{R}^{(d-1) \times (d-1)}$ is full rank, so

$$ \mathcal{U}^{\top}(\mathbf{x}) \nabla_{\perp} \log p_s(\mathbf{x}) \approx [\mathcal{U}^{\top}(\mathbf{x}) \mathbf{\Sigma}(\mathbf{x}) \mathcal{U}(\mathbf{x})]^{-1} \mathcal{U}^{\top}(\mathbf{x}) \mathbf{G}(\mathbf{x}) \mathbf{a}_{\theta}(\mathbf{x}, s). \qquad (G.44) $$

□

It is also possible to model the score, $\nabla \log p_s(\mathbf{x})$ directly by a neural network, $\mathbf{s}_{\theta}(\mathbf{x}, s)$ using the following score-matching loss:

$$ \mathcal{L}_{\mathrm{SSM}}(\theta) = \mathbb{E}_{s \sim \mathcal{U}[0,T]} \mathbb{E}_{\vec{\mathbf{x}}_s} \mathbb{E}_{\mathbf{v} \sim \mathrm{Rad}(d)} \left[ \frac{1}{2} \| \mathbf{a}_{\theta}(\vec{\mathbf{x}}_s, s) \|^2 + (\mathbf{v} \cdot \nabla)(\mathbf{G}(\vec{\mathbf{x}}_s) \mathbf{a}_{\theta}(\vec{\mathbf{x}}_s, s) - \mathbf{f}_S(\vec{\mathbf{x}}_s)) \cdot \mathbf{v} \right], \qquad (G.45) $$

$$ = \mathbb{E}_{s \sim \mathcal{U}[0,T]} \mathbb{E}_{\vec{\mathbf{x}}_s} \mathbb{E}_{\mathbf{v} \sim \mathrm{Rad}(d)} \left[ \frac{1}{2} \| \mathbf{s}_{\theta}(\vec{\mathbf{x}}_s, s) \|^2_{\mathbf{\Sigma}(\vec{\mathbf{x}}_s)} + (\mathbf{v} \cdot \nabla)(\mathbf{\Sigma}(\vec{\mathbf{x}}_s) \mathbf{s}_{\theta}(\vec{\mathbf{x}}_s, s) - \mathbf{f}_S(\vec{\mathbf{x}}_s)) \cdot \mathbf{v} \right] ds. \qquad (G.46) $$
However, for any $\alpha \in \mathbb{R}$, $\mathcal{L}_{\text{SSM}}(\theta)(s_{\theta}) = \mathcal{L}_{\text{SSM}}(\theta)(s_{\theta} + \alpha x^n)$, i.e. our loss function is insensible to the radial component of the score $(x^n \cdot \nabla) \log p_s$. Therefore, our MSGM framework does not provide estimation for the radial score $(x^n \cdot \nabla) \log p_s$. Moreover, the optimization problem parametrized by $s_{\theta}$ is ill-defined, and the loss should probably be regularized as follows:

$$ \mathcal{L}_{\text{SSM}}^{\text{reg}}(\theta) = \mathcal{L}_{\text{SSM}}(\theta) + \gamma \mathbb{E}_{s \sim \mathcal{U}[0,T]} \mathbb{E}_{\vec{x}_s} [(\boldsymbol{x} \cdot \boldsymbol{s}_{\theta})^2], \quad (G.47) $$

with $\gamma > 0$ large, says $\gamma = 10^6$.

## G.9 Girsanov theorem in the transport noise literature

Following the work done by Huang et al. (2021) for additive noise, we have relied on the Girsanov theorem (Oksendal, 1998) to prove the equivalence between score matching and ELBO maximization for MSGM. Girsanov theorem is widely used, we may cite here its recent uses in the transport noise literature. In a Bayesian context, Cotter et al. (2020a, 2023); González et al. (2025); Singh et al. (2025) introduce nudging in their particle filter. Also used with other type of noises, nudging biases the noise to make the solution closer to the observations. Similarly, in our case, the weighted score, $a_{\theta}(\vec{x}_t, T-t)$, biases the noise, $d\vec{B}_t/dt$, in our backward SDE to make its solution closer to the forward SDE solution (see equation 3.12). This noise change is the core of Girsanov theorem (see equation G.19). Resseguier (2023) also proposed to fit a parametric model for the transport noise by maximum likelihood estimation.

# H Comparison with diffusions on Riemannian manifolds

This appendix describes the similarities between MSGM on $\mathbb{R}^d$ and SGMs on manifolds. To introduce the subject, we first recall some theoretical elements related to Riemannian manifolds. The link with SGMs on manifolds also suggests a particular neural network architecture that we exploit in this work.

## H.1 Riemannian Manifolds and Differentiation

This section is devoted to a brief introduction to Riemannian manifolds and the associated differential calculus. For a more comprehensive discussion, we refer to Lee (2018). Let $\mathcal{M}$ be a smooth $n$-dimensional embedded submanifold of $\mathbb{R}^d$, where $n \le d$. For any $\boldsymbol{x} \in \mathcal{M}$ we denote by $T_{\boldsymbol{x}}\mathcal{M}$ the tangential (linear) space of $\mathcal{M}$ at $\boldsymbol{x}$. We denote by $g$ a Riemannian metric on $\mathcal{M}$, which assigns to each $\boldsymbol{x} \in \mathcal{M}$ an inner product

$$ g_x : T_{\boldsymbol{x}}\mathcal{M} \times T_{\boldsymbol{x}}\mathcal{M} \to \mathbb{R}. $$

In the case of a smooth embedded manifold in the Euclidean space, the induced metric is given by

$$ g_x(u, v) = \langle u, v \rangle_{\mathbb{R}^n}, \quad \forall u, v \in T_{\boldsymbol{x}}\mathcal{M}. $$

This makes $(\mathcal{M}, g)$ a Riemannian manifold. Let $\{e^{(1)}, \dots, e^{(n)}\}$ be an orthonormal basis of $T_{\boldsymbol{x}}\mathcal{M}$. Then, the orthogonal projection onto $T_{\boldsymbol{x}}\mathcal{M}$ is the linear operator $P_{\boldsymbol{x}} : \mathbb{R}^d \to T_{\boldsymbol{x}}\mathcal{M}$ that satisfies

$$ P_{\boldsymbol{x}}(v) = \arg \min_{w \in T_{\boldsymbol{x}}\mathcal{M}} \|v - w\|_{\mathbb{R}^d} = \sum_{i=1}^{n} \langle v, e^{(i)} \rangle e^{(i)}. $$

While the concept of Riemannian gradients can be derived for general manifolds, here we limit ourselves to the simpler presentation of embedded manifolds in the Euclidean space. In this setup, the Riemannian manifold can be defined as the classical gradient projected to the tangential space. In particular, for $f : \mathbb{R}^d \to \mathbb{R}$ smooth, its *Riemannian gradient* can be computed as

$$ \nabla_{\mathcal{M}} f(\boldsymbol{x}) = P_{\boldsymbol{x}}(\nabla f(\boldsymbol{x})), $$

where $\nabla f(\boldsymbol{x})$ is the Euclidean gradient. Furthermore, we want to define the Riemannian divergence in this framework. For a tangent vector field $f : \mathcal{M} \to \mathbb{R}^d$ with $f(\boldsymbol{x}) \in T_{\boldsymbol{x}}\mathcal{M}$, the *Riemannian divergence* is given as

$$ \operatorname{div}_{\mathcal{M}} f(\boldsymbol{x}) = \sum_{i=1}^{n} \langle \partial_{e^{(i)}} f(\boldsymbol{x}), e^{(i)} \rangle, $$
where $\partial_{e^{(i)}} f$ is the Euclidean directional derivative. Finally, the Laplace-Beltrami operator $\Delta_M$ can be defined as

$$ \Delta_M f = \mathrm{div}_M(\nabla_M f), $$

which generalizes the Laplacian to $M$.

In the special case that $M = r\mathbb{S}^{d-1}$, for a radius $r > 0$ then $n = d-1$ and

$$ T_x M = T_x r \mathbb{S}^{d-1} = \{v \in \mathbb{R}^d | \langle v, x \rangle = 0\}, $$

and $P_x(v) = v - \frac{1}{r^2}\langle v, x \rangle x$. Since $x^n = \frac{x}{r}$ we obtain $P_x(v) = (I - x^n(x^n)^\top)v$ and as a result

$$ \nabla_M f(x) = P_x(\nabla f(x)) = (I - x^n(x^n)^\top)\nabla f(x) = \nabla_\perp f(x). \quad (H.1) $$

Regarding the Riemannian divergence, we note that $x^n, e^{(1)}, \dots, e^{(n)}$ defines an orthonormal basis of $\mathbb{R}^d$. By Lemma D.2.1

$$ \nabla \cdot f(x) = (x^n \cdot \nabla)(x^n \cdot f(x)) + \nabla_\perp \cdot f(x). $$

For $f(x) \in T_x M$, we have $f(x) \cdot x^n = 0$. Thus:

$$ \nabla_\perp \cdot f(x) = \nabla \cdot f(x) - \underbrace{(x^n \cdot \nabla)(f(x) \cdot x^n)}_{=0} = \nabla \cdot f(x). $$

Differentiating the tangency condition $f(x) \cdot x^n = 0$ along $x^n$ leads

$$ 0 = \partial_{x^n} (f(x) \cdot x^n) = \langle \partial_{x^n} f(x), x^n \rangle + \langle f(x), \partial_{x^n} x^n \rangle. $$

Since $\partial_{x^n} x^n = 0$, we conclude that $\langle \partial_{x^n} f(x), x^n \rangle = 0$. Finally, expanding $\nabla \cdot f(x)$ in $x^n, e^{(1)}, \dots, e^{(n)}$ leads to

$$ \nabla_\perp \cdot f(x) = \nabla \cdot f(x) = \sum_{i=1}^{n} \langle \partial_{e^{(i)}} f(x), e^{(i)} \rangle + \underbrace{\langle \partial_{x^n} f(x), x^n \rangle}_{=0} = \mathrm{div}_M f(x). \quad (H.2) $$

In our setting $\mathrm{Im}(\Sigma(x)) = x^\perp = T_x M$. Hence, the right-hand side of the Fokker-Planck equation 3.4

$$ \mathrm{div}_M(\Sigma(x)\nabla_M f(x)) = \nabla_\perp \cdot (\Sigma(x)\nabla_\perp f(x)), \quad (H.3) $$

generalizes the notion of a divergence-form operator to the manifold setup.

## H.2 Conditional diffusions on scaled d-spheres

Several authors have recently developed SGM on Riemannian manifolds (De Bortoli et al., 2022; Huang et al., 2022; Benton et al., 2024) in order to generate data lying on a particular manifold. Clearly different, our goal is more classical: generating data in $\mathbb{R}^d$. However, each solution path of our forward and backward SDE lies on its scaled d-sphere $\|\vec{x}_0\|\mathbb{S}^{d-1}$. Clearly, d-spheres are particular cases of Riemannian manifolds and possibly the most studied. De Bortoli et al. (2022) describes diffusions of $\vec{x}^n$ in the d-sphere $\mathbb{S}^{d-1}$. The simplest one involves a Brownian motion on the d-sphere that converges to the uniform distribution on the d-sphere, $p_\infty^n$. Unfortunately, this appealing proposal does not directly apply to our framework: the Brownian motion on the d-sphere is not a solution of our forward SDE of $\vec{x}^n$. Indeed, in general, there exists $\vec{x}^n \in \mathbb{S}^{d-1}$ such that $\Sigma(x^n) = \sum_{k=1}^d (G^k x^n)(G^k x^n)^\top \neq I_{\mathbb{S}^{d-1}}$. So, the Fokker-Planck equation of De Bortoli et al. (2022),

$$ \frac{\partial}{\partial s} p^n(x^n) = \mathrm{div}_{\mathbb{S}^{d-1}}(\nabla_{\mathbb{S}^{d-1}} p^n(x^n)), \quad \forall x^n \in \mathbb{S}^{d-1}, \quad (H.4) $$

and our Fokker-Planck equation for the direction,

$$ \frac{\partial}{\partial s} p^n(x^n) = \mathrm{div}_{\mathbb{S}^{d-1}}(\Sigma(x^n)\nabla_{\mathbb{S}^{d-1}} p^n(x^n)), \quad \forall x^n \in \mathbb{S}^{d-1}, \quad (H.5) $$

do not coincide. However, the analyses from the SGM-on-manifold community on the finite-time distribution, its score, approximations, and score-matching losses choices could certainly facilitate the MSGM training process in the future.
In our case, the norm of solution being constant along path, we can write both the forward and the backward equations of the direction on the unit d-sphere from equation 3.1 and equation 3.12:

$$d\vec{\mathbf{x}}_t^n = \mathbf{G}(\vec{\mathbf{x}}_t^n) \circ d\vec{\mathbf{B}}_t, \qquad (H.6)$$

$$d\vec{\mathbf{x}}_t^n = \frac{1}{\|\vec{\mathbf{x}}_t^n\|} \mathbf{G}(\vec{\mathbf{x}}_t^n) \left( \mathbf{G}(\vec{\mathbf{x}}_t^n)^\top \nabla \log p_{T-t} \left( \|\vec{\mathbf{x}}_t^n\| \|\vec{\mathbf{x}}_t^n\| \right) dt + \circ d\vec{\mathbf{B}}_t \right), \qquad (H.7)$$

$$= \mathbf{G}(\vec{\mathbf{x}}_t^n) \left( \mathbf{G}(\vec{\mathbf{x}}_t^n)^\top \left( \|\vec{\mathbf{x}}_t^n\| \nabla_\perp \log p_{T-t} \left( \|\vec{\mathbf{x}}_t^n\| \|\vec{\mathbf{x}}_t^n\| \right) \right) dt + \circ d\vec{\mathbf{B}}_t \right). \qquad (H.8)$$

We note that $\|\mathbf{x}\|\nabla_\perp = \nabla_{S^{d-1}} = \partial_{\mathbf{x}^n}$ is the Riemannian gradient on the scaled d-sphere $\|\mathbf{x}\|S^{d-1}$. Therefore, using $p_s^\top$, the density of the couple of variables $(\|\vec{\mathbf{x}}_s^n\|, \|\vec{\mathbf{x}}_s^n\|) \in \mathbb{R}^+ \times S^{d-1}$,

$$\|\vec{\mathbf{x}}_t\|\nabla_\perp \log p_{T-t} \left( \|\vec{\mathbf{x}}_t^n\| \|\vec{\mathbf{x}}_t^n\| \right) = \frac{\partial}{\partial \mathbf{x}^n} \log p_{T-t} \left( \|\vec{\mathbf{x}}_t^n\| \|\vec{\mathbf{x}}_t^n\| \right) \qquad (H.9)$$

$$= \frac{\partial}{\partial \mathbf{x}^n} \log \left( p_{T-t}^{\top} \left( \|\vec{\mathbf{x}}_t^n\|, \|\vec{\mathbf{x}}_t^n\| \right) \|\vec{\mathbf{x}}_t^n\|^{1-d} \right) \qquad (H.10)$$

$$= \frac{\partial}{\partial \mathbf{x}^n} \log \left( p_{T-t}^n \left( \|\vec{\mathbf{x}}_t^n\| \|\vec{\mathbf{x}}_t^n\| \right) p_{|\cdot|} \left( \|\vec{\mathbf{x}}_t^n\| \|\vec{\mathbf{x}}_t^n\| \right) \|\vec{\mathbf{x}}_t^n\|^{1-d} \right), \qquad (H.11)$$

$$= \frac{\partial}{\partial \mathbf{x}^n} \log p_{T-t}^n \left( \|\vec{\mathbf{x}}_t^n\| \|\vec{\mathbf{x}}_t^n\| \right) \qquad (H.12)$$

$$= \frac{\partial}{\partial \mathbf{x}^n} \log p_{T-t}^n \left( \|\vec{\mathbf{x}}_t^n\| \|\vec{\mathbf{x}}_0^n\| \right) \qquad (H.13)$$

$$= \nabla_{S^{d-1}} \log p_{T-t}^n \left( \|\vec{\mathbf{x}}_t^n\| \|\vec{\mathbf{x}}_0^n\| \right) \qquad (H.14)$$

and finally

$$d\vec{\mathbf{x}}_t^n = \mathbf{G}(\vec{\mathbf{x}}_t^n) \left( \mathbf{G}(\vec{\mathbf{x}}_t^n)^\top \nabla_{S^{d-1}} \log p_{T-t}^n \left( \|\vec{\mathbf{x}}_t^n\| \|\vec{\mathbf{x}}_0^n\| \right) dt + \circ d\vec{\mathbf{B}}_t \right). \qquad (H.15)$$

In contrast, forward and backward SDEs of De Bortoli et al. (2022) read

$$d\vec{\mathbf{x}}_t^n = d\vec{\mathbf{B}}_t^{S^{d-1}}, \qquad (H.16)$$

$$d\vec{\mathbf{x}}_t^n = \nabla_{S^{d-1}} \log p_{T-t}^n \left( \vec{\mathbf{x}}_t^n \right) dt + d\vec{\mathbf{B}}_t^{S^{d-1}}, \qquad (H.17)$$

where $\vec{\mathbf{B}}_t$ and $\vec{\mathbf{B}}_t^{S^{d-1}}$ are Brownian motions on the d-sphere. They can be defined from Stroock's representation (Hsu, 2002, Example 3.3.2) as

$$d\vec{\mathbf{B}}_t^{S^{d-1}} = (\mathbf{I}_d - (\vec{\mathbf{x}}_t^n)(\vec{\mathbf{x}}_t^n)^\top) \circ d\vec{\mathbf{B}}_t, \qquad (H.18)$$

$$d\vec{\mathbf{B}}_t^{S^{d-1}} = (\mathbf{I}_d - (\vec{\mathbf{x}}_t^n)(\vec{\mathbf{x}}_t^n)^\top) \circ d\vec{\mathbf{B}}_t. \qquad (H.19)$$

The first main difference with MSGM is that the projection on the tangent plane, $(\mathbf{I}_d - (\mathbf{x}^n)(\mathbf{x}^n)^\top)$, (quadratic in $\mathbf{x}^n$) is replaced in our approach by $\mathbf{G}(\vec{\mathbf{x}}_t^n)$ (linear in $\mathbf{x}^n$). Accordingly the noise (conditional) covariance, $(\mathbf{I}_d - (\mathbf{x}^n)(\mathbf{x}^n)^\top)^2 = (\mathbf{I}_d - (\mathbf{x}^n)(\mathbf{x}^n)^\top)$ (projection property), is replaced by $\mathbf{G}(\vec{\mathbf{x}}_t^n)\mathbf{G}(\vec{\mathbf{x}}_t^n)^\top = \mathbf{\Sigma}(\vec{\mathbf{x}}_t^n)$. To make our diffusion coincide with equation H.16, we would have to consider

$$\mathbf{G}(\mathbf{x}) := \|\mathbf{x}\| (\mathbf{I}_d - \mathbf{x}^n(\mathbf{x}^n)^\top), \qquad (H.20)$$

which is Lipchitz continuous but nonlinear. As such, the noise covariance would be

$$\mathbf{\Sigma}(\mathbf{x}) = \|\mathbf{x}\|^2 \mathbf{I}_d - \mathbf{x}\mathbf{x}^\top. \qquad (H.21)$$

In general, we can hardly expect such a simple form from MSGM noise covariance. However, for the random tensor equation 6.1, we can show (see equation J.10) that:

$$2\mathbb{E}\mathbf{\Sigma}(\mathbf{x}) = \|\mathbf{x}\|^2 \mathbf{I}_d - \mathbf{x}\mathbf{x}^\top. \qquad (H.22)$$
In addition, our score involved in the backward SDE equation H.15 depends on the norm $||\vec{\tilde{x}}_t||$. The norm $||\vec{\tilde{x}}_t|| = ||\vec{\tilde{x}}_0||$ appears as a covariable – with prior distribution $p_{|\cdot|}$ – for the diffusion on the unit d-sphere. This is another major difference of our approach compared to SGM on manifolds. Besides, from this point of view, we can better understand how the direction and magnitude are re-coupled during MSGM generation. Along the reverse diffusion, the conditional score direction H.14 will focus along some orientations, counterbalancing the direction equiprobability of the latent space, i.e. reversing the "whitening" of the forward process. On different scaled d-sphere $||x_0||^{\mathbb{S}^{d-1}}$, the conditional score direction will be oriented differently, pushing along some orientations on some spheres and along other directions on spheres of larger radius. Accordingly, along the backward diffusion, the directions tend to align differently on different hyperspheres. The distribution of direction become more and more radius-dependent.

If data samples $\vec{\tilde{x}}_0$ are snapshots of a conservative dynamical system, all data points probably have the similar energy $E = ||\vec{\tilde{x}}_0||^2$, i.e. $\text{Var}(E)/\mathbb{E}[E]^2$ is small. All data points are on closed scaled d-spheres $\sqrt{E} \mathbb{S}^{d-1}$ and our approach becomes even closer to De Bortoli et al. (2022).

## H.3 Link with neural network architecture

Form equation H.14, we also note that

$$ \mathbf{G}(\vec{\tilde{x}}_t)^{\mathrm{T}} \nabla \log p_{T-t}(\vec{\tilde{x}}_t) = \mathbf{G}(\vec{\tilde{x}}_t^n)^{\mathrm{T}} \nabla_{\mathbb{S}^{d-1}} \log p_{T-t}^n(\vec{\tilde{x}}_t^n \mid \|\vec{\tilde{x}}_0\|), \quad (H.23) $$

justifying our neural network spherical architecture equation L.33

$$ \mathbf{G}(\vec{\tilde{x}}_t)^{\mathrm{T}} \nabla \log p_{T-t}(\vec{\tilde{x}}_t) \approx a_{\theta}(\vec{\tilde{x}}_t, T-t) = \tilde{a}_{\theta} \left( \frac{\|\vec{\tilde{x}}_t\|}{\|\vec{\tilde{x}}_t\|_{\epsilon}} \vec{\tilde{x}}_t^n, \log \|\vec{\tilde{x}}_t\|_{\epsilon}, T-t \right). \quad (H.24) $$

# I Analytic illustrations on simplified cases

## I.1 The two-dimensional case

We note here that for $d = 2$, we can find an analytic solution for our multiplicative forward SDE. Moreover, it corresponds to the Brownian motion on the circle.

Let us recall this forward SDE:

$$ d\vec{x}(t) = \mathbf{G}(\vec{x}(t)) \circ d\vec{B}_t = \sum_{k=1}^{K} \mathbf{G}^k \vec{x}(t) \circ d\vec{B}_t = \left( \sum_{k=1}^{K} \mathbf{G}^k \circ d\vec{B}_t \right)^k \vec{x}(t), \quad (I.1) $$

In dimension 2,

$$ d\vec{x}(t) = \alpha J \vec{x}(t) \circ d\vec{B}_t, \quad (I.2) $$

where $\vec{x} = \begin{pmatrix} \vec{x}_1 \\ \vec{x}_2 \end{pmatrix} \in \mathbb{R}^2$, $J = \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix}$ is the $\frac{\pi}{2}$-rotation.

$$ d\begin{pmatrix} \vec{x}_1 \\ \vec{x}_2 \end{pmatrix} = \alpha \begin{pmatrix} -\vec{x}_2 \\ \vec{x}_1 \end{pmatrix} \circ dB_t, \quad (I.3) $$

Then, in the complex plane, $\vec{x}^C = \vec{x}_1 + i\vec{x}_2 \in \mathbb{C}$, with $i = \sqrt{-1}$ and we have:

$$ d\vec{x}^C(t) = \alpha i \vec{x}^C \circ dB_t, \quad (I.4) $$

since

$$ d\vec{x}_1 + i d\vec{x}_2 = \alpha i (\vec{x}_1 + i\vec{x}_2) \circ dB_t = \alpha (-\vec{x}_2 + i\vec{x}_1) \circ dB_t. \quad (I.5) $$

The solution is the Brownian motion on the circle:

$$ \vec{x}^C(t) = \vec{x}^C(0) \exp(\alpha i B_t), \quad (I.6) $$
i.e.

$$ \vec{\mathbf{x}}(t) = R(\alpha \vec{B}_t) \vec{\mathbf{x}}(0) = \begin{pmatrix} \cos(\alpha \vec{B}_t) & -\sin(\alpha \vec{B}_t) \\ \sin(\alpha \vec{B}_t) & \cos(\alpha \vec{B}_t) \end{pmatrix} \vec{\mathbf{x}}(0), \quad (I.7) $$

i.e.

$$ \vec{x}_1(t) = \vec{x}_1(0) \cos(\alpha \vec{B}_t) - \vec{x}_2(0) \sin(\alpha \vec{B}_t), \quad (I.8) $$

$$ \vec{x}_2(t) = \vec{x}_1(0) \sin(\alpha \vec{B}_t) + \vec{x}_2(0) \cos(\alpha \vec{B}_t). \quad (I.9) $$

The key element of the proof was the possibility to write the forward diffusion with a single skew-symmetric matrix in equation I.2. Below we generalize this idea to larger dimension $d \ge 2$.

## I.2 Tensor built from a single skew-symmetric matrix

Here we assume that whole tensor $\mathbf{G}$ is built from the same dense skew-symmetric matrix $\mathbf{G}^1$ i.e.

$$ \mathbf{G}^k = \mathbf{G}^1, \quad \forall k \in \{1, \dots, d\}, \quad (I.10) $$

with $\mathbf{G}^1$ a skew-symmetric matrix. As explained in Appendix K.2.1, this tensor respect the condition A1 but not the A2. Nevertheless, this case and its analytic solution may be insightful.

### I.2.1 Matrix exponential

Here the full Brownian matrix $\mathbf{Z}$ can be simply factorized as

$$ \mathbf{Z}_s = \sum_{k=1}^{d} \mathbf{G}^k (\vec{B}_s)_k = \mathbf{G}^1 \sum_{k=1}^{d} (\vec{B}_s)_k. \quad (I.11) $$

It has the same distribution than

$$ \mathbf{Z}'_s = \sqrt{d} \mathbf{G}^1 \vec{B}_s', \quad (I.12) $$

with $\vec{B}'$ another single Brownian motion. The forward diffusion simplify to

$$ d\vec{\mathbf{x}}_s = \sqrt{d} \mathbf{G}^1 \vec{\mathbf{x}}_s \circ d\vec{B}_s', \quad (I.13) $$

with solution

$$ \vec{\mathbf{x}}_s = \exp(\mathbf{Z}'_s) \vec{\mathbf{x}}_0 = \exp\left(\sqrt{d} \mathbf{G}^1 \vec{B}_s'\right) \vec{\mathbf{x}}_0, \quad (I.14) $$

since $\mathbf{Z}'_s$ and $d\mathbf{Z}'_s$ commute.

### I.2.2 Diagonalization in the complex plane

$\mathbf{G}^1$ has pure imaginary eigenvalues and can be diagonalized in $\mathbb{C}$ on an orthonormal basis

$$ \mathbf{G}^1 = U_{\mathbb{C}}(i\mathbf{\Lambda})U_{\mathbb{C}}^H, \quad (I.15) $$

with $U_{\mathbb{C}}$ a complex unitary matrix, $\mathbf{\Lambda}$ a real diagonal matrix, and the superscript $H$ denotes the conjugate transpose. Then, the solution can be easily evaluate as follow

$$ \vec{\mathbf{x}}_s = U_{\mathbb{C}} \exp \left( i \sqrt{d} \mathbf{\Lambda} \vec{B}_s' \right) U_{\mathbb{C}}^H \vec{\mathbf{x}}_0, \quad (I.16) $$

For an even dimension $d$, and for all $j \in \{1, \dots, d/2\}$, there exists $\lambda_j \in \mathbb{R}$ such that

$$ (U_{\mathbb{C}}^H \vec{\mathbf{x}}_s)_{2j-1} = \exp \left( i \sqrt{d} \lambda_j \vec{B}_s' \right) (U_{\mathbb{C}}^H \vec{\mathbf{x}}_0)_{2j-1}, \quad (I.17) $$

$$ (U_{\mathbb{C}}^H \vec{\mathbf{x}}_s)_{2j} = \exp \left( -i \sqrt{d} \lambda_j \vec{B}_s' \right) (U_{\mathbb{C}}^H \vec{\mathbf{x}}_0)_{2j}. \quad (I.18) $$
For an odd dimension $d$, $G^1$ has at least one zero eigenvalue. Without loss of generality, we consider $\Lambda_{d,d} = 0$ and for all $j \in \{1, \dots, (d-1)/2\}$, there exists $\lambda_j \in \mathbb{R}$ such that

$$ (U_{\mathbb{C}}^{H} \vec{\mathbf{x}}_s)_{2j-1} = \exp \left( i \sqrt{d} \lambda_j \vec{B}_s' \right) (U_{\mathbb{C}}^{H} \vec{\mathbf{x}}_0)_{2j-1}, \quad (I.19) $$

$$ (U_{\mathbb{C}}^{H} \vec{\mathbf{x}}_s)_{2j} = \exp \left( -i \sqrt{d} \lambda_j \vec{B}_s' \right) (U_{\mathbb{C}}^{H} \vec{\mathbf{x}}_0)_{2j}, \quad (I.20) $$

$$ (U_{\mathbb{C}}^{H} \vec{\mathbf{x}}_s)_d = (U_{\mathbb{C}}^{H} \vec{\mathbf{x}}_0)_d. \quad (I.21) $$

### I.2.3 Real solution with sine and cosine

The diagonalization matrix, $U$, is complex but we can find a real unitary matrix, $U_{\mathbb{R}}$, to make $G_1$ block diagonal, and then expressing the solution with cosine and sinus as in equation I.8 and equation I.9:

$$ (U_{\mathbb{R}}^{\mathrm{T}} \vec{\mathbf{x}}_s)_{2j-1} = \cos \left( \sqrt{d} \lambda_j \vec{B}_s' \right) (U_{\mathbb{R}}^{\mathrm{T}} \vec{\mathbf{x}}_0)_{2j-1} - \sin \left( \sqrt{d} \lambda_j \vec{B}_s' \right) (U_{\mathbb{R}}^{\mathrm{T}} \vec{\mathbf{x}}_0)_{2j}, \quad (I.22) $$

$$ (U_{\mathbb{R}}^{\mathrm{T}} \vec{\mathbf{x}}_s)_{2j} = \sin \left( \sqrt{d} \lambda_j \vec{B}_s' \right) (U_{\mathbb{R}}^{\mathrm{T}} \vec{\mathbf{x}}_0)_{2j-1} + \cos \left( \sqrt{d} \lambda_j \vec{B}_s' \right) (U_{\mathbb{R}}^{\mathrm{T}} \vec{\mathbf{x}}_0)_{2j}, \quad (I.23) $$

For an odd dimension $d$, the real solution reads

$$ (U_{\mathbb{R}}^{\mathrm{T}} \vec{\mathbf{x}}_s)_{2j-1} = \cos \left( \sqrt{d} \lambda_j \vec{B}_s' \right) (U_{\mathbb{R}}^{\mathrm{T}} \vec{\mathbf{x}}_0)_{2j-1} - \sin \left( \sqrt{d} \lambda_j \vec{B}_s' \right) (U_{\mathbb{R}}^{\mathrm{T}} \vec{\mathbf{x}}_0)_{2j}, \quad (I.24) $$

$$ (U_{\mathbb{R}}^{\mathrm{T}} \vec{\mathbf{x}}_s)_{2j} = \sin \left( \sqrt{d} \lambda_j \vec{B}_s' \right) (U_{\mathbb{R}}^{\mathrm{T}} \vec{\mathbf{x}}_0)_{2j-1} + \cos \left( \sqrt{d} \lambda_j \vec{B}_s' \right) (U_{\mathbb{R}}^{\mathrm{T}} \vec{\mathbf{x}}_0)_{2j}, \quad (I.25) $$

$$ (U_{\mathbb{R}}^{\mathrm{T}} \vec{\mathbf{x}}_s)_d = (U_{\mathbb{R}}^{\mathrm{T}} \vec{\mathbf{x}}_0)_d. \quad (I.26) $$

Figure 5 illustrates the solution for $d = 4$ with 20000 realizations of $\vec{x}_T$ at large time $T = 100$, with $\lambda_1 = 1$, $\lambda_2 = 10$, $\vec{x}_0 = (1, 1, 1, 1)$, and $U_{\mathbb{R}} = I_4$. A rotation-invariant distribution, $p_{\infty}$, would induce rotation-invariant marginals and hence point cloud projections appearing rotation-invariant. This is clearly not the case here. This counter example shows that low-rank tensors as defined in equation I.10 cannot guaranty rotation-invariant latent distribution, and thus prevent the use of our simple eCDF-based sampling procedure.

Figure 7 illustrates the latent vector support for a random initial condition $\vec{x}_0 = N((1, 1, 1, 1), 0.01I_4)$. The supporting manifold is not one-dimensional anymore, but still depend on the initial direction distribution, $p_0^n$.

## I.3 Non-commutativity in the general case

For a general tensor $G$ in dimension $d > 2$, it is temping to look for a solution $\vec{x}_s$ of the forward SDE

$$ d\vec{x}_s = \circ dZ_s \vec{x}_s, $$

of the form $\exp(Z_s)\vec{x}_0$ with $Z = \sum_{k=1}^d G^k(\vec{B}_s)_k$. $Z_s$ being skew-symmetric, $\exp(Z_s)$ is unitary and such a solution would be reminiscent of the rotation form of equation I.6 and equation I.16 derived above. However, $\exp(Z_s)\vec{x}_0$ is not a solution of equation 3.1 in general, since $dZ_s Z_s \neq Z_s dZ_s$.

## J Rank and skew-symmetry conditions for random tensor $G$

In this appendix, we treat the case of random tensor $G$ as defined by equation 6.1. We will show that this tensor respects both assumptions A1 and A2 almost surely. Then, we will discuss the speed of contraction of the Fokker-Planck equation with this tensor.
<figure><img src="image_8.png" /><figcaption>Figure 5: Projection of samples, $\vec{x}_T$, sketching the support of the invariant measure, $p_\infty$, for a low-rank tensor I.10, $d = 4$ and $\vec{x}_0 = (1, 1, 1, 1)$. The top plot is in space $(x_1, x_2, x_3)$, the bottom plots are, from left to right, in space $(x_1, x_2)$, $(x_1, x_3)$, $(x_1, x_4)$, and $(x_4, x_3)$.</figcaption></figure>

## J.1 Proof of the rank condition

**Proposition J.1.** Let $M^k \in \mathbb{R}^{d,d}$ be iid random matrices with entries drawn independently from $N(0, 1)$. Define the skew-symmetric matrices $G^k = \frac{1}{2}(M^k - (M^k)^T)$ and for $x \in \mathbb{R}^d \setminus \{0\}$ define the (random) matrix

$$G(x) := [G^1 x, G^2 x, \dots, G^d x] \in \mathbb{R}^{d,d}.$$

Then, almost surely $\mathrm{rank}(G(x)) = d - 1$.

*Proof.* Let $x \neq 0$. Let $M$ be a random standard Gaussian matrix. Then, let $D = M - M^T$. Then, $D$ is Gaussian matrix with entries drawn from $N(0, 2)$, in particular

$$\mathbb{E}[D_{ij}D_{k\ell}] = \mathbb{E}[(M_{ij} - M_{ij})(M_{k\ell} - M_{\ell k})] = 2(\delta_{ik}\delta_{j\ell} - \delta_{i\ell}\delta_{jk}). \quad (J.1)$$

Consequently,

$$\mathbb{E}[(M - M^T)x x^T (M^T - M)] = -\mathbb{E}[D x x^T D]. \quad (J.2)$$
Moreover, as expected from the expression above, the initial norm, $||\vec{x}_0||$, scales the one-dimensional manifold supporting the invariant measure (not shown) and the initial direction, $\vec{x}_0^n$, has an influence at large time. Figure 6 shows the same example with $\vec{x}_0 = (\sqrt{2}, \sqrt{2}, 0, 0)$. The initial norm is the same but the initial direction is different. Therefore, the limit distribution, $p_\infty$, if it exists does depend on the initial direction, $\vec{x}_0^n$, making the latent sampling intractable.

<figure><img src="image_9.png" /><figcaption>Figure 6: Projection of samples, $\vec{x}_T$, sketching the support of the invariant measure, $p_\infty$, for a low-rank tensor I.10, $d = 4$ and $\vec{x}_0 = (\sqrt{2}, \sqrt{2}, 0, 0)$. The top plot is in space $(x_1, x_2, x_3)$, the bottom plots are, form left to right, in space $(x_1, x_2)$, $(x_1, x_3)$, $(x_1, x_4)$, and $(x_4, x_3)$.</figcaption></figure>

Now, for the covariance structure it holds

$$-(\mathbb{E}[D\boldsymbol{x}\boldsymbol{x}^\top D])_{ik} = \mathbb{E}[(D\boldsymbol{x})_i(D\boldsymbol{x})_k], \qquad (J.3)$$

$$= \sum_{j=1}^{d} \sum_{\ell=1}^{d} \mathbb{E}[D_{ij}x_j D_{k\ell}x_\ell], \qquad (J.4)$$

$$= \sum_{j=1}^{d} \sum_{\ell=1}^{d} x_j x_\ell \mathbb{E}[D_{ij} D_{k\ell}], \qquad (J.5)$$

$$= 2 \sum_{j=1}^{d} \sum_{\ell=1}^{d} x_j x_\ell (\delta_{ik} \delta_{j\ell} - \delta_{i\ell} \delta_{jk}), \qquad (J.6)$$

$$= 2\delta_{ik} \sum_{j=1}^{d} \sum_{\ell=1}^{d} x_j x_\ell \delta_{j\ell} - 2 \sum_{j=1}^{d} \sum_{\ell=1}^{d} x_j x_\ell \delta_{i\ell} \delta_{jk}, \qquad (J.7)$$

$$= 2\delta_{ik} \|\boldsymbol{x}\|^2 - 2x_i x_k. \qquad (J.8)$$
<figure><img src="image_10.png" /><figcaption>Figure 7: Projection of samples, $\vec{x}_T$, sketching the support of the invariant measure, $p_\infty$, for a low-rank tensor I.10, $d = 4$ and $\vec{x}_0 = \mathcal{N}((1, 1, 1, 1), 0.01I_4)$. The top plot is in space $(x_1, x_2, x_3)$, the bottom plots are, from left to right, in space $(x_1, x_2)$, $(x_1, x_3)$, $(x_1, x_4)$, and $(x_4, x_3)$.</figcaption></figure>

Hence $\mathbb{E}[Dxx^\top D] = 2(\|x\|^2 I_d - xx^\top)$. Consequently for any $k = 1, \dots, d$ it holds that

$$ \mathbb{E}[(G^k x)(G^k x)^\top] = \frac{1}{4} \mathbb{E}[(M^k - (M^k)^\top)xx^\top((M^k)^\top - M^k)] = \frac{1}{2}(\|x\|^2 I_d - xx^\top). \quad (J.9) $$

As a result, the matrix $G(x)$ has columns $G^k x \stackrel{iid}{\sim} \mathcal{N}(0, V)$ with

$$ V = V(x) = \mathbb{E}\Sigma(x) = \frac{1}{2}(\|x\|^2 I_d - xx^\top), \quad (J.10) $$

of rank $d-1$. Therefore, $\Sigma(x) = G(x)G(x)^\top \sim W_d(V(x), d)$ is a Wishart matrix.

Let $K = K(x)$ be a matrix

$$ K = [K_1, \dots, K_{d-1}], \quad (J.11) $$

with column vectors $K_i$ forming an orthonormal basis of the hyperplane $x^\perp$. Then by construction, we have

$$ \frac{2}{\|x\|^2} K V K^\top = I_{d-1}. \quad (J.12) $$

This means, that

$$ \frac{\sqrt{2}}{\|x\|^2} K G(x) = \left( \frac{\sqrt{2}}{\|x\|^2} K G^1 x \dots \frac{\sqrt{2}}{\|x\|^2} K G^d x \right), \quad \frac{\sqrt{2}}{\|x\|^2} K G^k x \stackrel{iid}{\sim} \mathcal{N}(0, I_{d-1}). \quad (J.13) $$

Therefore,

$$ \frac{\sqrt{2}}{\|x\|^2} K \Sigma(x) \left( \frac{\sqrt{2}}{\|x\|^2} K \right)^\top = \left( \frac{\sqrt{2}}{\|x\|^2} K G(x) \right) \left( \frac{\sqrt{2}}{\|x\|^2} K G(x) \right)^\top \sim W_{d-1}(I_{d-1}, d), \quad (J.14) $$
is a Wishart matrix, in particular $W_p(C, n)$ denotes the Wishart distribution with $n$ degrees of freedom. In the case $n \ge p$, such matrix is invertible almost surely (Muirhead, 2009, Theorem 3.1.4). In our case $n = d > p = d - 1$ thus almost surely

$$ \mathrm{rank} \left( \frac{\sqrt{2}}{\|x\|} \mathbf{K} \mathbf{G}(x) \right) = \mathrm{rank} \left( \left( \frac{\sqrt{2}}{\|x\|} \mathbf{K} \mathbf{G}(x) \right)^{\mathrm{T}} \right) = d - 1. \quad (J.15) $$

Now, since $\mathbf{G}(x)^{\mathrm{T}}x = 0$ we obtain almost surely that

$$ d - 1 = \mathrm{rank} \left( \frac{\sqrt{2}}{\|x\|} \mathbf{K} \mathbf{G}(x) \right) \le \mathrm{rank}(\mathbf{G}(x)) \le d - 1, \quad (J.16) $$

which yields the claim.

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

## J.2 Tensor renormalization

In practice, we renormalize the tensor $\mathbf{G}$ as follows:

$$ \mathbf{G} = \frac{\sqrt{d}}{\|\tilde{\mathbf{G}}\|_2} \tilde{\mathbf{G}} \quad \text{with} \quad \tilde{\mathbf{G}}_{ij}^k = \frac{1}{2}(M_{i,j}^k - M_{j,i}^k). \quad (J.17) $$

The normalization ensures that the trace of matrix defining the Itô term of our forward SDE – i.e. the term driving the exponential decreases of $\mathbb{E}\vec{x}_s$ (see the forward Itô SDE equation D.9) – is

$$ \mathrm{tr} \left( \frac{1}{2} \sum_k \mathbf{G}^k \mathbf{G}^k \right) = - \mathrm{tr} \left( \frac{1}{2} \sum_k \mathbf{G}^k (\mathbf{G}^k)^{\mathrm{T}} \right) = - \frac{1}{2} \sum_k \|\mathbf{G}^k\|_2^2 = - \frac{1}{2} \|\mathbf{G}\|_2^2 = - \frac{1}{2} d, \quad (J.18) $$

similarly to the trace of the matrix defining the Itô term of classical Ornstein Uhlenbeck forward SDE:

$$ \mathrm{tr}(-\mathbf{I}_d) = -d. \quad (J.19) $$

This normalization helps to better control the speed of convergence of the forward SDE without changing its skew-symmetry nor its rank.

## J.3 Mean speed of convergence with renormalized tensor

Note that in this case, for $(\mathbf{x}, \mathbf{y}) \in S = \{(\mathbf{x}, \mathbf{y}) \in \mathbb{S}^{d-1} \times \mathbb{S}^{d-1} | \mathbf{x} \perp \mathbf{y}\}$,

$$ \mathbb{E}\|\mathbf{G}(\mathbf{x})\mathbf{y}\|^2 = \mathbf{y}^{\mathrm{T}}\mathbb{E}\mathbf{\Sigma}(\mathbf{x})\mathbf{y} = \mathbf{y}^{\mathrm{T}}\mathbb{E}\left(\sum_{k=1}^{d}(\mathbf{G}^k\mathbf{x})(\mathbf{G}^k\mathbf{x})^{\mathrm{T}}\right)\mathbf{y} = \mathbf{y}^{\mathrm{T}}\frac{d}{2}(\|\mathbf{x}\|^2\mathbf{I}_d - \mathbf{x}\mathbf{x}^{\mathrm{T}})\mathbf{y} = \frac{d}{2}. \quad (J.20) $$

So, we can expect exponential convergence of the Fokker-Planck equation with the speed

$$ \mathbb{E}[\alpha(\mathbf{G}, d)] = (d-1)\mathbb{E}\min_{(\mathbf{x}, \mathbf{y}) \in S} \|\mathbf{G}^{\mathrm{T}}(\mathbf{x})\mathbf{y}\|^2 = (d-1)\mathbb{E}\|\mathbf{G}^{\mathrm{T}}(\vec{\mathbf{x}}_0)\mathbf{y}_0\|^2 = \frac{1}{2}d(d-1). \quad (J.21) $$

Therefore, the convergence gets faster when the dimension increases.

However, the tensor $\mathbf{G}$ is normalized (see equation J.17), so the evaluation of the convergence speed is modified. We note first that:

$$ \mathbb{E}\|\tilde{\mathbf{G}}\|^2 = d \mathbb{E}\|\tilde{\mathbf{G}}^1\|^2 = \frac{d}{4} \mathbb{E}\|M^1 - (M^1)^{\mathrm{T}}\|^2 = \frac{d}{2} \sum_{ij} (\mathbb{E}(M_{ij}^1)^2 - \mathbb{E}M_{ij}^1 M_{ji}^1), \quad (J.22) $$

$$ = \frac{d}{2} \sum_{ij} (1 - \delta_{ij}^2) = \frac{1}{2} d^2 (d - 1). \quad (J.23) $$
So, we obtain an estimate by Cauchy-Schwartz and Jensen's inequality

$$ \mathbb{E}\|\mathbf{G}(\boldsymbol{x})\boldsymbol{y}\|^2 = \mathbb{E}\left\|\frac{\sqrt{d}}{\|\tilde{\mathbf{G}}\|}\tilde{\mathbf{G}}(\boldsymbol{x})\boldsymbol{y}\right\|^2 = \mathbb{E}\left[\frac{d}{\|\tilde{\mathbf{G}}\|^2}\|\tilde{\mathbf{G}}(\boldsymbol{x})\boldsymbol{y}\|^2\right], \quad (J.24) $$

$$ \le \mathbb{E}\left[\frac{d}{\|\tilde{\mathbf{G}}\|^2}\right] \mathbb{E}\|\tilde{\mathbf{G}}(\boldsymbol{x})\boldsymbol{y}\|^2, \quad (J.25) $$

$$ \le \frac{d}{\mathbb{E}\|\tilde{\mathbf{G}}\|^2} \mathbb{E}\|\tilde{\mathbf{G}}(\boldsymbol{x})\boldsymbol{y}\|^2, \quad (J.26) $$

$$ = d \frac{d/2}{d^2(d-1)/2}, \quad (J.27) $$

$$ = \frac{1}{d-1}, \quad (J.28) $$

and finally we obtain the following bound

$$ \mathbb{E}\alpha(\mathbf{G}, d) \le (d-1)\mathbb{E}\|\mathbf{G}^\top(\tilde{\boldsymbol{x}}_0)\boldsymbol{y}_0\|^2 = 1. \quad (J.29) $$

# K Going beyond the rank condition for MSGM scalability

The dense tensor of Appendix J imposes a computational complexity as $O(d^3)$. To scale up the method, we shall consider sparse tensor **G**. However, the rank condition A2 makes it difficult to find sparse tensors. Therefore, we here open the discussions to a weaker set of assumptions.

## K.1 Weaker assumptions

We recall here the two main assumptions of the paper

Skew-symmetry: For any $k \in \{1, \dots, d\}$, the matrix $\mathbf{G}^k = (\mathbf{G}_{i,j}^k)_{i,j}$ is skew-symmetric. (A1)

Rank condition: For any $\boldsymbol{x} \in \mathbb{R}^d \setminus \{0\}$, $\mathrm{rank}(\mathbf{G}(\boldsymbol{x})) = d-1$. (A2)

Note that the Fokker-Planck equation 3.4, Proposition F.2, Proposition F.1, and Theorem 3.4.1 require only the assumption A1. So, the backward SDE, ODE and score-matching loss are general enough and do not prevent the use of sparse tensor **G**. In contrast, our current proof of the asymptotic results Theorem 3.1.1, Theorem D.4.1, and Theorem 3.3.1 rely on the restrictive assumption A2, and unfortunately, it seems difficult to find a sparse tensor **G** matching this assumption.

### K.1.1 Rank condition almost everywhere

Therefore, we discuss here a weaker set of assumptions where the noise rank condition A2 is verified for almost all $\boldsymbol{x} \in \mathbb{R}^d$ only. This set of assumptions will yield a definition of a sparse tensor in Appendix K.2.2 providing satisfactory numerical results in practice.

Skew-symmetry: For any $k \in \{1, \dots, d\}$, the matrix $\mathbf{G}^k = (\mathbf{G}_{i,j}^k)_{i,j}$ is skew-symmetric. (A1)

Rank condition almost everywhere: For almost all $\boldsymbol{x} \in \mathbb{R}^d$, $\mathrm{rank}(\mathbf{G}(\boldsymbol{x})) = d-1$. (A3)

The assumption A3 means the set $A_G = \{\boldsymbol{x} \in \mathbb{R}^d | \mathbf{G}(\boldsymbol{x}) < d-1\}$ has zero Lebesgue measure, i.e. $\int_{A_G} dx = 0$. Obviously, the assumption A2 implies the assumption A3.

The right-hand side of the Fokker-Planck equation 3.4 is a function of $\nabla_{\perp} p_{\infty}$ only. Hence, under the weaker assumptions A1 and A3, rotational invariant distributions are still invariant measures of the Fokker-Planck equation. Following the proof of Theorem D.2.1, we saw that the invariant density is characterized by $\|\mathbf{G}(\boldsymbol{x})^{\top} \nabla_{\perp} p_{\infty}(\boldsymbol{x})\| = 0$ almost surely. So, if $\mathbf{G}(\boldsymbol{x})$ has rank $d-1$ for almost all $\boldsymbol{x}$ with respect to the Lebesgue measure, then this requires $\nabla_{\perp} p_{\infty}(\boldsymbol{x}) = 0$ almost everywhere. Therefore, the invariant measures of Fokker-Planck equation 3.4 must be rotational invariant almost everywhere.

However, the existence and uniqueness of a classical solution (Lemma D.4.1) and the convergence guarantees to the invariant distribution (Theorem D.4.1) need more careful analysis. We only sketch some challenges involved, the full analysis will be carried out in a follow up work.
If the diffusion process enters the area of points $x$, such that $\text{rank}(G(x)) < d-1$, one has to make sure that the Diffusion process is not trapped in such an area, even if it has a measure zero. In particular, let $D \subset S^{d-1}$ be defined as

$$D = \{ \boldsymbol{x}^n \in S^{d-1} \mid \text{rank}(G(\boldsymbol{x}^n)) < d-1 \}.$$

Then, we call $D$ a trap set. If the process $\vec{\boldsymbol{x}}_s^n$ once ever entering $D$ with positive probability, it cannot leave $D$ again, i.e.

$$\mathbb{P}(\vec{\boldsymbol{x}}_s^n \in D, \forall s \geq s_0 | \vec{\boldsymbol{x}}_{s_0}^n \in D) > 0.$$

Hence, in this case, convergence to the correct invariant measure has to ensure that the trap set is not invariant under the diffusion-controlled dynamic. Such an analysis then sufficiently can be implied by Hörmander / bracket-generating conditions, i.e. hypoellipticity analysis. Based on this, the asymptotic results Theorem 3.1.1, Theorem D.4.1, and Theorem 3.3.1 must be adapted. In this case we expect the convergence rate to the invariant distribution to be slower compared to the exponential convergence rate obtained in the case of strong rank condition, see also Appendix K.3 for a related discussion.

Although a detailed analysis of this research question is out of the scope of the current manuscript, we want to stress its relevance related to the *scalability* of the proposed method for the high-dimensional case. The standard construction via a random dense tensor $G$ poses scalability problems. On the other side, sparse tensors provide a tool to enable such scalability provided that they satisfy the (weaker) rank conditions. While the non-local sparse tensor discussed in Appendix K.2.3, satisfy the strong rank condition, the local sparse tensors from Appendix K.2.2 only satisfy the weak rank conditions. Still, the latter have been applied in our numerical investigation for the high-dimensional test cases with Particle Image Velocimetry measurements as discussed in Appendix M.7 yielding first very promising results.

### K.1.2 Itô term rank condition

Now, we discuss another weaker set of assumptions where the noise rank condition A2 is replaced by an Itô drift rank condition. Although attractive, a detailed analysis in Appendix K.2.1 will lead us to consider this set of assumptions as insufficient for the MSGM sampling procedure.

Skew-symmetry: For any $k \in \{1, \dots, d\}$, the matrix $G^k = (G_{i,j}^k)_{i,j}$ is skew-symmetric. (A1)

Itô term rank condition: the matrix $S := \frac{1}{2} \sum_{k=1}^{d} (G^k)^T (G^k)$ is full rank (A4)

From Lemma D.1.1, we note that $S = -\nabla\nabla \cdot \Sigma$. Lemma D.1.2 gives the Itô forward diffusion which can be expressed with $S$. The assumption A2 is not needed for these lemmas. These results are true as long as the assumption A1 is verified. Taking the expectation of the Itô diffusion, we get:

$$\frac{d}{ds} \mathbb{E} \vec{\boldsymbol{x}}_s = -S \mathbb{E} \vec{\boldsymbol{x}}_s. \qquad (K.1)$$

Instead of controlling the convergence of the full distribution $p_s$, the assumption A4 controls the convergence of the mean only. It leads to the following property justifying our assumption choice.

**Proposition K.1.1.** Let the assumption A1 holds. Then, the following assertions are equivalent

* The assumption A4 holds.

* $\mathbb{E}[\vec{\boldsymbol{x}}_s | \vec{\boldsymbol{x}}_0] \underset{s \to \infty}{\longrightarrow} 0$.

* $\text{Var}[\vec{\boldsymbol{x}}_s | \vec{\boldsymbol{x}}_0] \underset{s \to \infty}{\longrightarrow} \|\vec{\boldsymbol{x}}_0\|^2$.

*Proof.* $S$ is positive semi-definite, so it is diagonalizable in an orthonormal basis, and from equation K.1, $\mathbb{E}[\vec{\boldsymbol{x}}_s] \underset{s \to \infty}{\longrightarrow} 0$ if and only if $S$ is positive definite, i.e. the assumption A4 is verified.

Besides, by assumption A1, the norm $\|\vec{\boldsymbol{x}}_s\|$ is conserved along the diffusion, so

$$\text{Var}[\vec{\boldsymbol{x}}_s | \vec{\boldsymbol{x}}_0] = \|\vec{\boldsymbol{x}}_s\|^2 - \|\mathbb{E}[\vec{\boldsymbol{x}}_s | \vec{\boldsymbol{x}}_0]\|^2 = \|\vec{\boldsymbol{x}}_0\|^2 - \|\mathbb{E}[\vec{\boldsymbol{x}}_s | \vec{\boldsymbol{x}}_0]\|^2 \quad (K.2)$$

which converges to $\|\vec{\boldsymbol{x}}_0\|^2$ if and only if $\mathbb{E}[\vec{\boldsymbol{x}}_s | \vec{\boldsymbol{x}}_0] \underset{s \to \infty}{\longrightarrow} 0$. $\square$