Multiplicative Diffusion Models: Beyond Gaussian Latents
Authors: Robert Gruhlke
Valentin Resseguier
Merveille Talla
Venue: ICLR 2026 Poster (2026)
Keywords: score-based diffusion
generative modeling
multiplicative noise
non-Gaussian latent variables
conservative dynamics
heavy-tailed distributions
Fokker–Planck equation
Primary area: generative models
TL;DR: Generative Modeling with Multiplicative Score-Based Diffusions

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 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 data (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.

Published as a conference paper at ICLR 2026
MULTIPLICATIVE DIFFUSION MODELS: BEYOND GAUSSIAN LATENTS
Robert Gruhlke', Valentin Resseguier", Merveille Merveill Freie Universit at Berlin, Berlin, Germany
T
PINRAE, OPAALE, Rennes, France
Inria centre at Rennes University, Rennes, France
r. n.gruulkefuu -er.lliiee--bbrii......
(valentin, resseguier, merveille wtt.:llleeaaaaaaiiteee ves
ABSTRACT
We introduce new class of generative models based on multiplicative score-driven diffusion. In contrast classical diffusion models that rely on additive Gaussian noise, our construction is driven by skew-symmetric multiplicative noise. It yields conservative borward-backward dynamics inspired by the principles of physies. We prove that the forward process converges exponentially fast 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 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 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 principled alternative to current bosee generative models, with strong potential for domains where rare but critical events dominate.
INTRODUCTION
Mathematically equivalent (Song et al., 2021), diffusion models and score-based generative models demonstrate impressive skills and among the current state-of-the off-the-art for the generation of two- and three-dimensional images. Unconditioned sampling scores can be easily modified to conditioned sampling scores address various inverse problems (Rybchuk 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., [0044] Stamatclopoulos & Sapsis 2025). Those generative AI challenges may 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 al. 1991 Klyatskin 1994. Piterbarg & Ostrovskil
1997
Mikulevicius & Rozovskii 2004| Mémin 2014 Holm 2015: Resseguier et al. [022 [hen
et al. 2023) and, more generally. from slow-fat flot-fast&o systems with multiplicative noise (Majda al.
1999 Franzke et al., 2005 Gottwald & Melbourne, [013] Gottwald & Harlim] 2013). Transport noise
models may be understood generative models based stochastic fluid dynamics rather than fitted

---

Published as a conference paper at ICLR 2026
neural networks. As other generative models, they suit particularly well to Bayesian inverse problems (Cotter et al. 2020b a. Resseguier et al. 2022] Dufée et al. (022))
Here, 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 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 process is multiplicative. We call it Multiplicative Score-based Generative Model (MSGM). Involving random rotations around the origin, greatly differs from previous diffusion models and opens new research path. The key aspects are summarized in Figure
type of diffusion model where the noising
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-Plaan Then, we separately analyze the norm and direction of the diffusion. The norm is steady, whereas the direction follows 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.
Deep theoretical analysis of MSGM:
equation of forward diffusion and its invariant measures.
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 -an-Gassiian 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 numerical procedure mimic the heavy-tail distribution with MSGM. We add 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.
forward SDE:
dæ
G( dB,
MSGM (multiplicative noise)
larget spicce
d'
G le (ao ( fe,, T dt odB:)
latent spacee
L.)
SGM (additive noise)
dæ dt + V2 ( T t)at cdB.)
forward SDE: dæ
d. v2 oBB
Figure 1: Illustration of multiplicative score-basee boore generative modeling (ours) compared to additive score-based generative modeling.
2

---

Published as a conference paper at ICLR 2026
2 ADDITIVE SCORE-BASEL GENERATIVE MODEL
2.1 FORWARD AND BACKWARD SDES
Diffusion models 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ē. 2 .s + vddB. 81
(2.1)
where 2, Rd is distributed according to some density ps for 0, s + B, is d dimensionnl Brownian motion, and Tg distributed according to the dataset interest. is an process: the continuous-time version of first-order autoregressive (AR) model and the distribution P. convergences to a standard Gaussian density exponentially for 8 + oc, e.g. in total variation or Wasserstein distance. We can then define for E [o, I] the backward equation
dæ = x + 2V ( )t + Vd.
(2.2)
with + B, another d- dimeennsional Brownian motion and æo pr (identifying the density pr with its distribution). Then for any s E [0, T). æт-2 and w have the same law ps. In practice, when an approximate score V log Pr-t is available we initialize equation 2.2 with standard Gaussian distribution 2o ~ NO, Ia) and integrate the backward SDE from 1 to / = T (i.e. from 8 = T to 0), ideally letting zy become statistically similar to our dataset of interest.
2.2 A NEURAL NETWORK TO FIT THE SCORE
e.g., In practice, the score Pr is approximated by surrogate model, (T e.g. fitted artificial neural network (ANN). Alternatively, one can work on ag (xx T - (2,0 T t) (Huang et al. 2021). For large-dimensional problems, [Snn et al. (2020) proposes to learn this neural network by Sliced Score Matching (SSM). Here, ag is obtained by minimizing the loss function
CSGM SSM (0)
1-"
-T
Ey Rad((d) E. s.
| |1l0o +)] (v. (( a., 8) - ).)) v)] ds. (2.3)
where ||-|| is the Euclidean norm, Rad(d) denotes the d-dimensional Rademacher distribution and E= is the expectation along each path realization x., Section A details the most common score matching losses and their link to the concept of the Evidence Lower Bound (ELBO).
MULTIPLICATIVE SCORE-BASED GENERATIVE MODEL.
Rather than relying on additive SDE equation 2.1 propose 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 I, and the backward process E., respectively.
3.1 FORWARD SDE
Instead of considering a forward SDE with additive noise, we rely on multiplicatiie noise model
dē G( dB...
(3.1)
where d >>3 2, G Rdd Rdxd is linear and o stands for the Stratonovich notation. The readers unfamiliar with this notation may interpret the Stratonovich noise + Bd 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:
} 0 ss as) G ((*z)) ( B. B.-da)-
(3.2)
3

---

Published as a conference paper at ICLR 2026
For deeper understanding, Section BasiiinnnBrre recalls some important notions of stochastic calculus, including the Stratonovich notation and the relationship to Itō calculus. Let be represented by a third-order
tensor [G E Rd.dd and define the random matrix Z, = Σi= G² (B, Then, equation 3.1 can be written more explicitly as:
dæ. = (.
Σ
(odB.)& odZ,w, 2 2 (Z Z Zo-do) (3.3)
k=1
k=1
where the time increments the random matrix, odZ, ~ 2 (Z,+ds ds-d are uncorrelated.
Additionally, will impose two assumptions valid throughout the paper:
Skew-symmetry : For any l E (1....,d) d), the matrix G = (G )&s skew-symmetric. (A1) Rank condition: For any e RR rank(G((x)) = 1.
(A2)
In particular equation 3.1 does describe a geometric Brownian motion, as the noise Z, is not diagonal. It includes zeros on the diagonal due to the skew-symmetry of G A geometric
Brownian motion would necessitates Z, = diag B A strategy to obtain tensor G that matches assumptions A1 A2 will be discussed in Section By linearity, the skew-symmetry of all G (assumption AT implies the symmetyy of the whole multipicaative noise matrix odZ,. This structure inspired by transport noises in fluid dynamies (Kraichnan 1968 Piterbarg Ostrovskii
1997, Resseguier et al. 2021). In this analogy, V. 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 a. Here, we might address problems outside the scope of fluid dynamics, though maintaining the noise skew-symmetry (assumption AT) and thus the energy conservation, as discussed in Section 3.2
With assumption A2 the noise spreads in a large space: Im G( = x. ensures sufficient variability in the noising process and, in turn, a tractable distribution for æT when T becomes large. Theorem 3.1.1. Let the assumptions A7 and A2 hold. Then, the Fokker-Planck equation of equa- tion 3. reads
① = V (x) V [p(()))
po(x)
x E Rd
(3.4)
with conditional noise covariance Exx) G(z)G(a) and v_ denoting the orthogonal projection of nabla V on the tangent plane x i.e.
V1 == (Ia-2"(2")1))
(3.5)
for x" == x/|[>|| with E D = Rd {0} and 0 otherwise, the unit vector orthogonal to the d-sphere. Moreover, any stationary density Poo of equation 3.4 inatiinn-nvariint on IRd.
The proof is detailed in Section D.2 In order to highlight the connection to diffusion models on Riemannian manifolds, we note that V is the Riemannian gradient on the unit d-sphere: gd-1 {i Id ||>||| 1}, Section H. For possible extension of the considered dif-
fusion equation 3.1 to the case of non-zero drifi, we refer to section D.6
3.2
DYNAMICS OF NORM AND DIRECTION
We now consider for any s 2 0 and w. # O, the spherical decomposition
n 0 with F == / ~ sd-1
=
(3.6)
First, we note that the norm, skew-symmetry of odZ, implies that dž = od H. is orthogonal to v. and hence:
|| remains constant throughout the noising process. Indeed, the
d 11 2σ, odæ = O, Vs 2 0.
(3.7)
Consequently, | = || The vector H. moves randomly on | the d-sphere of radius o Therefore, the distribution of the norms of the latent variable is exactly the distribution of
4

---

Published as a conference paper at ICLR 2026
the norms of the points of the dataset. We refer to Section D33 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 partieular, 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, will fit the distribution of the log-norm, Flog l.l., with eCDF. We refer to Section C for details and an overview of sampling from one-dimensional distributions.
- We now focus on p" the distribution of the direction, E., in particular as 00. For better readability, we postpone the full discussion and the main theorems Section D.4 Lemma D introduces the Fokker-Planck equation of the direction on gd-1 and its unique invariant measure, ph Then, Theorem D.4.2] shows the exponential convergence of the initial distribution pi pi= the uniform distribution on the unit sphere s Consequently, since g is compact, this implies convergence in total variation of p p'c and convergence in distribution of æ to x (($-11)..
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 x² distributions with degrees of freedom (see Section E.2) This property differs from the usual SGM. SGM latent variables are Gaussian, leading to X distributions for the norms of the latent variables regardless of the data set. According to Lafon 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 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 Section E.6 More generally, Section 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 Section E.I) Moreover, the norm and direction are independent from each other, which will drastically facilitate the sampling procedure. These results hold for latent vectors xo ~ Pao. In practice, integrations of forward and backward diffusions only possible over a finite time T. However, the following theorem states that the law of the solution, TT, will become close to Pso exponentially as fast as T -> +00. So, we can confidently rely on finite-time integration.
Theorem 3.3.1. Let assumptions Al and A2 hold. Let 0o po E CD and let P\J be the (radial)
density of @oll. Then, the Pokker-Planck equation 3.4 has a unique solution pa E (D)
for all s > 0. Moreover, the Fokker-Plance equation has the stationary distribution
P1.1 (1|==l) 1
Pro (x)
(3.8)
Th In particular. || and x. are asymptotically independent s - +00. Moreover, there exists (G(G, d) > 0 such that
Ip. Poo (B I cxp(-as))(po Poo (12=
(3.9)
The proof and details on a are given in Section D.5 and a specific case is discussed in Section J.3
[gd-1| is the volume of the scaled d-sphere The factor ||a|||-d in equation 3.8 is expected. Indeed, T T-
|| i.e. it corresponds to the uniform distribution on the scaled d-sphere ||
We will now consider the practical question on how draw samples from the latent distribution with density Pso from equation 3.8 It is of product structure between the radial and the directional
component. So, we sample the norm Ros and the direction Fx separately and multiply them. The
5

---

Published as a conference paper at ICLR 2026
norm Reo can be sampled from one-dimensional approximation of the data norm (see Section and the direction oe is uniform. So, sample xc. ~N(0,, Ia) and set
-.V
2 =
+n R 20 with
on
2
-N 200 xl
+N
Rs = f(1
-V
(3.10)
and f(r) = exp log (F (F2((a(12)))
-c,
44 >0.0
(3.11)
-. If x = 0, we set x = 0. Proposition E.3.2 shows that this procedure leads to samples with the correct distribution, up to the approximation of the log-norn CDF Fing I-le ~ 1-k. Moreover, the direct map Flng l-le can transform latent vector, Iy, into a Gaussian one (see Section E.4 This transformation may be useful for future applications like inverse problems time evolution fittings.
3.4 REVERSE ODE/SDE AND SCORE MATCHING
From the Itô forward SDE (see Lemma D.I.2). we know that the Stratonovich reverse SDE writes
dz. = G( loggr- + odB:)
(3.12)
and the reverse probability flow ODE is given
dæ
dt ((..) (G (G(Z:))VV log pr prrppp--r))) ((.))
(3.13)
The corresponding derivations are formulated in Proposition the appendix using Anderson 1982; ; Song et al. (202I). Following Huang et al. [2021) we directly model G ) log Pr x;;) by neural network as T t). Additionally,
and Proposition F.2 and are proven in
incorporate spherical input layer, see Section L.4.1 We fit the parameters 0 by sliced score matching (SSM) (Song et al., 2020), because in the multiplicative case we do not have analytic formula for the conditional score 1, ()) and because of the better scalability of SSM to high-dimensional problems that would like to address in the future. To this end, minimize the loss function:
CSSM () 10 E E (od-)) ( +|12 (v G(Z.)ao s)) v (3.14)
where (ad(d) denotes the Rademacher distribution in Rd The following theorem states that even in our multiplicative case, score matching is equivalent to maximize the ELBOO E In line with Benton al. (024)); [Ren al. (2025), this theorem generalizes the result of Huang 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 Section
Theorem 3.4.1. Let assumption
holds. Then, there exists a constant C such that
Po(x) > Eoo (x) = C
E- m
(ll V (G(z))ao( 8)) 2o ds. (3.15)
We proof this theorem in Section G The first term C E (oggr(@T)) = x is a constant w.r.t.
0. So, 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 itself equivalent to explicit score matching (ESM), denoising score matching (DSM), and SSM (Huang et al. 2021). ). Note that formally replacing by V2, 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.
WORKFLOW
Algorithm summarizes the proposed MSGM procedure. Here we make use of color to highlight the differences compared to SGM. For more details, refer to Section

---

Published as a conference paper at ICLR 2026
Algorithm 1: MSGM (Multiplicative Basee Generative Model).
1 }M Input: tensor G, one-dimensional distribution model fy. data { f Im=1- te. time horizon T, time steps Nz. and Ný. time scheduler g. score model dg. initial ANN parameter fo. iterations Niter
Training stage-
1: 7" < fit
(f,, <1 mm 1I) M
(Fitting of Jax see Section
2: for n = 0 to Niter - do
3:
2o ~ 1
M Em=1 8_mm
(Random mini-batch from dataset)
20 U(T(te: T))
4:
(Sample uniform gridded time
5:
I. + SRKA (s,, L N Tu, 20,0, gG)
(Forward diffusion integration via Algorithm |Slicing directions
6:
V~ Rad(d)
7:
((0m) + SSM (s, a; gao.,, v)
(Score-matching loss, from equation 3.14 {e.g. via ADAM)
8:
0nt1- opttmiizee update ((O.))
9:
end for
10: 0" ONar
(Set final ANN parameter]
Generative sampling stage-
-N
11: To ~ *N(0,, Ia);
(Sample strong white noise"
12: Zo = f.- (1
(Sample weak white noise, see equation 3.11 (Reverse diffusion integration via Algorithm
3 I
13: TT <- SRK4 N oo gGgg,
+
14: return 7", 0", aT
RELATED WORKS
Combining machine learning and mechanistic approaches now common approach. We may cite
physics-informie informed neural networks (PINNs) (Raissi et al. [2019] Lu & Xu (0224,, 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., [024] Kim et al. (022). 2). and chaos from energy-based models
(Fournier & Pierfrancesco 2025) among others. Here, we shall focus score-based generative
models. Bastek et al. (2024) add the physical equations inside their score matching loss. Holzschuh
ct al. (2023) fit score correct a backward physical equation but does not propose any generative
model. To denoise corrupted images, several authors (e.g. [hou et al. 2014 [haa 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 al. (0255) do, but their
framework simplifies to SGM by considering the pixel-wise logarithm of images. Guth et al. (2022): ;
Lempereur & Mallat (2044) encode a target multiscale structure (c.g. turbulence) by a hierarchy of
normalized wavelets conditioned by the larger scales. [Cen & Vanden-Eijjnen (2025) adapt the noise to that multiscale structure in a stochastic interpolant context. Lobbe et al. 2023]2222 replaced the
Gaussian process involved in the transport-noise equations (Kraichnan) 1968) Piterbarg & Ostrovskii
1997 Resseguier et al., 202T) by Shrodinger bridge (De Bortoli et al. (021). They inserted a SGM
inside transport noise dynamics, whereas inserted dynamics similar to transport noise inside
SGM. Following general Bayesian approaches, some of the litcrature on transport- relies on
the Girsanov theorem fit drift modification evaluate likelihood (Cotter al. 2020a Singh
2021 to SDEs inspired by et al. 2025) Section G.9. By extending the ELBO of (Huang et al.]
transport noise, Theorem 3.4. justifies our fit of the backward SDE drift.
Several authors have recently proposed Langevin equations (Amaudon et al.
2019
Luesink & Street
2025 Ayala al. 2025) and SGM (De Bortoli et al., 2022; Huang al.
2022 Lou et al. 2023
Benton et al. 2024)) on Riemannian manifolds in order to generate data lying on particular manifold.
Clearly different, our goal is more classical: generating data in Rd 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 ||2_||//[1 - De Bortoli
cet al. (2022) describes diffusions and SGM in the d-sphere gd-1 A detailed comparison is given in
Section
7

---

Published as a conference paper at ICLR 2026
Regarding extremes generation, using variational autoencoders (VAE), Lafon et al. (2223 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évv and Student-f) latent distributions. Our approach automatically makes the tails of the latent distribution fat when necessary. It learns from the distribution of the data norm, Pl.- Similarly, the diffusion proposed by Dharmakeerthi et al. (2055) adapt to data but through a nonlinear drift and an additive noise. [i et al. (2024); Price et al. (225);; 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 suggest that our method probably overcomes this limitation. To represent the directionality of extremes, many authors decompose norms and directions of extreme events (Engelke al., 2019; Palacios-Roolee 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. 22019; Lafon et al.) (023). Build on random rotations, MSGM naturally suggests such 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 Section H.2)
EXPERIMENTS
For numerical experiments, we choose define tensor G in a simple way. We sample d random matrices, keep only their skew-symmetric parts, and normalize:
dd G with G.y = (MM, - MY.) and M snd N/ 1).
G
(6.1)
In Section J we show that this random tensor G respects conditions Al and A2 almost surely, Section K proposes alternative tensor definitions with sparse structures that allow high-dimensional applications. Following Section K.2.2. MSGM can generate images as in Section M.7.2 For the test cases below, we also checked in Section M.6.2Jand Section 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, Tca, with scale parameter ry:
7/T odd ~ PCa with PCa(x) = x² +2
(#ca):
(6.2)
It is an extreme of fat-tailed distributions with a power-law tail: (Ca(x) |x|-2 for large |z]. Real problems often involve both correlation and dimensionality d > 2. So, we correlate Cauchy variables, 2n Axca. with fixed matrix, A, initialized with i.i.d. coefficients A,,, 1).
Figure@ 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 Section E.5 and Section 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 of--oh--thhh the--rt art SGM can be found in Section M.6.
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 (Georgeault & Heitz, 2026). Vorticity quantifies the local rotation speed of fluid and is known have point-wise distributions with tails fatter than Gaussian ones (Wilezek & Frieddich 2209). We focus on a benchmark fluid flow: a wake flow at Reynolds number 3900 created by circular cylinder embedded in a mean stream (Parnaudeau et al. 2008). Each vorticity sample
8

---

Published as a conference paper at ICLR 2026
test gen.
test gen.
10
10
10
10
10
10
10
10
10
10
-10
10
-10
-10
10
-10
10
-10
10
-10
-10
10
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 vector of 4 correlated Cauchy variables. On the diagonal, log- oig-histogram and logarithm of the pdf KDE estimation are superimposed.
gen.
gen.
-3
3
-3
3
Figure 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 the pdf KDE estimation are superimposed.
is evaluated at d = 16 spatial points to ensure low dimensionality. We use limited training data (1024 data points) to make rare events even more rare and learning more challenging. Section M.7 provides a deeper description of this experimental dataset. Figure 3 highhlghht 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, theoretically suggested by Section E.5 and experimentally verified in Section M.7 In particular, the MSGM latent distribution seems have Laplace tails and to be more accurate the low-data regime (see Figure 4). Additional we carried out high-dimensional experiments with d = 1024 in Section MM..2] based on sparse tensor G developed in Section K More details data, preprocessing. illustrations, and other numerical experiments are given in Section M
CONCLUSION AND DISCUSSION
We have proposed a new type of diffusion model with multiplicative noise. After a theoretical analysis of this ansatz, algorithm is specified to mimic fat-tailed distributions, surpassing SGM in this task.

---

Published as a conference paper at ICLR 2026
109
10
MMD
10-1
MSGM
train data
10
10-2
28
212 2.6
2
224
26
2
2
effective number of iterations
Training set
(a) Single MMD value for the correlated Cauchy distribution (d = 4) as function of the number of effective ADAMS iterations.
(b) Mean and 80% C1 of MMD for the vorticity measurements
distribution (d 16) as function of the number of training samples.
Figure 4: Convergence behaviors of the Maximum Mean Discrepancy (MMD) for two different test. cases. 10d 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 be
integrated numerically since we do not know an analytic solution for large-rank tensors (see Section [ for the solution with low-rank tensors). It implies either 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, 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 diffusions d-spheres. In addition, random matrix theory and free probabilities (Biane 1997: Delyon Yao, 2006; Demni 2008 Lévy] 2008 Delyon, 2010 Demni & Hmidi [012] 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 unitary Brownian matrix, converging for large dimensions to a free multiplicative Brownian matrix. Both theories could facilitate the sampling and the score evaluation the MSGM forward diffusion. Moreover, dense third-order tensor G prevents image processing and other rige-dimensional applications, related to, say, turbulent fluid dynamics. In fact, dimensions d of such problems are very large sasee-typon typically d O(10%) or more. A dense tensor G as we use in our numerical experiments has d coefficients, and the memory and computational costs would become prohibitive these cases. To address this issue, Section K proposes several sparse tensors and alternative to assumptions And images in dimension d 1024. Furthermore, we are currently developing physics- bessec sparse tensors G. Here, MSGM forward SDE is the spatial discretization of stochastic partial differential equation involving transport noises (Kraichnan 1968 [Piterbarg & Ostrovskil 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.
Section M.7.2 shows first MSGM generated
In addition to the improvements discussed above, many paths remain to be explored. First, our theoretical results could be generalized other multiplicative diffusions. We have considered dense linear maps Φ + G(x) with Im G(=) x- for any æ * 0. We believe that sparse linear maps of Section K and non-linear Lipchitz-continuous maps can yield similar theoretical results as long as that image condition is fulfilled for almost every x E 12d (see Section K.1). The non-linear case would include particular sphere-wise diffusions of De Bortoli al. (2022) (see Section H.2). Second, we could address dynamical system forecasting. With the Gaussianization of MSGM latent vectors (see Section 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 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 al. 2025).
10

---

Published as a conference paper at ICLR 2026
ACKNOWLEDGMENTS
We thank Nizar Demni for proofreading, Erwin Luesink, Gilles Carron, and Erwan Brugallé for comparison with Brownian motions on hypersheres, Thomas Optiz, Philippe Naveau, Etienne Mémin, Dominique Heitz, and Théo Voillemin for discussions, Jocelyn De-Goer-De-Here and colab. for sharing GPUs. This work has been funded by the joint INRAE-Inria PhD grant, the ERC project $56408-STUON and by DFG under Germany's i Excellence Strategy (EXC-2046/2, project ID 390685689).
REFERENCES
Ronald 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):3133326, 1982.
Hassan Arbabi and Themistoklis Sapsis. Generative stochastic modeling of strongly nonlinear flows with non- non-gaussian statistics. SIAM/ASA Journal on Uncertainty Quantification, 10(2):555-583, 2022.
Alexis Arnaudon, Alessandro Barp, and So Takao. Irreversible langevin meme 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 Amold. 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 Whitaker. Nonlinear ensemble filtering with diffusion models: Application to the surface quasigeostrophic dynamics. Monthly Weather Review, 153(7):1155-1166) 155-1169, 2025.
Jan-Hendrik Bastek, WaiChing Sun, and Dennis M Kochmann. Physics-informed diffusion models. arXiv preprint arXiv:2403. /4404, 2024.
Maxime Beauchamp, Ronan Fablet, Simon Benaichouche, Piere 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 Arnad Doucet. From denoising diffusions 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:1 12588, 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. Attaari, 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.
11

---

Published as a conference paper at ICLR 2026
Yanlai Chen, Sigal Gottlieb, Lijie Ji, and Yvon Maday. An eim-degradation free reduced basis method via over collocation and residual hyper ceduction-based error estimation. Journal of Computational Physics, 444:1 444:110545, 2021.
Yifan Chen and Eric Vanden-Eijnden. Scale-adaptive generative flows for multiscale scientific data. arXiv preprini arXiv:2509.02977 2025.
Lucia Clarotto, Denis Allard, Thomas Romary, and Nicolas Desassis. The spde approach for spatio-tempore datasets with advection and diffusion. Spatial Statistics, pp. 100847, 2024.
Colin Cotter, Dan Crisan, Darryl Holm, Wei Pan, and Igor Shevchenko. Data assimilation for quasi- geostrophic model with circulation Bieulltion-preserving stochastic transport noise. Journal of Staristical Physics, 179(5):1 186-1221, 2020a.
Colin Cotter, Dan Crisan, Darryl D Holm, Wei Pan, and Igor Shevchenko. A particle filter for stochastic advection by lie transport: case study for the damped and forced incompressible two- dimensional euler equation. SIAM/ASA Journal on Uncertainty Quantification, 8(4): 1444:144444492 2020b.
Colin Cotter, Dan Crisan, and Maneesh Kumar Singh. Data assimilation for the stochastic camassa- holm equation using particle filtering: 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.0008; 2024.
Valentin De Bortoli, James Thornoon, Jeremy Heng, and Arnaud Doucet. Diffusion schrödinger bridge with applications to score-based generative modeling. Advances in Neural Information Processing Systems, 34: 14:11795-17709, 2021.
Valentin De Bortoli, Emile Mathieu, Michael Hutchinson, James Thornton, Yee Whye Teh, and Arnad Doucet. Riemannian score- score-based generative modelling. Advances in Neural Information Processing Systems, 35:2406 -4222, 2022.
Bernard Delyon. Concentration inequalities for the spectral measure of random matrices. Electronic Communications 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 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 ensemble kalman filters. Quarterly Journal of the Royal Meteorological Society, 148 1099 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 nuper-resolution for two-dimensional spatiotemporal turbulence. Computer Methods in Applied Mechanics and Engineering, 433:1 17478, 2025.
12

---

Published as a conference paper at ICLR 2026
Mare Anton Finzi, Andres Potapezynski, Matthew Choptuik, and Andrew Gordon Wilson. A stable and scalable method for solving initial value pdes with neural networks. In The Eleventh International Conference Learning Representations. 2023.
Gerald B. Folland. Real Analysis: Modern Techniques and Their Applications. John Wiley & Sons, New York, 2 edition, 1999. ISBN 978-0-471-31716-6
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-Fijnden. 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 rype. Prentice- Hall, Englewood Cliffs, NJJ, 1964.
Mevlana C Gemici, Danilo Rezende, and Shakir Mohamed. Normalizing flows on riemannian manifolds. arXiv preprint arXiv:1611.02304 2016.
Philippe Georgeault and Dominique Heitz. Non-time-resolved PIV dataset of flow over a circular cylinder at Reynolds number 3900, 2026.
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 lan Melbourne. Homogenization for deterministic maps and multiplicative noise. Proceedings 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. 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.10S5 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, Sciences,471 471 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 manifolds. Number 38. American Mathematical Soc., 2002.
Chin-War Wei Huang, Jae Hyun Lim, and Aaron C Courville. A variational perspective siffusion-based generative models and score matching. Advances in Neural Information Processing Systems, 34: 22863-22876, 2021.
Chin-We Huang, Milad Aghajohari, Joey Bose, Prakash Panangaden, and Aaron C Courville. Riemannian diffusion models. Advances in Neural Information Processing Systems, 35:2750-2761, 2022.
13

---

Published as a conference paper at ICLR 2026
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 Machine Leamiing pp- 4523 4532..
PMLR, 2021.
Aapo Hyvärinen and Peter Dayan, Estimation of non non-normalized statistical models by score matching. Journal of Machine Learning Research, 6(4), 2005.
Priyank Jaini, Ivan Kobyzev, Yaoliang Yu, and Mareus 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:1 10841, 2022.
Peter E Kloeden, Eckhard Platen, Peter Kloeden, and Eckhard Platen. Stochastic differential equations. Springer, 1992.
Valerii Klyatskin. Statistical description of the diffusion of passive tracer in 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 wrXXv:2300.10987, 2023.
John M Lee. Introduction to Riemannian manifolds, volume 2. Springer, 2018.
Kookjin Lee and Kevin Carlberg. Model reduction of dynamical systems on nonlinear manifolds using deep convolutional autocncoders. Journal of Computational Physics, 404:108973, 2020.
Etienne Lempereur and Stéphane Mallat. Hierarchic flows estimate and sample high-dimensional probabilities. arXiv preprint arXiv:2405.0334 2024.
Thierry Lévy. Schur-weyl duality and the heat kereel 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, (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):2024331,, 2025.
Aaron Lou, Minkai Xu, Adam Farris, and Stefano Ermon. Scaling riemannian diffusion models. Advances in Neural Information Processing Systems, 36:80291 -00000 2023.
Yulong Lu and Wuzhe Xu. Generative downscaling of pde solvers with physics-guided diffusion models. Jourmal of Scientific Computing, 101(3):1-23, 2024.
Erwin Luesink and Oliver 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 wxXiv:2288.11977 2022.
14

---

Published as a conference paper at ICLR 2026
Andrew Majda, Ilya Timofeyev, and Eric Vanden-Eijnden Models for stochastic climate prediction. Proceedings of the National Academy of SSiencces 96626))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, 1088(2):1191144) 2014. doi: 10. (001080003001192220010
Remigijus Mikulevicius and Boris Rozovskii. Stochastic Savier-Stoker equations for turbulent flows.
SIAM Journal on Mathematical Analysis, 35(5):1250-1311 2004.
Robb J Muirhead. Aspects of multivariate statistical theory. John Wiley & Sons, 2009.
Philippe Naveau and Johan Segers.
Multivariate extreme value theory. arXiv preprint
vrXiv:2412.18477, 2024.
Bernt Oksendal. Stochastic differential equations. Spinger- Verlag, 1998.
Fátima Dalacios-Rodriguee Gwladys Toulemonde, Julie Carreau, and Thomas Opitz. General- ized 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. Hcavy-tailed diffusion models. arXiv preprint arXiv:2410.1417.1 2024.
2024.
Philippe Parnaudeau, Johan Carlier, Dominique Heitz, and Eric Lamballais. Experimental and numerical studies of the flow over 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 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 Engineering, 18:215-2611, 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:1 1631, 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-base data assimilation. Advances in Neural Information Processing Systems, 15:40521-40544 2023.
15

---

Published as a conference paper at ICLR 2026
Alex Rybchuk, Malik Hassanaly, Nicholas Hamilton, Paula Doubrawa, Mitchell J Fulton, and Luis A Martinez-Tossar Ensemble flow reconstruction in the atmospherie 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- diffusion model for multiplicative noise removal in texture-rich images and its fast explicit diffusion solving. BIT Numerical Mathematics, 1319)
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 arXlv:2510.02730, 2025.
Maneesh Kumar Singh, Joshua HooeeCColling Colin J Cotter, and Dan Crisan. Data assimilation using a global girsanov nudged particle filter. arXiv preprint arXiv:2507.1768; 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 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.
Stanislas Strasman, Antonio Ocello, Claire Boyer, Sylvain Le Corff, and Vincent Lemaire. An analysis of the noise schedule for score-bases based generative models. arXiv preprint arXiv:2402.04650, 2024.
Michael E. Taylor. Partial Differential Equations III: Nonlinear Parabolique Equations, volume 117 of Applied Mathematical Sciences. Springer- Verlag, New York, 2nd edition, 2011. ISBN 078-1-4419-7048-..
Surya T Tokdar, Sheng Jiang, and Erika I. Cunningham. Heavy-tailed density estimation. Journal of the American Statistical Association, 1 119(545): 2024.
Howard G Tucker. A generalization of the glivenko-cantelli theorem. The Annals of Mathematical Statistics, 30(3):828-830, 1959.
Michael Wilezek and Rudolf Friedrich. Dynamical origins for non- non-gaussian vorticity distributions in turbulent flows. Physical Review E E-Statistical, Nonlinear, and Soft Malter 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.
Yieun 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-26) 2600, 2014.
16

---

Published as a conference paper at ICLR 2026
Appendix
Table of Contents
Losses for score matching
18
B Stochastic Calculus and Stratonovich Integrals
19
B.1 Itō Integrals and SDES
19
B.2 Stratonovich Integrals and Chain rule
19
Conversion between It Stratonovich forms
20
Fokker Planck equation
20
C
Sampling from 1D) distributions
20
D
The Tokker-Plance equation and its invariant measures
21
DJ
Itō form of the forward SDE
21
D.2
Fokker-PP Planck equation Theorem
................................
22
D33
Distribution of the norms
23
D.4
Fokker Planck equation of the direction
25
D.5
Proof of Theorem 3.3.1 Convergence of Fokker Planck equation
27
D.6
Beyond pure Stratonovich noise
31
E
Latent distribution
31
E.1
The invariant measures define white noises the weak sense
31
E.2
Condition of Gaussianity for the latent vector
32
E.3
A tractable algorithm to sample latent vectors
33
Caussianization the latent vectors
34
shorter distance between latent data distribution
34
Relevance of MSGM latent space heavy-tail distributions,
35
F
Backward diffusion
38
G
Proof of Theorem 3.4.1 equivalence between ELBO and score matching
39
Statement of the theorem
39
Notations correspondence
40
G.3
Marginal density from Feynman- Kac representation
40
G.4
Change of measure and Jensen's inequality
40
G.5
Girsanov theorem
41
G.6
ELBO evaluation
41
From ELBO to our SSM Loss
42
Remark score parametrization
42
G.9
Cirsanov theorem in the transport noise literature
43
H Comparison with diffusions on Riemannian manifolds
43
Riemanmian Manifolds and Differentiation
................................
44
Conditional diffusions on scaled d-spheres
................................
45
Link with neural network architecture)
46
Analytic illustrations simplified cases
46
The two dimensional case
46
1.2
Tensor built from single skew-symmetric matrix
47
Non-commutativity the general case
52
Rank and skew-symmetry conditions for random tensor
52
17

---

Published as a conference paper at ICLR 2026
Proof of the rank condition
52
2
Tensor renormalization
53
Mean speed of convergence with renormalized lensor
53
K
Going beyond the rank condition for MSGM scalability
54
K.
1
Weaker assumptions
54
K.
2 Sparse tensors
56
3
Discussion about local and local structure
59
L
Numerical scheme
60
L.1
Numerical integration SDES
60
Scheduling
60
Loss evaluation
63
Neural network architecture
63
M
Details about numerical experiments
65
M.
Test cases
65
M.2
Data preprocessing
................................
65
M.3
Comparison strategy
66
M.4
Swiss roll
66
M.5
Anisotropic Gaussian distribution)
67
M.6
Multivariate Cauchy distribution
68
M.7
Vorticity field from Partiele Image Velocimetry measurements
Summarized comparison MSGM and SGM
88
N.I Theoretical aspects
88
Empirical aspects
88
A LOSSES FOR SCORE MATCHING
This section presents some classical score-matching losses. SGM backward SDE, the score V log (x) is replaced by a fitted neural network s0 (x, T minimizing some losses, like denoising, implicit, or slicing score-matching losses. Alternatively, one can work on 0g(x; (x, T 1) 2200( T -) (Huang et al. 2021). It leads to this SGM backward SDE:
t). This fitting is performed by
dw. x, 22(ao( l)dt + dB.),
(A.1)
A typical loss to learn this neural network is denoising score matching (DSM)
CDSM = ['' E (l (ael
(2V
(A.2)
where ||-|| is the Euclidian norm. By integration by part, we can show that DSM is equivalent to Implicit Score Matching (ISM) (Hyvärinen & Dayan 2005)
LISM (
E x. (1
1) V (V2a0 ds.
(A.3)
A reference V Vlogp. p. is not needed anymore. However, the divergence term may be untractable for large-dimensional problems. Using the Hutchingson trick, [Sng et al. (2020) shows that this loss equivalent trackable version the Sliced Score Matching (SSM)
LSSM 1'"
E- Ev-A(O, ( (2 s) + (v. ()(22o v)) do.
(A.4)
18

---

Published as a conference paper at ICLR 2026
Score matching is equivalent to maximizing the Evidence Lower Bound (ELBO) both in discrete time
(Lu
2022) and in continuous time (Huang et al. 2021). Indeed, denoting & the ELBO, Huang
al.
2021
shows that:
&o (a: = log Pol aT xo = a
(A.5)
f B 2 a ( + (Viag( s)) 20
ds.
The first term does not depends the neural network parameters 0. The expectation of the second term over æ following the dataset distribution is LIsM- So, maximizing the ELBO equivalent to minimize the ISM Table I of Huang et al. (021) summarizes the classical score matching losses.
B STOCHASTIC CALCULUS AND STRATONOVICH INTEGRALS
This appendix provides a concise overview of essential stochastic calculus concepts from Oksendal (1998); Kunita (1997) relevant to our work, especially the Stratonovich interpretation of stochastic differential equations (SDEs).
B. ITÓ INTEGRALS AND SDES
Let (2,FFP) be a probability space equipped with a filtration (Fi)ezo satisfying the usual conditions, and let (Bi)22 be standard m-ddimensional Brownian motion. Given adapted process X, Raxm satisfying appropriate integrability conditions, the Itô integral of X with respect to B is defined as the mean-square limit:
X, dB, = lim Σ
(B - B..),
(B.I)
[m|-0
where the sum is taken over a partition 1I of []. An SDE interpreted in the Itô sense reads:
dX, ((,,, di + ()),,))dB
(B.2)
where f f: Fdd R4 Rd is the drift, and G Rd x R. + Raxm is the diffusion coefficient.
B.2 STRATONOVICH INTEGRALS AND CHAIN RULE
Unlike the Itô integral, the Stratonovich integral is defined using a symmetric discretization:
|0 X. o dB. = lim Σ
X
X (Brui B..).
(B.3)
2
Bonneen
Stratonovich SDE written as:
dX, = (X,, dt o( dB,
(B.4)
key advantage of the Stratonovich formulation is that satisfies the classical chain rule. For any smooth function Ø Rd - R. we have:
dø V ) (s( Vx(X,) o)
(B.5)
Moreover, in multiscale deterministic or stochastic equations, if the fast component is 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 Stratonovich integral (Amoold] 1974). Note that is not always true for nonlinear dynamics (Gottwald & Melbourne 2013) Gottwald et al. 2015). Accordingly, the readers may interpret the Stratonovich noise Bd formal representation a process with a short correlation time that nevertheless respects the classical rules of differential caleulus, in particular, the chain rule.
19

---

Published as a conference paper at ICLR 2026
B.3
CONVERSION BETWEEN ITÓ AND STRATONOVICH FORMS
Given Stratonovich SDE, it is always possible to convert it to the equivalent Itō form:
dx.- fs(X 1 m E,,x,,)) (G,((, di + G)
(B.6)
j=1
where Ga( t) is the j-th column of the diffusion matrix G( 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, f(X,,t) dt t) dB,,
(B.7)
induces a time-evolution equation for the probability density p(æ, 1) of X,. This is known as the Fokker-Planck equation, given by:
0.pla, () - (f(x, - + v (V (E((,())((,(1)))
(B.8)
where Σ(x,t) G(x, t is the diffusion tensor. The Fokker Plakek-Planck equation describes the deterministic evolution of the probability density associated with the stochastic process.
C SAMPLING FROM ID DISTRIBUTIONS
Let us denote by PI. the distribution of the norms, || In MSGM, it is also the distribution of @ (see Proposition D.3.2). This distribution arbitrary, but a one-dimensional distribution. So, 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 regularized log-norm: log ||=]|e- := log(]|æ=] +e) with small, typically
-01 e 10 6 From a data set of the log- normnorms of M training samples, (log Fr Ile); = (log x Ile),
we define eCDF Fog l.l. as
Fing l.l : R [0. 1],
(C.1)
1 M R + M (R 1 II.)"
Σ
(C.2)
Then, we approximate the distribution of the norms, Plog (R)RR, by the empirical one, Piog (R)RR == Frog I-le (dR). In particular, we can sample a new norm of latent variables, || from a uniform one- nne-dimensional variable u ~ L4(0, 1) as follows
| = æoll = exp(log @olle) - = exp (Pime 1-l. (u))
(C.3)
(4)
cCDF is an efficient tool, but it cannot generalize the distribution 1j outside the training set (| II), For better generalization, instead, one could use a one-dimensiinal kernel density estimation or fitting of parametric distributions. In the of one-dimensional distributions with fat tails, classical kemel density estimation (KDE) suffers from bias in the tail estimation or peaks due to sparse data in the robust approaches that, in general, do not require the
tails. In that case, could consider existence of moments of the target distribution Tokdar et al. (2024).
In this paper, rely the eCDE.
20

---

Published as a conference paper at ICLR 2026
D THE FOKKER-PLANCK EQUATION AND ITS INVARIANT MEASURES
D.1 ITó FORM OF THE FORWARD SDE
Define the conditional noise covariance (xx) as
[(x) G( G(z))(z)) Bl(dā,)) ( , = x].
(D.1)
We begin with lemma.
Lemma D.1.1. Let the skew-symmetry assumption Al hold. Then,
d
1 V (E(x))
G k Fa.
(D.2)
2
2
Proof. Let explicitly state the matrix divergence. For k I..... d define (x)x) (x)] : c= then we decompose Σ as
[x)) G(z))T Σ (x).
(D.3)
Then, taking the divergence of the j-th column of Σ* {x)), we obtain
(z)()) v (x)
Σ
(z))
(D.4)
Ozi
1=1
d
Σ
G C
(D.5)
Oeg
p.g=1
Σ
Giy (ōs + xxp0ig) 2
(D.6)
p
d
Σ
G
P G G x.
(D.7)
p,4=1
[tr(G*
G + G (G*x))
(D.8)
By skew-symmetry trace( G*) = 0 and consequently
d
. (I(x))
(V
(x)))
= Σ G*
1=1
k=
Lemma D.1.2. (Forward SDE Itō) Let the skew-symmetry assumption equation A7 hold. Then, the Itô form of the forward SDE equation 3./ of z, is given by
dæ ( () ds G( .dB
(D.9)
Proof. Using the standard Stratonovich-to--t- formula (e.g. Kunita 1997. it holds
dē. = Gd B.). + dd...
(D.10)
1/4 Σ G H.. (B.)). + dB
(D.11)
1 Σ
2 2 ds G(d.)dB.
(D.12)
( {) ds + G x dB
(D.13)
where the last equality comes from Lemma D.I.T
21

---

Published as a conference paper at ICLR 2026
D.2 BORKER-PLANCE EQUATION AND THEOREM 3.1.1
Let mo ~ P2. with Pæ. E C² (TR.) be twice continuously differentiable. Let ps denote the density
of the distribution of x. For each E TRd we define the normalized vector xx" == | for +
and 0 otherwise, which is orthogonal to the d-sphere S Furthermore, let V be the orthogonal projection of the gradient V in the tangent plane x: Ta Sd-d the tangent space on the Riemannian manifold gd-1 the point æ defined for f E C² (TR:)
Li = vf - (x" f) x".
(D.14)
Lemma D.2.1. It holds for any smooth vector field f that
V (xr () f) V
(D.15)
Proof. Let f be smooth vector field, then
2
(a"("")
x (2r" a V)Si (lx", v)))
(D.16)
holds for each j L,.... d that
(x" V)="11
Σ
(x=")
(x"))
1
Σ
(x))
I2
(D.17)
Ox:
|
||>||
||
i=
xx//||<11
Σ
//
(D.18)
1=1
0.
(D.19)
Consequently,
(a:"
(x f] = (x" V)r"" f (x.". V)S x= =(x" V) x="
(D.20)
Using the decomposition = x (x""")) (I (>"' ="'x") + V1, equation D.16 and equation D.20 we conclude
V
ʃ
x x
V.) (x" v)S x" V f f- VV)"" 0 + V S.
Define the conditional noise covariance [(x) as
{xx) = G G(z) = El(dz (dz.) T = xx]
(D.21)
We can now state and proof Theorem 3.1.1
Theorem D.2.1. Let the assumptions A7 and A2 hold. Then, the Fokker-Plance equation of equa- tion 3.7 reads
{
a (x)) V1 (1 L(x))))
or e ER.
(D.22)
Po
P
Moreover, any stationary density Poe equuation D.22 is rolation-invariant on RRd
Proof. From the Itō SDE equation D.9 the Fokker-Planck equation describing the evolution of (ps) (p=)=20 is given as
0s
0
Ps
V
2
v (E()))
v. (E()))((a))
(D.23)
V.
2 V. (E(z))p.(z)
v (E())) + (
(D.24)
V.
(D.25)
22

---

Published as a conference paper at ICLR 2026
The skew-symmetry condition in assumption A1 implies for any R that
((=") E())x"
Ga (x (T") 0.
(D.26)
k=1
Combining this with the result of equation D.16
it holds that
V (F())Vp(<z) = (x" <= [(x)Vps(=)) + V1 Vpx))p (D.27)
The decomposition V ="(x") (1 x" ("")) ""(x")<V V and equation D.26 yields
__ (E(x))V)) = V (E(x)V + (E(x)=" (x') (D.28)
V (P(x)V
(D.29)
Hence, by linearity
8 P. V. v 1 VL
(D.30)
0s
We shall now explore the of possible invariant densities Pro the Fokker-Planchl equation D.22 We will show that Poo is stationary if and only if it is rotation-invariant.
Let Poo rotation-invariant, i.e. P = then is a stationary solution of the Fokker- Plankkk equation. The set of rotation-invariant measures is not empty, e.g. containing the isotropic normal distributions (0, Is).
Conversely, let Pao be stationary solution of the Fokker-Planck equation, in particular
V ( (+)) (x)( = O.
(D.31)
Integrating over the test function & Poo gives a necessary condition for pso be invariant measure:
Po Po
dx,
(D.32)
V V ( E(x)V PPo
(D.33)
Rd
=
G (x) (x)
_P ViPo (x) dæ.
(D.34)
Hence, for a.e. æ E Rd it holds that ViPao (x) E ker k (x)). By assumption A2 this kernel has dimension 1. Moreover, GT (æ) and by definition of V we have that V Poo L.. That means
Peo (a)) ker (T)))
span n = {0).
(D.35)
We conclude that V ViPoo (x) = almost everywhere on RR i.e. the measure is rotation-invarian
D.3 DISTRIBUTION OF THE NORMS
In this section, along the noising process whereas in SGM the norm dynamics is random with mean going to va asymptotically.
precisely that the nomn of the MSGM SDE solution remains constant
D.3.1 NORM DYNAMICS IN SGM
The dynamics of the SGM diffusion norm is stochastic. The following proposition states that the norm the SGM latent concentrates around its mean, Vd, for large dimension d.
Proposition D.3.1. If F. is an Ornstein Uhlenbeck process then
E
|
1 20
2s 1 + -29)d - d.
(D.36)
83000
and
Vd v 01))
E211 E
[
+ Vdl, + =
(D.37)
with both EK3 and T.TT2 bounded for large time and TET? independent the dimension d.
23

---

Published as a conference paper at ICLR 2026
Proof. To get the dynamics of the squared norm mean in SGM, we can take the expectation of the following Itō equation
d d| 2τ. dæ.
+ o , > 2( 1|, d)ds + 2V22, dB., V 0. (D.38)
Thus,
|| c - + - e 20) d. E
(D.39)
80000
To obtain the full nom dynamics from equation D.38 we note that + 28 has finite variations. Accordingly,
d(e (( d)) = 2/ 2a a. dBB..
(D.40)
and temporal integration and the analytic expression of the Ornstein Uhlenbeck process yields:
/-" 22(0-5) .d.
11.1 = E
11 T + +2V2
dB
(D.41)
+ 2/22 -x0 1 "(-o--')
= E 1
+2v2 1/¹ 10- (2000-''-0 dB.. dBs.dB....
(D.42)
|| 20
vd 1. te "K.,
(D.43)
with the martingales
88 1//0 / e -2(2--5'-8") 8") dBs dB
I. Val
(D.44)
K, = /8 t0 / c dd.
(D.45)
K₂ corresponds the martingale part of the Ornstein Uhlenbeck solution projected on 2g. is well known that EK? is bounded for large time EI? may be less known and we shall evaluate it below:
1.
E12 =
d
E
Σ
e't d(Br), d(
(D.46)
8
E
Σ
1 est d( d(Ba)) d(ss)) d(Bsg)rs
d
P1:P2
(D.47)
88 1 s + pp 5101 - 52) 5(8"
8
=
d
4s
E
Σ
ds"
JU
0
P1:P2
(D.48)
8c 48 10 2 ds'ds"
(D.49)
4e 40 (( 22d (620 -1dds',
(D.50)
4s 4s 1) + 2( - 1)),
(D.51)
(D.52)
84000
24

---

Published as a conference paper at ICLR 2026
D.3.2 NORM DYNAMICS IN MSGM
For MSGM, the norm follows totally different dynamics. We recall that the skew-symmelry of odZ, implies that dæ, = odZ,æ v, is orthogonal to æ and hence:
d 1 2x, odæ. = O, A3
(D.53)
Consequently, 2., moves randomly on | | the d-sphere of radius Toll. and the increments dæ, are tangent to the d-sphere. In particular, we obtain the following result.
Proposition D.3.2. Let the skew-symmetry assumption hold. Let æo be random variable. Then, for all 2 0 the distribution of | equals the distribution of ||00|----
Therefore, the distribution of the norms of the latent variable is exactly the distribution the norms of the points of the dataset. Moreover, T. = if and only 20 = 0. As a consequence, we can exclude all points exactly equal zero from dataset, treat them aside, and hence consider, without loss of generality, that T f 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 Sd-1 i.e. the distribution of & in particular as s + Do.
=
D.4.1 MAIN RESULTS ON THE DISTRIBUTION OF DIRECTIONS
We in Section D.3 that Ex moves randomly the d-sphere of radius || and that the increments, dē. G x ) sB,,, are tangent to the d-sphere. If the rank condition, assumption verified, then the support of the noise distribution G( dB, coincides with the d' 1-dimensional tangent space, i.. it will likely explore all local directions around z.. With time, the support of the = solution distribution will gradually cover the whole d-sphere, i.e. every direction v, will become equiprobable. Lemma D.4.I illustrats and precises this claim.
Lemma D.4.1. Let assumptions Al and hold. Let a initial density pe E (Zd())))) and E(x") G( G(x")), Then, the Fokker-Planne equation
05 "=" = V VI''")VIPP' V
x" E
(D.54)
has a unique density solution p" E (2( all 0. Moreover, there is a unique invariant measure p'a that Fokker-Plaane equation, i.e. the uniform distribution on the d-sphere USS111) with density
Poo
(xx"
Vxr" gd-1
(D.55)
with gd-1 2n 4/2 /T (g) the volume of the d-sphere 8d-1 and I the gamma function.
Lemma D.4.1 is a consequence of Theorem 3.1.1 as shown in Section D.4.2 Note that in this case V
ga- the Riemannian gradient on
see Section H.1
Given the unique invariant measure of Fokker- Planck equation formulated on gd-1, we can also show that we have exponential convergence of the initial distribution p to ph the uniform distribution on the unit sphere gd-1
Theorem D.4.1. Let assumptions and hold. Then, there exists a a G, d) with
||p". p 11 I exp(-as)]lp ph 21 11200--13
(D.56)
The convergence rate a given as aG, d) = 1) min (m.)
S (( v) s x gd-1 l
(D.57)
Consequently, since 1d-1 is compact this implies convergence in total variation of ph to ph and convergence in distribution of v, with xo u(d-11). The full proof is detailed in Section D.4.3
25

---

Published as a conference paper at ICLR 2026
D.4.2
PROOF OF LEMMA D.4.
I.4.
Proof. Existence and Uniqueness:
2 (x)). L is a parabolic type operator according to Friedman 1964) since x: + (ar) is positive definite by assumption A2 on gd-1, Indeed, for any E T T the tangential (lincar) space of Sd-1 at x,
Consider L(ph) = V - Va)VVii
x y 2 0. y' E()yy= G
(D.58)
with equality if and only if G (z)y = 0. Then, the rank condition A2 implies y = 0 as previously in equation D.35 Consequently, the associated spatial operator Lo defined by
Lop" VT ( (L=(<)V (x))
(D.59)
an elliptic operator on Sd-1, compact manifold without boundary such that the semi-group e is strongly continuous on (Sd-)))), 8 2 O. As pe E (S(S--1), according to chapter 1, proposition
1.1 in Taylor (2011). there exists unique solution pd C2 (2( for [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
>
Invariant measure: Repeating the lines in the proof of Theorem 3.1.] that pic is distribution.
given in Section D.2 it follows The only rotation-invane distribution on the d-sphere the uniform
D.4.3 PROOF OF THEOREM D.4.I DIRECTION
LIMIT BEHAVIIR OF POKKER-PLANCE EQUATION OF THE
Theorem D.4.2. Let assumptions A7 and A2
hold. Then, there exists a a d) with
llp. Poc 11:2000-) I exp(-as)]lp P'c d2ll? 11/2
(D.60)
The convergence rate given as a(G. d) = (d-
min (x.W)s
S= 5=((x, y) E s x s (lyy). x
(D.61)
+ Proof. Let p. denoting the density of E. Define ea p" 1o with ph = 1- being the uniform distribution on gd-1 Then, by linearity of the Fokker-Planer equation, e, satisfies
d V V (x))
(D.62)
Since y and ph are densities 14-1 we have fga-1 dæ for all 2 Consequently, since gd-1 is compact manifold without boundary, Poincaré inequality holds, i..
1 |
| 112 1) S d-
ga- e laa ">
(D.63)
with
V Et (ÿ)
y=u Projy, VE((z) y=
= (z)
(D.64)
Consequently, integration by part on gd-1 leads to
1 d
| dl
n 2 12
(34-1)
=
Iga-1 E (x(x
dæ,
(D.65)
2
Ige-1 =-fa- V
(D.66)
We will now bound y [x)yyT from below for any y E x and x E Sd-1, in particular with a bound independent of æc. Since Σ(x) is symmetric, it is real diagonalizable with eigen-basis denoted as vi (æ), va(x) E g and eigenvalues (z)...... >x(x). By construction [xx)
0, hence we can set va(x) := x /||æ=|| and ^a(x) = Moreover, by the rank condition A2 A o 0 for o d. By orthonormality of the eigenvectors en (æ),..., vd-1 (x) then span the tangent plane x at æ on
sd For any 1.....dd. .d we have that
G GT(z)y/) [y]
(x) = (, [(x)),()))) = GT |
(.))
min
(D.67)
26

---

Published as a conference paper at ICLR 2026
The polynomial (a. y)
P(x.3) 3) (G yl) the compact S {(xx)) y)
gd-1 x
g x x-y) attains its minimum P", which from the rank condition satisfies P(x. y) 2 p > 0 for
all (æ, 3) S. As consequence >x(x) p for = I...., d and a E sd-1, which implies for all E x- that
y Σ x
GT
(=) P* 11ll
(D.68)
Therefore, combining equation D.63
equation D.64 and equation D.66 we obtain
lef l1 (80--) $ - - V_ e 1 "_ve"lll
& - 1'
(D.69)
2 dt Then, by Gronwall for a p*(d- 1) > we conclude that
lp -Po 1880-11 ||e"| 12 I |10 110
exp(-ash
(D.70)
D.5 PROOF OF THEOREM 3.3.1 CONVERGENCE FOKKER-PLANCK EQUATION
This section is devoted to the analysis of the Fokker- Planckk equation in the whole domain R 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 Σ (z) for 0, we exclude the origin in the analysis.
Theorem D.5.1. Let D Rd (0) for d 1. Let assumptions Al an hold. Let zo po (D) and let lId be the (radial) density of @ol. Then, the Bokker-Plandee equation
0 0s Pa(x) V (LL((())VVLL()))
MED. D,
(D.71)
has unique solution ps the stationary distribution
(D) n (D) for all 0. Moreover the Fokker-Planck equation has
Pl. (1| 1
Poo (x)
(D.72)
|
Sd-1
In particular. || and Φ, are asymptotically independent for + to0. Moreover, there exists aGG d) > 0 such that
lpp. - - (22 S eep((-as)]lp Poo 12
The convergence rate a given as
a(G,d) d) (d (d-1) min (x.
GT
S ((x, y) E sd-1 gd g
(D.73)
Proof. the following, we will frequently use the change of variables in polar coordinates in IR". More precisely, writing x rw" with r > 0 and xt E S"-1, the Lebesgue measure decomposes as
dæ n- dr do (x=",
(D.74)
surface measure on Sa-1 This change of variables is justified by the change-of-variables theorem; see Folland (1999 Theorem 2.49) and the discussion on polar coordinates.
where da(x") = dx denotes the
We will proof existence, uniqueness, regularity, invariant property and convergence separately.
||| r) be the start value of the FP equation D.54 on gd-1 of Lemma D.4.1 This gives rise to a smooth unique density solution p" (a>" | = r) for & > 0 and any r > 0.
Existence: Let (o( [|]"|]]]]]
Moreover, p" (z -0ll r) p. (z -= r) since d -. 0. Now define
(= 2. (z: क [|]]0)),
27

---

Published as a conference paper at ICLR 2026
where we denote x" ||||| Ps is density since
Ps (z) dæ /
(x") 1 ([|=]]) ||||| da,
(D.75)
24
Ps (x
||| ||
= r)m. (r), dr da"""
(D.76)
R4
(x ||| 2 =. (r) dr dæ".
(D.77)
B 84
pl)(r))
(x | = r) dx" dr,
(D.78)
R+
1.
(D.79)
We have V V gd-1 and V does not act on radial functions. Besides, E(x) = G (a: (
| with G linear so [xx) = [|||ك| |
= | (x=)) Hence
(L)((
|| ((===)) = (l]( --- L æ=1l = ||a=])),
V (ax क..
(D.80)
(D.81)
[= -'d [ V p= || || ||=))),
(D.82)
and
1
V(x)V = (())p) [|)=||
Vga 1p" (x² ||
(x"||
|||x]]]
(D.83)
[
|,>|>>>><>>>
Jp
Vge-ip"
[ ||=||)
(D.84)
(1==))),
E(x"
- ||=||| (D.85)
|
TA (E(x)VLp,(t)
x Vga-
|
[))
(D.86)
(1
|) x - V gd-1
Ve="") (x: E.1l --==-))
(D.87)
2
(x
| = ||>>|| )
(D.88)
Os
2 i
(p.
(|
| ||=||)
(D.89)
Os
θ
2 0s (x)),
(D.90)
i.e. p.(()) V (E V (=())) Then, Ps solves the Fokker-Plance equation on Rd
Uniqueness: Assume there exists another ` solving the FP on Rd verify
(x(x) (x. x ((Ix=))
Since d|| = 0. by marginalizing p. (integrating on Sd-1), we have the uniqueness of the radial
density (r) = -sa-- Ps (dæ"" = P| (r)r--d = p). (r)r
28

---

Published as a conference paper at ICLR 2026
Since dl| 2 0, we have 0
0s P2, ([|=])) 0. Therefore,
Os Pe(x) = 0s
(pl (2)
(D.91)
0
2
(1ael))))
Ml. (x" | + p1.,
Øs (2.0(()))))
(D.92)
05
0
(1.x(1
Pl
[=|))
(D.93)
Os
In addition,
0
Os (,(x))
VT [E (x)(
(D.94)
V1 (E
(D.95)
|2.<(|==|||<=|| VLLL"\ [=)))
(D.96)
[2)0([= (E(x") V [x)))
(D.97)
and finally,
(1)0()))))) 2
as P1 (1 (x="((===11)] Vga-1 (E(x") Vga="VVa (L.
0. (D.98)
Then, Pl (1.a(x" |||=]| is solution the de Fokker-Planck equation on the sphere, for any x such that = {r Rt (r) > 0}. If || | A. then 02 (|]+ = 0 and
||=|| E ==
Px(x) P. (x"1| = > ppp (1=
For a
E A. then (1 * and (-1+(=1|==]] is solution of the Fokker Planck equa-
tion D.54
According lemma D.4.I] the Fokker Planck equation [D.54] has unique density
solution p" (ax" moll ||=|]). Hence, for any æ such that (()=l))))) >> we have
(x = p. (x: [1l | = | (2e))))1( (|==
is true for any æ. So, Ps (x: the unique solution of the Fokker Planck in Rd..
Regularity: By definition of the marginal density, have
/ gii
/
Po Pr. d) do d-1
Po(
8d
with (rr, 0) = r0. According to the assumption 2o ~ Po E C² (Rd {0}) and compactness $d-1 conclude that P|. ( ()
Since [| = r) is C², we have that p. (x"
(x"
x
= r) is smooth by Lemma D.4. for any
s > Consequently, Pa (xx) smooth on D for any 0.
Invariant distribution: The distribution
PIi (1
Poo (x)
Sd-1j
radial function in || The operator V does not act on radial functions and 3 is in the kernel
of V such that V )
Hence
0
Is (1)
Vi
(E (x)), D,
(D.99)
Pl. (||==|)
1
VI.
(D.100)
|
gd-1
(1 Pl. |0
= V.
V1
(D.101)
| d-
|
(D.102)
29

---

Published as a conference paper at ICLR 2026
Therefore, the Fokker-Planck distribution Poo is the stationary. In addition, Poo is a density since
P1i
(1
Pro (az) da
dx,
(D.103)
x d-
[Sd-1]
11
Rd
Ra
PL. (r)
p
dr da""
(D.104)
pd
Sd-1|
R4 ga-1
1
(r.)(r)
dr dx
(D.105)
Sd-1
R. g4-1
/ P..()d dr
1 |dd-1
dx"
(D.106)
R+
=
(D.107)
Convergence:
Hence, obtain for p. - Ps that, we can bound the speed of convergence
P/i ((22
Ipp. - Poo (001122(2))
P=(2"" || = -|==||pp] [|<=||]]]||]
dæ.
2 1
|
Bd
(D.108)
/
Pl.
xe 11)
p" (x=" [ol==las [[|=||<]]))
[ 2dd-2 dæ,
(D.109)
g
2
(r)
(s | =r)
d0
pd-1 dr,
(D.110)
Sdd-1 1
R3
Pl. (r))
1
Po( - || =r) [Sd
pdd dr,
(D.111)
12(-1-1)
R4
|poco || = r)
(r)
< eepp(-as
rd-1 dr,
1.2 gd
(D.112)
/
1
p
2c 1)
exp(-ass)
Po (zx" || -|<==])
day
(D.113)
Sd-11
ac
2dd2
Rd
/
P/. (1|=0|
exp(-aas)
Po (a0(==1] 20ll ((
dx.
xx
S
(D.114)
exp(-as))1pp Poo Plolli..eo 1
(D.115)
where in the inequality used Theorem D.4.2 The upper bound is finite since Poc E (RR)). In 1-d
order to see this, we will show that the function (r) is oo). Recall that
())(()) = /
(r( 0)) do pd- /
Po (o(((r))) de.
3d
gd
30

---

Published as a conference paper at ICLR 2026
Then, application of Jensen inequality leads
/
pd-1 /
p
P-1(7)2
dr
pr( 0)) d0
nr-ddr, dr.
(D.116)
R
R
84-1
yd- dr, (o( 0)) do
(D.117)
x S
Po( ()) 2))2, d0 dr,
(D.118)
R4 9--1
Sd
1| (.(2)))))
(D.119)
Consequently,
Pl. [1+ |])
PPo l1 (RR))
(D.120)
||2|| 1
12(29)
(r)
drd0 v [gd-1
p ( < co.
(D.121)
1d-1
pd-
8d 14
D.6 BEYOND PURE STRATONOVICH NOISE
A possible extension to the described diffusion in equation 3.1 would be to add drift term, i.c. considering
dæ Aa ds + G( dB.,
with skew-symmetric matrix A E Bddd Then, the associated Fokker-Planck equations will additional involve an advection term (Aæ) V Psr which can used to improve the speed of convergence of the dynamics.
E LATENT DISTRIBUTION
The latent vectors 20 Poo 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 Gaussian. We also show that the MSGM latent distribution is always closer than the SGM latent distribution to the data distribution. Finally, focus the case of Cauchy data distribution, where SGM leads to singularity, unlike MSGM.
E. THE INVARTANT MEASURES DEFINE WHITE NOISES IN THE WEAK SENSE
In additive SGM, latent vectors E. ~ Poo are Gaussian white noise in the strong sense, i.e. for any i + j. the coordinates (#oc); ) and (#oo) are centered, independent, and identically distributed. In contrast, the latent vectors MSGM are white noises in the weak sense, as stated by the following proposition. For any i # j. the coordinates (#xo) and (x o0 are uncorrelated but neither Gaussian nor independent, in general.
Proposition E.1.I. Let the assumptions Al and A2 hold, El 2 < +oc, and Pro be stationary density of the Fokker- Planckk equation D.22 Then, Foo Peo white noise in the weak sense, i.e. Ex x0 = 0, E(x .))? < +00 independent i, and E ( (o) ( = Vi.i E [l....,dl with i * j.
31

---

Published as a conference paper at ICLR 2026
Proof. From Theorem D.2.1 Poo is totation-variant So there exist function h RR" IR such that for any a E Rd, Poo x
h[[[=]]) Then,
B( = / zgh da -J- Sen. (L x1 I = 0,
(E.1)
since the function 24 - 2. hv || ||] even.
Similarly, for i + in {1....,. d}
S me maach((1aal)))
E B(Ix0)
(E.2)
Jse. In (f sch(l doseddo
(E.3)
0.
(E.4)
Moreover, for in {1,.... dd. we have
+00 >
E| 30
(E.5)
E 20)
(E.6)
E(
(E.7)
[ (1e (f (fat((l=e))d
(E.8)
which does not depends
Remark 1. Since E[ >0) does not depend on i, easily evaluate it from Theorem 3.3.1 and
Proposition D.3.2
E (x
11 112
lim
21 d E 13
(F.9)
d
and thus
1 2 La.
E0 [x 00 xo)
(E.10)
Therefore, monitoring the covariance of T. and its distance ooll are convenient proxies of the forward SDE convergence.
E.2 CONDITION OF GAUSSIANITY FOR THE LATENT VECTOR
Proposition E.2.1. Let assumptions Al and A2 hold and P1,22 be the density of @0ll.. Then, the latent distribution Poo is Gaussian if and only if P/, is a scaled x² distribution with degrees of freedom, denoted a² X2 with a > O.
we know that Poo is rotation invariant, i.e. a function ||=]]] If this distribution is Gaussian, it has to be of the form (0,0, a Ia) with > 0. Then, || || || a² x2- Reciprocally, if there exists >> such that P1,12 a x²x3 then PJ = axa. where we denote by axa the distribution of positive random variable X = aa R such that R x3. From Theorem 3.3.1 we know that
Proof. From Theorem D.2.1
1-d
| 2
Pro (x) = PP. [[)))
Paxa
(1]=|])
|
(E.11)
is the distribution (0,00 a²I written in spherical form. So. the latent distribution P.o Gaussian.
32

---

Published as a conference paper at ICLR 2026
Remark 2. Isotropic Gaussian data Ty N, 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 po
po(x)) = 100 ||0||
||a||||
PI. (1a1(11a)
||0|||
||
(E.12)
||==||
The latent distribution would be Gaussian as long as the distribution of norn is P1-1 = Paxa But the conditional distribution of direction can be any valid conditional distribution on the d-sphere. For instance,
)
Pó
|
-
with e = 01,0,....00)
(E.13)
a:
is a valid candidate even though is not Gaussian (since its support is R x {03d-1).
E.3 TRACTABLE ALGORITHM SAMPLE LATENT VECTORS
With the following proposition, if know the distribution of norms, PI. we can sample latent vectors from Poo.
-N Proposition E.3.1. Let To ((0, and wo F-1 (z)) with (x) == f([[==]) + 0 and O otherwise,
ø if
Tan
(T)==FTT
Vr >r>0
(E.14)
for the generalized inverse CDF F. PI. and Fx2 is the CDF of the x² distribution with d degrees of freedom. Then 2o Poo.
= Proof. Since Zo ~ NN(, Ia), we know that || -N [1
-- Do 1| x²-1. i.c. F (12) ~ u(0, 1) and
-N
-N then R f(] To II) FTj' FFya
P.I- In addition, the normalized vector is zo and the direction
11a)
U(s-11). The norm
are independent. Therefore, R and follows the correct distribution.
are also independent. We can conclude that To
In practice, do exactly know the distribution of the data norm PI. So, we do have access
to Fj.1 f. Instead, approximate the distribution of log || @ with ||æ]|e = |
+ e, denoted
Plog I.l,, by a model Prog I.l., or equivalently Flog by model Fiog (see Section similar sampling procedure for the latent vectors, replacing F by our approximatin.. We obtain
We perform
samples of an approximate latent distribution p.x. as stated by Proposition E.3.2
-N Proposition E.3.2. Let To (o, and to = = ( x + 0 and otherwise,
with (x)(x) = f([[=]) T if
f(r) : exp (Fiet (Fa(aa) -, Vr 0,
(E.15)
for the generalized inverse of the approximated CDF FooRlil. associated the approximated PDF Piog Iol and F-2 the CDF of the x² distribution with d degrees of freedom. Then To ~ Poos where
|| 1-d
Poo is the empirical approximation of Poe. that is Poc = Prog I-1. (log xxlle)
Va E TRd.
ga
Proof. We can follow the same proof that for Proposition replacing FJ.1- f. P/.,- and Pso by Flogl.1- l. Prog lile and Poo respectively.
33

---

Published as a conference paper at ICLR 2026
E.4 GAUSSIANIZATION OF THE LATENT VECTORS
If needed, we can easily build second latent space with standard Gaussian vectors. As stated by the following proposition, for any (non-Gaussian) latent vector Ty, we can create a Gaussian vector T
-.V
N TT = Ry with -" ///1/ and R = f-
(E.16)
N If T is zero, we just set 2 zero.
Proposition E.4.1. Let TT Pea, where Poo the empirical approximation of Poo, that Poo (2:) 2
N Vx E 12, and Fr F (TT) with F(x) ([|=])
Pioglele (log [ and 0 otherwise,
if x +
|
(7())
V (Flogl1
Vr > O,
(B.17)
for the approximated CDF FloR|-le associated to the approximated PDF pion l. and Fxa is the CDF of the x2 distribution with d degrees of freedom. Then Ey (o).
-
Proof. We can follow the same proof that for Proposition E.3.1 replacing FiT f. N2, Ia), PM and Poo by Fing l.le f Piog le Poo. x", and (0, Ia) respectively.
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 the SGM latent distribution in KL-divergence.
Proposition F.5.1. Let the assumptions Al and A2 hold, P1,2 be the density of |0oll 2 Pso and pod (0, Ia) be the MSGM and the SGM latent distributions respectively, then
DrI(peellpo) & DRI(PNllpo)
(E.18)
with equality if and only PIP is a x² distribution with d degrees of freedom.
Proof. We recall that the MSGM latent pdf is
(z) = P1.1 (||==1) 1 | Sd
(E.19)
and the data distribution reads
pa( PJ. (1|==1) d- Po (a=
|||||=0||| ||
(E.20)
Let denotes Px3 the x² distribution with d degrees of freedom
Px3 (11==)) CN (x) PON
x c Pa |1 = ||=:|])-
(F.21)
CCN F (20) with F(x) := f- ((|=|) | i A and otherwise, and
It the distribution of 2u
f-((r) = (F (Fi,(r)), (FF((F))())))
0,
(E.22)
and (R) SSo (r) the edf associated PI.-
34

---

Published as a conference paper at ICLR 2026
We have
0 < (xI.polllo"" )
(E.23)
/ po(2)) log
po
dx.
(E.24)
Po
/ po(x)) log
P|.1 <1|=>
[]
(|x=1) dæ,
(E.25)
Px2
= / po(a)) log
PI. (||==])
ae u 1-
dxc
(F.26)
||
[8d-1|
Pxd (1x=1)
/ po(2r) log
Poo (a)
dac
(B.27)
1x (a:)
/ po(x)) log
Po (z)
Poc. (a:
dx,
(E.28)
1::0 (x) (x
/ po(2)) log
Po (z)
po(x)) Peo (x)
dæ
po( log
dx,
(E.29)
1x (z)
D PKi(Poll) (12))
DK (xIl(PolPPa))
(E.30)
with equality if and only if pa PN i... P|.
Px3
E.6 RELEVANCE OF MSGM LATENT SPACE FOR HEAVY-TAIL DISTRIBUTIONS.
This appendix provides 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 B.5.1 In particular the derived inequality in Proposition E.5. becomes meaning less both sides are not finite. However, as we will see for example of heavy tail distribution such as the (product) Cauchy distribution, this not the case. To this end we will show in Section 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 Dection/E.6.2.
We note that the analysis can be extended to broader class of heavy tailed distributions and more general SGM latent spaces such as general Gaussian distributions.
E. INFINITE KL DIVERGENCE BETWEEN CAUCHY DISTRIBUTION AND STANDARD GAUSSIAN
Let d(x) /2w
22 be the density of the standard Gaussian A I). and let po( πI m(1+2²)
be the product density of univariate Cauchy distributions. Then, the following holds. Lemma E.
(Ko(pollo) =
Proof. Let L > and define the set
M {w
E ddd
21 >
L. 1111
j= d).
Then for M and == x --2d-I
C 1
(o( 71 + xt)
nd
1
1+1
n 12d-
1 1 x7
1+
1=2
Moreover, for any Rd, it holds that
(af
E
2/1/2
(xx) = (2T)
2/2
122
I (2x)--122
4/2 =2//2
Consequently, for L large enough,
po(
po(x) log
>
+ Olog
((x)
2
35

---

Published as a conference paper at ICLR 2026
where O refers Landau-symbor of big-O notation. Together, for M and L large enough, there exists C > 0 such that
po( log (a)
C> 0.
(bx))
2
Consequently,
/ po(2)) log
Po(x) (Φx)
/
dx 2
Cd ... dzd
dzi = cc.
1ky151,522
By Lebesgue decomposition, we conclude DRI 00.
E.6.2 FINITE KI. DIVERGENCE BETWEEN CAUCHY DISTRIBUTION AND ITS RELATED foo.
Let 2 and again consider the product of Cauchy densities
po(
II
TT x?) ==1
Let 2% Pp and let pR be the density of R equation 3.8 consider the density
@0ol Then, motivated by latent space distribution
Pn
Pn( (+
Poo (x)
E Rd (O))
(E.31)
|
1d-11
Then, the following holds Lemma E.6.2.
DKL(pollPoo) < c0.
Proof. It holds that
(x Poo (x)
log
lo)))( - log (o])x()) + (d - 1) log |||||||||| + log S
Hence,
DKL(Pollpso) llog (l.llo@@0)) po(2))] [l, ( ((]l + ((d 1)1E [ log (
where E denotes the expectation with respect to probability measure Pndæ. We will show, that each term separately is finite. We start with the first term, followed by the the third. The finiteness of the second term turns out to be consequence of the finiteness of the second term.
First termiIt holds that
log (20)) -d lo
log(1 1 (20)3)
for Z0 Po- Since coordinates (#o); are iid is enough check the marginal integrals for finiteness. In particular it holds that
log 11 + x²) dx log(4) < co.
m1 x2)
Consequently
E [olooooo20lll) < xc.
(E.32)
36

---

Published as a conference paper at ICLR 2026
Third temm: For the second term, let R || and M ((zo)d. Then, almost
surely
M MSRS VaM
7 log M & log R & log M + log d.
2
Consequently,
Dog RJ - E, MJ] - log d.
Thus if we show Epllog M] then lpllog R. < oo as well, since both expectation only differ up to a finite factor. Using the CDF of (Io)1 e.g. for (20)) holds that
P(
F(-1) arctan t, 82 0.
Consequently, since ( (Fo) are iid, the CDF FM of M satisfies
Fxr(t) = P(M < = 2 arccant) t 1 0.
Hence, the density Sas of M is given (for d > 2) as
(t = ot (t) = d (2 (arctant) tjd-1 1
1 12
Now, by integral splitting we find that
ac
/ log log (t)dt + log(t),
Epl M] log logt (f)
(E.33)
By noting that for & t. (t) s Otd-1 for some C > 0 and
[M(0) log())) & - (-
1 σ² a> a>0,
the first integrant of equation E.33 is finite using a = d. For t 2 arctan r/2 and hence (t) < C" 1+2² for some C' O and the second integral of equation E.33 is finite since
[ log(t)
1 dt = 1 < oc.
It follows that [oglog M] is finite.
Second term: Recall that
pn(r)= = f
(do(() = dd1 / (o()
gd-1
g
Since for r0 with 0 (0 ......., (a....,,),, with 0 < 1 and using the fact that
1
po(r0)
( + 22)
conclude
(r
क 2
Td + r2)d
Therefore,
-
(R( 2 d-1
18
=: C2
nd(1 + 2
[1 + 72yd
37

---

Published as a conference paper at ICLR 2026
Hence,
logpR(r) 2 log Ca (d - logr d d + +
which yields
EpllogpR(1) [oll]] 2 log Ca (d - 1) [log RJ dB, log(1 RR)].
(E.34)
d
For the third term in equation E.34 it holds that 1t2 = 2 2o) Now for M S 1 we have
since Rt < ddM < vd that log + 1t²) log(1 vd) independent of R. For M 2 log(1 + Ii²) & log(1 + dM²2) I log(dM² + dM²) = log(2d) + 2 log [M). Since we already showed that Ep log R) is finite, we conclude that the lower bound equation E.34 is finite. For the upper bound, note that
po(r0) S Td
Vr > O.
Thus,
PR(P) I yd-1 |3d-1| =: Card-1 nd
for some Ca > O. So
log logpn(r) S log Ca + (d - 1) log
And finally, since Epllog R| 0c
Ep Blg [R) < log Cd + (d - 1) log R 00
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 symmeryyyyrry assumption A7 hold. Then, the Itô form of the reverse SDE associated to the forward SDE 3.7 given the SDE
da, (V (z, + G G( x
logpr- Pr V
G(z,)dB
(F.I)
In the Stratonovich form, it reads:
dx. G( ( ( + odB,
(F22
Proof. From Anderson (182)); Song et al. (202T) and the forward SDE (see Lemma D.1.2 we
know that the Itō reverse SDE with negative ds writes
dz = (V )) (s - (V D) G(2.) TV log p. (2.)ds
+G v dB
(F3)
(V 3)( 1logg + G( ((d.)dB, :
(F.4)
Replacing the decreasing E [0. T by Brownian motion B, we obtain the Itô backward SDE with positive dt
with increasing E [0, T and using another
da. (V [) G G( G ( ) TV log ( + G G( dBB
(F.5)
38

---

Published as a conference paper at ICLR 2026
Then, Lemma D.I.1 and the standard Stratonovich-to-Itô formula (e.g. Kunita 1997) Stratonovich form of the backward SDE:
yields the
dæ. = d (.... B.): + VV E) G(a,) TV log (.)
+ G( dB:,
(F.6)
- Σ) (t (V E)(a,)
+G( G (V dtt odBi)
(F.7)
G
logpr- odB
(F.8)
Proposition F.2. (Backward probability flow ODE) Let the skew-symmetry assumption Al hold. Then, the reverse probability flow associated to the forward SDE is given by the ODE
dæ:
2 G( log loggr-rrtee pr
(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æ. 1 (V () - (V (d) G ( G G TV logp.( dtt, (F.10)
G((
G G TV Io (
(F.11)
Replacing the decreasing E [0. T) by T with increasing t E [0, and using another Brownian motion B, obtain the Itô backward SDE with positive dt
dæ =
(F.12) 12)
L
G PROOF OF THEOREM 3.4.I: EQUIVALENCE BETWEEN ELBO AND SCORE MATCHING
This appendix derives score-matching soure-matching-basee ELBO for MSGM training. In this work, we focus the simple forward multiplicative SDE equation 3.1 Nevertheless, we here derive slightly more general theorem, where include a possibly non-zero Stratonovich drift fs.
G.1 STATEMENT OF THE THEOREM
Theorem G.I.I. Let us consider the forward SDE
dæ = fs( fs( + G( dB.,
(G.I)
where assumption A7 holds. Then, have
(= > (x)
(G.2)
with the following ELBO
Exc = E log pr
(epp(((1)) (T) 20 x
(G.3)
E t. fla0 s) V. (G(2.) 8)) fs fs( 2o ds.
Proof. Here, review the work of Huang et al. (2021 SGM and generalize some of their results to derive an ELBO and justify score matching for MSGM. Note that Benton et al. (0044) : Ren et al. (2025) proposes very general SGM framework with associated ELBO and score matching losses. The MSGM ELBO and thus the above theorem can be understood as particular case of their work. The explicit dependence in 0 is omitted for readability.
39

---

Published as a conference paper at ICLR 2026
G.2 NOTATIONS CORRESPONDENCE
The forward and backward processes are denote Y. and X, in Huang et al. (2021) and æ; and w; in this paper. The forward Itō equation Huang al. 2022 is denoted:
dY, = f(Y,, s)ds + gY s)dB,
(G.4)
Lemma D.1.1 gives the forward Itō equation MSGM. It yields the following notation correspon- dence:
g = G(æ),
(G.5)
= E(()) D(x) =3(())
(G.6)
f(x) 1(V [T) fs(Ø).
(G.7)
And the backward equation is:
dæ, = ((z + G
(G.8)
with drift
(ex = ff(x)) + ( P(g- = -ss(xx) + 2 V (x)
(G.9)
(G.10)
where we would arrive at the approximate backward SDE of Figure if we replace V log (w,)
by so a T T-a also parametrized ag G We note that case, f&z 0, the drift reads (V Σ) + EV log Pr- 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.
202T treats the general of multiplicative noise. It states that
Po(x))
E
pr xT
((Fr)ss exp
(-V I( s) v 1( (T, T T s)]ds -ss)dss 2o
(G.11)
T [f" V (μ 2V E) T -s)ds) 200
E
pr
exp
(G.12)
where
dz
-i T s)ds G( s)dB'
(G.13)
jix,t) u(x,t) - (V [x))
(G.14)
and B', is Brownian motion.
Remark our case, pixx,) p(x, (V ETT E(z) 1 [V
-- E)T is twice the Itō to Stratonovich correction of the backward SDE equation G.8 (see Lemma D.I.2. 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. 1.14 verified by &s including twice the Iōō to Stratonovich correction.
Remark 4. Note that -= -
G.4 CHANGE OF MEASURE AND JENSEN'S ' INEQUALITY
From the Feynman-Kac representation equation G.12 and Jensen's inequality, obtain ELBO as in Huang et al. 202T).
Let (2,F, F. P) be the underlying probability space for which B' is Brownian motion. Suppose Q is another probability measure on (52, F) equivalent to IP (i.c., they have the same measure zero sets). We can hence apply the change-of-measure
(-1" V. (u )(V E).) T -s)ds) 2o
dIP dQ
(z(
E
P1 2z)) exp
(G.15)
40

---

Published as a conference paper at ICLR 2026
Then, apply Jensen's inequality:
dF dQ +
/-¹
-T V (n 2V [ T s)ds 2 a
Po E
log
log pr (TT)
g00
(G.16)
Compared to Huang et al. 2022 we have the additional term (V Σ)T that is, -2 Σ)
G.5
GIRSANOV THEOREM
Huang et al. 2022 apply the Girsanov theorem to the following forward SDE equation (17) of Huang
al. (2021)):
dæ. = (- Ga ds GdB., 5.
(G.17)
since the Itō Stratonovich correction 1 (m (i) ( E) zero in Huang et (2021).
However, it is not the case in MSGM and here we use the Girsanov theorem to this forward SDE instead:
dæ = (-je + Ga)ds - GdB
(G.18)
The Girsanov theorem Oksendal 1998 solving
Theorem 8.6.3) states the following. Let B an Itō process
dB. -aw, dB',
(G.
for w E l and Bu = 0 where satisfies the Novikov's condition. Then B is Brownian motion
with respect to Q and
E L s(w,s)- -d', -
- (/ ||aww, |lalw xo = x (G.20)
E
log
Co = x
dQ
En [lalw, 8) x ar] ds,
(G.21)
since 7 ~ 1o' s(w,s) dB'', is martingale and thus E [ aaw, 8) dB' -0 Oksendal 1998
Theorem 3.2.1).
G.6 ELBO EVALUATION
Equation G.2l enable to evaluate the ELBO 8o given by equation G.16 To evaluate the divergence term, we note that:
(4 - (V (x)
-fs(x) + 1 (V (x) + E(x)so(x. - -+ (V-E)
(G.22)
ff(()) + (Z))
(G.23)
Then, the ELBO simplifies to:
dP
(x) =
E
log
2o = æ
+ log pT xr) 2o x
(G.24)
E.
(u [ Σ)T 2 ds.
log Pr
xT) zu x:
(G.25)
E m
flla0( s) + V (G fs((z)) =
x ds.
41

---

Published as a conference paper at ICLR 2026
We recall that in our case, fs cancels out. The first term E log pr (T)) æ is a constant w.r.t.
to O. So, if when maximizing the ELBDO, this term has effect on the optimization procedure. Therefore, even with our multiplicative noise, the minimization of the ELBO corresponds precisely to Implicit Score Matching (ISM), which 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 Al and zero Stratonovich drift, i.c. f. 0. The theorem states the
(x 0) > (2) == C(xo) - Lo (xo10)
(G.26)
with C being a constant with respect to the parameters be learned. More precisely,
C(ro))
E log (o)pp(()))))) (T))
(G.27)
Loo (2)o0()) /'' R [ s)1 V ( (G( fs(x))) 2o = ds.
(G.28)
Then, average the data Zo to obtain the following lower bound for the likelihood of the dataset:
Pag Eino & ( () Rez C(ar)) E.g, fox (2ro
(G.29)
Our objective is to find the neural network parameters 0, that try maximize the likelihood of the data set, Eo, Po (0). Since Exo ((zo) a constant with respect to 0, we maximize - -Ez,, (20 0) Let us explicit the terms above with the Hutchinson trick, Evw [VVT] Ia (Song al. 2020
Ego C(xo)
= Egr logpr PT (T))
(G.30)
Exg Coo (zol0)
E_ t
ao
fs T.. s) V. -)
ds,
(G.31)
//" E.
2 a
a 3] v (G
ag a (G(
8)) 7 do,
(G.32)
E
aaooo. x. .) + (G as s.) s))
(G.33)
= TE
| a )) + (Ry~Rad(d) G [. x))
(G.34)
= TE Ead(d)))) [2laa(E,, 5)1l + (v. G) s))
(G.35)
= (0).
(G.36)
Therefore, maximizing the ELBO. Eyng Eso Ex (mo
F. C(xo) Euno Loo (ro) (0), equivalent to
minimizing our practical score-matching loss, (fsM(
G.8 REMARK ON THE SCORE PARAMETRIZATION
Following Huang et al, 2021 we directly model G G( log (x) by neural network ag (æ s). If needed, the projected score, V Ps, be retrieved directly from ag as shown below. Note that the full score,
V log = V logp, + (x=
loggss Ps:
(G.37)
involves radial term, (x=") log that cannot be directly estimated in MSGM.
42

---

Published as a conference paper at ICLR 2026
Proposition G.1. We assume that assumptions A7 and A2 hold, and that we have approximation, ag, of the scaled score and an orthonormal basis u2 ...., ud ua( ) of x that we concatenate in (x) = [u)())..............).......))) u. ud(()) E (gdx(d-1) Then,
(m) E (1a))
(G.38)
approximates the projected score
U (oggs)))
(G.39)
Proof. Since Rd -Rx" 2 x we have La = ===="<=")> obtain
UU)) + Using "())x" = 0, we
UT (T(T)(()))((Z
(G.40)
= VTo
(G.41)
(G.42)
= (Taa)
(G.43)
(a)) ( is full rank,
UT V U
(G.44)
It is also possible to model the score, V (. directly by neural network, so 8) using the following acore-matching loss:
LssM(O)
=
E5-20/0,, Rud(d)
f1lao(E... (.. s) (v V)( -) v/.
(G.45)
= Ee Ev-Rad(d)
[ [1ls0( x(.)) + (t fs( v) ds.
(G.46)
However, for any a E (0) LssM ax"), i.e. our loss function insensible to the radial component of the score (xr" . V) log ps. Therefore, our MSGM framework does not provide estimation for the radial score (x* V) log ps- Moreover, the optimization problem parametrized by 89 is ill-defined, and the loss should probably be regularized follows:
SSM (0)
(Ss + t. #, [(xx.s))3]
(G.47)
with 7 > O large, says = 106.
G.9 GIRSANOV THEOREM IN THE TRANSPORT NOISE LITERATURE
Following work done by Huang al. 2021 for additive noise, we have relied on the Girsanov theorem OKsendall 1998) 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 al. (2020a 2023) González et al. (025); Singh 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 case, the weighted score,
ag (z t). biases the noise, /BB/dt, in our backward SDE make its solution closer the forward SDE solution (see equation 3.12). This noise change is the core of Girsanov theorem (see equation p. Resseguier (2023) also proposed fit 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 2d and SGMs on manifolds. To introduce the subject, 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.
43

---

Published as a conference paper at ICLR 2026
H.1 RIEMANNIAN MANIFOLDS AND DIFFERENTIATION
This section is devoted to a brief introduction to Riemannian manifolds and the associated differential calculus. For more comprehensive discussion, refer to Lee 2018). Let M be smooth dimensional embedded submanifold of Rd, where S d. For any a: E M we denote by T.M the tangential (linear) space of M at æ. We denote by a Riemannian metric on M, which assigns to each a E M an inner product
9N : T.M T.M- R.
In the case a smooth embedded manifold in the Euclidean space, the induced metric is given by (z, v) (u, v)R=, Vu, E M.
This makes gM, Riemannian manifold. Let {e("), ..... e(n)) be orthonormal basis of T. Then, the orthogonal projection onto T.M is the linear operator Pz : RR + T.M that satisfies
(v) arg min ||u w|le Σ (c, (0.e
wET.M
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 Fd E smooth, its Riemannian gradient can be computed as
VMI(2) (a)))
where ((()) is the Euclidean gradient. Furthermore, we want to define the Riemannian divergence in this framework. For tangent vector field f . + Rd with f(x: E T&M, the Riemannian divergence is given as
div ( = (0 (())
f=1
where 0at f is the Euclidean directional derivative. Finally, the Laplace-Beltrand operator AM can be defined as
AMf dV.)m(V
which generalizes the Laplacian to M.
In the special case that M = rgd-1, 1 for radius > 0 then n = d 1 and
T.M- = T.rsa-1 = for R (u, x) = 0},
and (v( = 2a (u, Since x= x we obtain (v( (I xx"x") ((==")))r y and as result
-uf(=))PP((f(aa) (1 x² (z) (f) vif(x).
(H.I)
Regarding the Riemannian divergence, we note that æ", (1) (n)) defines an orthonormal basis of Rd. By Lemma D.2.1
V f(x) (x" (x f(x)) + V f(x).
For f(x) E T,M, we have f(x) x. 0. Thus:
f(x) V f(xr) (xr" (f((()) xx) V f(x))
Differentiating the tangency condition f(x:) xx" = along xx leads
dr (f(z) x²) (0)) + (f(x), 0g"
Since Oz"x" 0, we conclude that f(x))
0. Finally, expanding f(x) in
x" e(1) e(() leads to
V f(xx)
V f(x) = Σ (0 f(x)
(0) + (0x) e()) = div (f(())
(H.2)
im
In our setting Im (m(E())) -ي TyM. Hence, the right-hand side of the Fokker-Plande equa- tion
3.4
div (V(( = V (F())V[(()),
(H.3)
generalizes the notion of fovergence--orm operator to the manifold setup.
44

---

Published as a conference paper at ICLR 2026
H.2 CONDITIONAL DIFFUSIONS ON SCALED d- d-SPHERES
Several authors have recently developed SGM on Riemannian manifolds (De Bortoli et al. 2022 Huang al. 2022 Benton al. 2024 in order generate data lying a particular manifold.
Clearly different, our goal is more classical: generating data IR However, each solution path
of our forward and backward SDE lies on its scaled d-sphere | [Sd-1 Clearly, d-spheres are particular cases of Riemannian manifolds and possibly the most studied. De Bortoli et al. (2022)
describes diffusions of æ in the d-sphere gd-1, The simplest the d-sphere that converges to the uniform distribution on the d-sphere, P Unfortunately, this appealing proposal does not directly apply to framework: the Brownian motion on the d-sphere
involves a Brownian motion
not a solution of our forward SDE of æ Indeed, in general, there exists æ =
E gdd1 such that
(x)))
Ex k=
G (0""")(e (G* + Igs-t. So, the Fokker-Planck equation of De Bortoli et al.
2022
i px") = divga- (Vs--p(" Vx" gd-1
(H.4)
Ds
and our Fokker-Plance equation for the direction,
0s p" (xx) = divga-: (E="" Va ES
(H.5)
do not coincide. However, the analyses from the SGM-on-manifon community on the finite-time distribution, its score, approximations, and score-matching losses choices could certainly facilities the MSGM training process in the future.
In our case, the norm of solution being constant along path, write both the forward and the backward equations of the direction on the unit d-sphere from equation 3.1 and equation
dæ; dæ: a
G G( dB.,
(H.6)
1
G(Zz)
(G( lo (1 dt + odB.),
(H.7)
xell
G G(
G G(z"")
(I log (11=,1(zz,) di + od
(H.8)
We note that | _ = Vga-1 Oen is the Riemannian gradient on the scaled d-sphere || Therefore, using p.o. the density of the couple of variables 11, "²") RR+ x sd-1,
[ log -r (1x,,1(z2) )2""
|| t-11
log pT
(H.9)
log
PT-
[1:02 |1,
(-1 1111-
(H.10)
Or
n log
PT-c
x
||
(1
|
dxm
(H.11)
i
log 1
(671
[
(H.
i
log logp7_.
(+n n 22
||
(H.13)
dr
Vga
logpt
(H.14)
and finally
dæ
G(
G Φ )
TV log
nn xx
||
di + odB:
(H.15)
In contrast, forward and backward SDES of De Bortoli et al. (2022) read
1d
n dæ
dB,
(H.16)
gd
(-1) dl + B,B [1
dæ
= V gi log togp"-c
(H.17)
45

---

Published as a conference paper at ICLR 2026
where B, and B
Brownian motions the d-sphere. They be defined from Stroock's
representation
Hsu 2002 Example 3.3.2) as
gd-1
21 ( oB dB, = -Ia
(H.18)
dB, = (Ia (a, (1, dB,.
(H.19)
The first main difference with MSGM is that the projection on the tangent plane, (La (4a- (x")), =n (quadratic in æ") is replaced our approach by G x (linear x"). Accordingly the noise
(conditional) covariance, (Ia (x") = (Ia (a" (a=))) (projection property), is replaced
f-12 by G G G az, ) E(x,) ). To make our diffusion coincide with equation H.16 we would have consider
G( = | (Ia xx (x="
(H.20)
which is Lipchitz continuous but nonlinear. As such, the noise covariance would be
Σ(x)
|
(H.21)
In general, we can hardly expect such simple form from MSGM noise covariance. However, for the random tensor equation 6.1 we can show (see equation 0.10) that:
2E ax)
|| Ia- .x..
(H.22)
In addition, our score involved in the backward SDE equation H.15 depends on the norm || The norm || | appears covariable with prior distribution PI.1 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,
better understand how the direction and magnitude 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 || ||8-11 the conditional score direction will be oriented differently, pushing along some orientations on some spheres and along other directions spheres of larger radius. Accordingly, along the backward diffusion, the directions tend align differently on different hypershperes. The distribution of direction become more and more radius- dependent.
+ If data samples æn snapshots of conservative dynamical system, all data points probably have the similar energy 20. P, i.e. Var(E)/EB22 is small. All data points are on closed scaled d-spheres VE gd-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
G ( logPr- (zz)) = G(
T log 1o
6- ||
(H.23)
justifying our neural network spherical architecture equation
L.33
G G( TV Volgrr- T - āe
xr log
TTolle.T
1 .eelle:
(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, corresponds to the Brownian motion the circle.
46

---

Published as a conference paper at ICLR 2026
Let us recall this forward SDE: (z
K G (t) dB, (t) dB, *k
Σ
G dB, z(t).
(I.1)
k=1
In dimension 2.
dz((t aJæ(() dB,,
(1.2)
where x=
11 2
E RR? J
(
0
-1
is the -ootation.
1
d
I
= a
22
o dB t
(1.3)
22
I1
+0
Then, in the complex plane,
= 1 + 132 E C, with i V-1 and we have:
C
(1) = aix +0 +1
dz
od
(1.4)
since
11
dž + = 2 2z2) dB, = a0-22+ + tiž dB,
(1.5)
The solution is the Brownian motion on the circle:
t1 (0) explaiBB)
I (1) = x
(1.6)
i.e.
R(aB.)=(0)- = cos(aB,) -1 sin(aB,))
sin(aB.)
2(1)
x(0),
(1.7)
+1
cos cos(aB;))
i.e.
-6 (1(1) (0) cos(aB,) T0)(0) si
-1
(B.,). -1
(I.8)
=1 (t) = (0) sa + (0) cos (...
(I.9)
The key element of the proof was the possibility write the forward diffusion with single skew- symmetric matrix in equation I.2 Below we generalize this idea to larger dimension > 2.
1.2 TENSOR BUILT FROM A SINGLE SKEW-SYMMETEE MATRIX
Here assume that whole tensor G is built from the same dense skew-symmetric matrix G i.e.
G² GO', vk E {I.....d} d),
(I.10)
skew-symmetric matrix. As explained in Section K.2.1 this tensor respect the condition Al not the A2 Nevertheless, this case and its analytic solution may be insightful.
with G'
1.2.1
MATRIX EXPONENTIAL
Here the full Brownian matrix can simply factorized as
Z. Σ c (B.)x G² Σ (B.)..
(I.1I)
It has the same distribution than
= V
(I.12)
with B another single Brownian motion. The forward diffusion simplify to -
dē =d x, dB
(I.13)
with solution
v, = = = exp (VdG'E.) 20.
(I.14)
since Z and dZ', commute.
47

---

Published as a conference paper at ICLR 2026
1.2.2 DIAGONALIZATION IN THE COMPLEX PLANE
Gh has pure imaginary eigenvalues and be diagonalized in C orthonormal basis
G = UCCiAAUU
(1.15)
with Uc complex unitary matrix, real diagonal matrix, and the superscript H denotes the conjugate transpose. Then, the solution can be easily evaluate as follow
v. = exp (ivanD.) U æ0.
(I.16)
For an even dimension d, and for j E d/2]), there exists X3 E I such that
exp (V (Ut 20)23-1
(U
(I.17)
(U 2.)2j exp
(-ivax;B.)
( 2o)2..
(I.18)
For an odd dimension d, Gt has least one zero eigenvalue. Without loss of generality, we consider A and for all j E {I,..., (d - 1)/2), there exists Xj E RR such that
(U x 22j21
exp
(ivā (v
(U zol2j-1,
(I.19)
(UL 2.)25 exp
(. (U 2o)2j.
(1.20)
(U Dz)a (U wo)a.
(I.21)
1,2.3 REAL SOLUTION WITH SINE AND COSINE
The diagonalization matrix, U, complex but we can find a real unitary matrix, Up. to make G, block diagonal, and then expressing the solution with cosinus and sinus in equation L8 and equation] L.
(U 22)-1111 = cos (Vd (U (Uf - sin (Vax,B') (UT
(I.22)
(U sin (V (UH + cos (vaxςD)) (UZ
(I.23)
For an odd dimension the real solution reads
- sin (vaxzB') (Ug
(f cos Vd x, (UT
(I.24)
(U sin (V (Ug +oo (vaxzE')) (Ug
(I.25)
(U = (o)
(I.26)
Figure illustrates the solution for with 20000 realizations of ØT large time 7 100,
with A 12 10, 2o = (1, 1, 1), and UR I4. A rotation-invariant distribution, Poo. would induce rotation-invarian 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 1.10 cannot guaranty rotation-inven inntiraninvariantt latent distribution, and thus prevent the use of our simple eCDF-based based sampling procedure.
Figure 7 illustrates the latent vector support for random initial condition Zo ((111, I, 1), 0.0144). The supporting manifold is not one-dimensional anymore, but still depend on the initial direction distribution, r8.
48

---

Published as a conference paper at ICLR 2026
1
0
x
4011
+2
X1
MWM
M 0
1
1-
X2
-1
0
1
-1
X1
X1
X1
X4
Figure 5: Projection of samples, Tr. sketching the support of the invariant measure, Pso. for low-rank tensor [.10. d = 4 and Zo 1, 1). The top plot is in space #2, 23),, the bottom plots are, form left right, space (I1, 22)) (x1 23), (x1,444. and (I4,, 43).
49

---

Published as a conference paper at ICLR 2026
Moreover, as expected from the expression above, the initial norm, | scales the one-dimensional manifold supporting the invariant measure (not shown) and the initial direction, F has influence at large time. 6 shows the same example with Zo = (V2, v 0, 0). The initial norm the same but the initial direction is different. Therefore, the limit distribution, Pae if it exists does depend = on initial direction, Tg. making the latent sampling intractable.
x
0
-1
2
-2
0
-2
X1
2
0
*
0
x
-2
-2
-2
-
-2 0
-2
-2 0
-2
X1
X1
Figure 6: Projection of samples, ær, sketching the support of the invariant measure, Poo. for low-rank tensor 1.10 d 4 and T (v2, / 0, 0). The top plot is in space (21,42,27)) the bottom plots are, form left to right, in space (#1. 22), (I1, 1z). (#1, 24), and (I4.23).
50

---

Published as a conference paper at ICLR 2026
1
0
-1
-2
2
0
-1
+2
X1
1
1
1
X2
0
0
0
0
-1
0
- OOR
-
-
X1
X1
X1
X4
Figure 7: Projection of samples, Tr, sketching the support of the invariant measure, Pso. for low-rank tensor 1.10 d = 4 and Zo 1,, 0.01 The top plot is in space (I1, aT)) the bottom plots are, form left to right, in space (21 22), (21, a3), (æ1, 24), and 23).
51

---

Published as a conference paper at ICLR 2026
1.3 NON-COMMUTATIVIT IN THE GENERAL CASE
For general tensor G in dimension d > 2, it is temping to look for a solution æ. of the forward SDE
dæ od at
of the form exp( Z.) 2o with Z - -E (B.)ss. Z. being skew-symmetric, exp( Z.) is unitary and such a solution would be reminiscent of the rotation form of equation I.tiiion[_Jjam and equation I.16 derived above. However, exp( (Z is not solution of equation 3.1 in general, since dZ,Z, & Z.dZ,
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. We will show that this tensor respects both assumptions [A and A2 almost surely. Then, we will discuss the speed of contraction of the Bokker-PPance equation with this tensor.
J.1 PROOF OF THE RANK CONDITION
Proposition J.1. Let M E Rdd be iid random matrices with entries drawn independently from (0, 1). Define the skew-synmetric matrices G (M* (M)) and R (0) define
the (random) matrix
G(x) (G' GZar.... G E Rd.d
Then, almost surely rank (ann(((()))) = d - 1.
Proof. Let x + O. Let M be random standard Gaussian matrix. Then, let D = M MT. Then, D Gaussian matrix with entries drawn from (0, 2), in particular
Did M[M M, Max)] 2(6118ee 8zedza).
(I.1)
Consequently,
( MT M M)] -E(Da D).
(1.2)
Now, for the covariance structure it holds
[B[Dax D)),s
E((Da))
(J.3)
B(D DRere)
(J.4)
T B ( Dae.:
(J.5)
t=1
2
1 81e05x):
(J.6)
d
x eeg -2 2 d
28
Σ
2
2
(J.7)
#=1 1=1
1=1
20 || 2x 1
2r.
(J.8)
Hence D(Dxx D) = 2[[1=] La T. Consequently for any k = 1..... d it holds that
1 E[M' (M*)) Jax ((M*)) -M 2/2 (1 xx.)).
1
3G x)
(J.9)
As result, the matrix G( has columns
G
itd ~
N(0, V) with
V = V(x) = =EE(x) = 2 (||=| I xx.);
(J.10)
52

---

Published as a conference paper at ICLR 2026
of rank d - 1. Therefore, ((x)
(a(V(()), is a Wishart matrix.
Let K = (x) matrix
K K Kd-1],
(J.11)
with column vectors K; forming an orthonormal basis of the hyperplane xt. Then by construction, have
2
KVK = Id-1.
(J.12)
2
This means, that
L
KG(x) = |
KGK æ...
v2 KG"
L2 KG (O,,oa-))
(I.13)
-
|=
|=
Therefore.
<K(T) 2
V2 30
K
= (1m KG()) (1a) W
(J.14)
Wishart matrix, in particular (C n) denotes the Wishart distribution with degrees of freedom. In the case n 2 such matrix is invertible almost surely (Muirhead 2009 Theorem 3.1.4). In our case d >> p d thus almost surely
rank
22 KG())) |
rank (1 (G()))) d
(I.15)
Now, since G(x))xx = we obtain almost surely that
d rank
V2 KG(x)) s rank(G(())) <d [m]]
(J.16)
which yields the claim.
J.2 TENSOR RENORMALIZATION
In practice, renormalize the tensor G follows:
Vd G with G = (M; My.))
G
(J.17)
|G||2
The normalization ensures that the trace of matrix defining the Itô term of our forward SDE i.c. the term driving the exponential decreases of EE (see the forward Itô SDE equation D.9
tr 1/4 E -r I )
1IG = GIG = = 1d,, (J.18)
-
&
similarly the trace the matrix defining the Itô term of classical Ornstein Uhlenbeck forward SDE
tr (-) = -d.
(J.19)
This normalization helps to better control the speed of convergence of the forward SDE without changing its keew-symmetry nor its rank.
J.3 MEAN SPEED OF CONVERGENCE WITH RENORMALIZED TENSOR
Note that in this case, for (z, y) E (( E x gd-1|=
d E| G(*)))) = ) y Σ (ah (G = y. d
" xx d
2 (J.20)
k=1
So, we can expect exponential convergence the Bokker-Plancee equation with the speed
Blo(G, d.) (d - 11) min
[GT(x)vl = (d 1) T 2ol) / d(d 1). (J.21)
(m.g)eS
53

---

Published as a conference paper at ICLR 2026
Therefore, the convergence gets faster when the dimension increases.
However, the tensor G is normalized (see equation J.17 so the evaluation of the convergence speed is modified. We note first that:
d
E||G|| = d d ==
M - TTP
-(B(M))2-EEMM (1.22)
d Σ (1
1 -
k3
(d - 1). 2
(I.23)
2
So, we obtain an estimate by Sauchy-Schwarte and Jensen's 's inequality
vd
d
(G))))))
E
G())a|
G())))l
(1.24)
GG|
GG|
E
E
(J.25)
||G||2
d
(x)))l)
(J.26)
E G|
2
d/2
(J.27)
ddd - 1)/2'
(J.28)
1'
and finally we obtain the following bound
Ea(G, d) < (d- (d-1))||(TT2))) 1.
(J.29)
K GOING BEYOND THE RANK CONDITION FOR MSGM SCALABILITT
The dense tensor of Section J imposes a computational complexity Odd). 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 weaker set of assumptions.
K.1 WEAKER ASSUMPTIONS
We recall here the two main assumptions of the paper
Skew-symmetry. For any & E {L....,d}, d), the matrix G² (G is skew-symmetric. (A1)
Rank condition: For any E d ((0), rank(G(æx) = d 1.
(A2)
Note that the Fokker-Plaan Planck equation 3.4 Proposition [.2] Proposition 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 [.1.I] Theorem D.4.1] and Theorem 3.3.I rely the restrictive assumption A2 and unfortunatcly, it seems difficult to find a sparse tensor 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 E Rdd only. This set of assumptions will yield a definition of a sparse tensor in Section K.2.2 providing satisfactory numerical results in practice.
Skew-symmetry. For any k (1....,d d), the matrix G
(i, is skew-symmetric.
(A1)
Rank condition almost everywhere: For almost all æ E Rdd rank (rakkGG(:)) G(x:) = d 1.
(A3)
The assumption A3 means the set AG fo E Rd G(x)) < - 1} has zero Lebesgue measure, i.c. SAq dr = Obviously, the assumption A2 implies the assumption A3
54

---

Published as a conference paper at ICLR 2026
The right-hand side of the Fokker-Planck equation 3.4 is function of V ViPs only. Hence, under the weaker assumptions [A] and rotational invariant distributions are still invariant measures of the Fokker-Planck equation. Following the proof of Theorem D.2.I] saw that the invariant densities, Pac, are characterized by [G 1 (=)] = O almost surely. So, if ((x) has rank d- 1 for almost all with respect the Lebesgue measure, then this requires V LP (2:) 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 guaranties the invariant distribution (Theorem D.4.I) need more careful analysis. We only sketch some challenges involved, the full analysis will be carried out in follow up work.
If the diffusion process enters the area of points æ, such that rank rank(((()) < d 1, one has to make sure that the diffusion process is not trapped in such an area, even if it has measure zero. In particular, let D C gd-1 be defined as
D (æ" gd-1 rankkG((=") 13--
Then, we call D a trap set. If the process V. once ever enters D with positive probability, it cannot leave again, i.c.
= F(#. E D, Vs 2 80 2 E D) >
Hence, in this case, convergence the correct invariant measure has ensure that the trap set is not invariant under the biffusionn-contrllldd dynamic. Such a analysis then sufficiently can be implied by Hörmander/ bracket-geeerrtiing conditions, i.e. hypoellipticity analysis. Based on this, asymptotic results Theorem (31.0)mmm..1111 [31.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, also Section K.3 for related discussion.
Although detailed analysis of this research question is out of the scope of the current manuscript, want to stress its relevance related to the scalability of the proposed method for the high-dimensional case. The standard construction via 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 Section K.2.3] satisfy the strong rank condition, the local sparse tensors from Section K.2.2 only satisfy the weak rank conditions. Still, the latter have been applied in our numerical investigation for the nigh-dimensional test cases with Particle Image Velocimetry measurements as discussed in Section M.7 yielding first very promising results.
only
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 Section K.2.1 will lead us to consider this of assumptions insufficient for the MSGM sampling procedure.
Skew-symmetry. For any l (1, (1....,d) d), the matrix G (G),)),, )i,j is skew-symmetric.
(A1)
1 d
Σ (G (G*)T is full rank
Itô term rank condition: the matrix = 2
(A4)
From Lemma D.I.I note that S -V Σ. 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 Al is verified. Taking the expectation of the Itô diffusion, we get:
d Ex, ds
SEE.
(K.1)
Instead of controlling the convergence of the full distribution ps. 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 A7 holds. Then, the following assertions are equivalent
The assumption A4 holds.
55

---

Published as a conference paper at ICLR 2026
E x Øol 0.
5-200
Var x
Tul 6-000
Proof. S is positive semi-definite, so it is diagonalizable in an orthonormal basis, and from equa-
tion K.1
Ex - 0 if and only if S is positive definite, i.c. the assumption A4 is verified.
50000
Besides, by assumption A1
the norm | conserved along the diffusion,
Var
[
| 2:
|zoll |B,b,|20||""
(K.2)
which converges to | || if and only if E [F0] 0.
We highlight the fact that the assumption A4 is weaker than assumption A2 as stated by Proposi-
tion K.1.2 Is is actually strictly weaker assumption since the tensors defined in Section K.2.1 and
Section 2.2 respect assumption A4 but not assumption A2
Proposition K.1.2. Let the assumption
A7 holds. Then, the assumption A2 implies the assumption
A4
Proof.
If the assumption A2 holds, then, Theorem 3.3.1 implies that 2.
00 = -| Fo00
The asymptotic latent direction, oo is independent of the initial condition Do and has zero mean. Therefore.
E 2 0
E
[
B ooo m
[E|
[
| E [ 0, (K.3)
and by Proposition K.I.I
the assumption A4 holds.
K.2 SPARSE TENSORS
Here we propose several possible choices of sparse tensors.
First, we will consider simple low-rank tensor in Section K.2. and show that it makes the latent distribution untractable. Then, we will introduce sparse local tensor in Equation (K.8), which is adapted to MSGM and leads to good generative skills in practice. Finally, we propose sparse nonlocal tensor in Equation (7.17) that involves more Brownian motions but meets the original assumptions Al and A2 of our paper.
K.2.1 LOW-RANK TENSOR
simple choice of tensor with d² ( non-zero coefficients would be take d times the same dense random skew-symmetric matrix G i.e.
G M
G (M,, M),),
(K.4)
std (0,
(K.5)
Section I.2 provides an analytic solution for the forward diffusion in this case. Such a solution would be strong advantage for our learning procedure, bypassing the need for numerical integration of the forward diffusion, and enabling denoising score matching instead of sliced score matching. However, Proposition I below shows that there is rank deficiency, probably inducing the existence of rotation-invariaan latent distribution, Poe. preventing MSGM sampling in practice. Indeed, numerically illustrated in dimension d = 4, the analytic solution of Section shows latent distribution intractable in practice. The latent distribution is not notation-invarinn and does depends on the initial direction distribution, P. It seems to be direct consequence of the rank deficiency.
We conclude that low low-rank tensors as in equation K.4 are not a suitable choice for MSGM. Moreover, it suggests that assumptions AT and A4 as in Section K.1.2 are not sufficient for MSGM.
56

---

Published as a conference paper at ICLR 2026
Proposition K.2.1. If G is defined from equation K.4 and equation K.S rank(G
then, for any R. < 1. Assumption [A] is verified, assumptions [2 and A3 are not for d > 3, and
assumption A4 is verified almost surely if and only if the dimension d even. Moreover have
- (G' and did-1) ld.
Proof. The tensor defined by equation K.4 and equation K.5 obviously matches the skew-symmetric condition AT
For odd dimension d. G - like all skew-symmetric matrix is singular. Thus S is singular and even the weak condition A4 is not satisfied.
For even d the polynomial p: Rdxd + R. M M det( 2 (M M M) is non-zero since there exists invertible skew-symmetric matrices. Therefore, the {M ıRd det(M M 0} forms a proper algebraic variety with zero Lebesgue measure, Hence, since the Gaussian distribution is absolutely continuous w.r.t. to the Lebesgue measure, it holds
P(det(G") = 0) = 0
(K.6)
and so G is invertible with full rank with probability 1. Thus
d s -=
G'
(K.7)
positive definite. Therefore, A4 is verified for even dimension d.
However, for any d 3, neither conditions
A2 nor condition
A3 is satisfied. Indeed, for any < 1. This is expected since the diffusion involves
E (R ( -rank x; single Brownian motion (see Section
Since the entries in M are independent standard normal Gaussian, we have V(G!))= ((V(M))) -
v(M).) = Then, (G*)
Σ C.ty Gky- Hence for i k
==1
d-1 2
B[ = - BIG
371
since G For + Glyy and Gky involve independent entries of M leading B[G Gkyl 0.
As consequence
E[S] [ (G)) d(d-1) la Rd.d.
K.2.2 LOCAL SPARSE TENSOR
Let us define tensor with only 2d = O(d) non-zero coefficients.
{
if id
-1 if 1[d] = =j=k k
igijj.Ed
(K.8)
otherwise.
with [d] stands for modulo d. It is built from subset of the canonical basis of skew-symmetric matrices, keeping only matrices with most non-zero values close to the diagonal. ensures a strong sparsity and local structure for
Gh...
The skew-symmetry assumption A is obviously fulfilled from the definition rank condition assumption A2 is not in general. Fortunately, the assumptions A3 and A4 still hold. In
K.8 However, the strict
particular, the Itō term matrix simplifies as shown by the following proposition.
We implemented this version of sparse tensor. For small dimension applications in Section M.6.2 and Section M.7.1 it has been found to provide numerical results as good as the dense tensor implementation (see Figures 29 and 43). For large dimension applications, dense tensors complicate or even prevent MSGM applications. There, we obtained satisfactory results with local sparse tensor (see Figure 49 in Section M.7.2
57

---

Published as a conference paper at ICLR 2026
Proposition K.2.2. If G is defined from equation K.8 then, for any E (R)(0)), rank
1 Moreover, we have S = Ig and the assumptions A1
A3
and A4 are verified.
Proof. For any (0))), we have
322
0
- 21
23
G(z) [G' ((,,,..,,GT.]
x]
22
(K.9)
.dd
O
-dd-1
31
To simplify notations, all the indices in this proof will be defined modulo d. For instance, 2++1 for i = d stands for T11-
For any y E Rd G(z) (E Rd) if and only if, for all & d, Es+134 2:4141 and 1 #c+1 Di. Finally,
y: = =
Ti+1 Ij
3; Φ1
y1-
(K.10)
Therefore, y Re. Reciprocally, we can verify that Ri Ker(G((:)). We conclude that (x))) = d (im(KerGG(z))) = = -
To evaluate the matrix S, we note that
G e ex+1e
(K.11)
with (e4) the canonical basis of Rdd Then,
1
G*
(K.12)
2
d
=
Σ
(exe] - Ex+1e
(K.13)
2
k=1
1
Σ
(0
ctes ext10 +0),
(K.14)
2 k=1
(Ia -Ta)
(K.15)
2
= la.
(K.16)
K.2.3 NON-LOCAL SPARSE TENSOR
We also propose another tensor with d(d-1) O(d) non-zero coefficients.
{
if - kid =j
G CKK'
-1 ifi kid =
1 i d, 1 K 1-11
(K.17)
otherwise.
where [ is the least integer greater than or equal to d-1 symmetric matrices. ensures a relative sparsity and encodes non-local structure for - k...'' x.
2 It is the canonical basis for skew-
Here, the sparse tensor G is of size d x x ddd-d11)) 1)/2 instead of dl x d d. Our theoretical framework differs slightly. The forward diffusion involves d(d 1)/2 one-dimensional Brownian motions. Consequently, the neural network, az s), approximating the scaled score, G(æ) V log Ps (x), has d(d 1)/2 coefficients. The size of the neural network parameters can increase and may complicate the training procedure. An alternative could be to work with neural network, (z s). which approximates the true score, V log (xx), having coefficients only.
This choice of tensor meets all the assumptions, including Al and A2 proofed below. However, because of the additional implementation complexity mentioned above, we postpone its numerical evaluation MSGM for future work.
58

---

Published as a conference paper at ICLR 2026
Proposition K.2.3. G is defined from equation K.17 then, for any Rd rank (ank(G(z))
1 Moreover, we have S = Ig and the assumptions A1
A2 A3
And are verified.
Proof. For any [0), have
G(x) = [G 1. Gid,
x
(K.18)
We already know that Im (x)) since G())
Now we assume that y E x and
define
| yax xy
(K.19)
Applying on æ, we get:
Qx
(y||=11|
zy x) y.
(K.20)
|
Besides, (Gk,k') )e, is the canonical basis of skew-symmetric matrices and Q is skew-symmetric so d)d-1)
there exists a E such that
Q =
Oge G
(K.21)
M=1 k=
and thus
['2
Qu =
ag
kk..k' x= G(x)a E Im())))).
(K.22)
#= k=1
We conclude that Im (G())) x and rank (G()))()) d - 1.
To evaluate the matrix S, we note that
G...'' c extrek
(K.23)
with (ex)s the canonical basis of Rdd and defining again all the indices modulo d.
s
(K.24)
2 ==
k=
(exe cx+xee])
(K.25)
k=1 k=1
dd- d
(0- exel - -etwee + +0,
(K.26)
2 ==1 k=1
1 (Ia + Ia),
1
= 2
(K.27)
Lg.
(K.28)
K.3 DISCUSSION ABOUT LOCAL AND NON LOCAL STRUCTURE
The random tensor Section J and the large sparse tensor of Section K.2.3 may be interpreted as non-local since G changes coefficients of æ which are not sorted next each other in For large dimension d, we believe that this can accelerate the convergence in comparison with local tensors, like the sparse tensor of Section K.2.2) discretized version of transport noise SPDEs. Indeed, for local dynamics the randomness may take time to spread by going from one coefficient to the next whereas in non-local dynamics the randomness spread directly in the whole state space
59

---

Published as a conference paper at ICLR 2026
at each time step. Our preliminary numerical results (not shown) seems to confirm this intuition. Moreover, the stronger theoretical properties of those non-local tensors rank condition A2 and thus exponential convergence of the distribution also tends to confirm conjecture. However, diffusion models in large dimension strongly rely on the powerful skills of convolutional neural intrinsic local structures. Accordingly,
networks (CNN), which have up to attention layers may be difficult for CNN to learn how to denoise non-local noising process. More theoretical and experimental works would be needed confirm this intuition. This of the scope of this already lengthy paper and is currently under investigation by the authors.
L NUMERICAL SCHEME
L.1 NUMERICAL INTEGRATION OF SDEs
L.1.1 STOCHASTIC RUNGE-KUTTA METHOD FOR STRATONOVICH SDES
We consider the Stratonovich stochastic differential equation (SDE):
da: fs((,, 201) dt G1,, dBB,,
(l.1))
where fs R2 + Rd is the drift, G x Rd gdxm is the diffusion term, and B, is an d-dimensional Wiener process.
The following Runge-Kutta (RK) method (Kloeden et al. 1992) approximates the solution #n+1 R xttn+1) over the interval [tns tx+1], with time step At in+1 tn and Wiener increment AB, B B
K1 fs(tn, Fn.
At + G(tn, Fn.)
AB
(L.2)
K
At
K
K₂ fs
in +
En +
At + G
In +
Un +
A ABp;
(L.3)
2
2
2
2
K₂
At
K₂
K3 fs
In +
In +
A/+ G
In 4
Un +
ABm:
(L.4)
2
2
2
2
K4- fs(tn
At, At, F.
Kg) At G(t,,
At, Fn K3))
AB
(L.5)
In+1 In
K + 2K + 2K3 + K.).
(L.6)
This method leverages the structure of Stratonovich SDEs and their differential geometry properties. It is particularly well-suited SDE equation 3.1 with skew-symmetric noise and no Stratonovich
drift.
L.2 RENORMALISATION
Both our forward SDE equation 3.1 and backward SDE equation 2.2 preserve the solution norm || However, even the above Runge Kutta discretization can break this symmetry. To enforce numerically, we normalize after each time step.
The final integration scheme summarized in Algorithm 2 Here, we highlight the differences compared to the classical RK4 in color. Note that the optional of normalization in line 10 of the Algorithm relevant only for MSGM but not for SGM.
L.2 SCHEDULING
In order to enable both a sufficient statistical convergence of the forward SDE at time s = T and a convenient time step, implemented time scheduling for both SGM and MSGM. We first recall the basic principle of scheduling in continuous time, then propose a method for MSGM, and finally discuss the theoretical consequences.
60

---

Published as a conference paper at ICLR 2026
Algorithm 2: SRK4 for conservative Stratonovich SDEs with renormalization.
Input: Integration time T, number of time step Nr. initial condition zo, drift fs. diffusion G
1: At Nr
(Time step)
2: for n = 0 to Ny 1 do
3:
A NO, Al Ia)
Wiener increment)
4:
la - nAi.
5:
K1 + fs(In, Enl
At + G(,,, Ln) A
6:
K₂ - fs + At In + K At + G(a
At
In + Ki ABn
7:
K3 ft + 2 Fn + K2 At G(tn
At
In + K2
A AB,,
2
8:
K1 + ft Atn Un + K3) At + G(tn
An In
K3) AB,
9:
In+1 f In + B (K1 + 2K2 + 2K3 + K4)
(Classical RK4 blend]
10:
In+1 [2++1| xn+1
(Optional step: Enforce | || = | [D)
11: end for
12: return NNr
(Approximation of ær
L.2.1 USUAL SCHEDULING
In continuous time (Song et al. 2021), convenient way is to make change of variable, replacing
the time by
=(s) S gd')ds..
(L.7)
with
1 (8mm 8m (3M - 9m )
(s) - 3(s)
3(0) -
(L.8)
and BM > 3m O. Note that other schedulings are also possible (Strasman 2024 and can possibly be optimize to better adapt to the problem at hand. We first describe the hyperparameters values chosen in our numerical experiments and then explain how scheduling affects SGM and MSGM theories.
C Since we built our code from an existing one (tttps: //tp::://httttbb Cooom// fode Huang et al. ())T)), by default we choose the values provided there for SGM scheduling: 8.. = 0.1 and 3M = 20. We expect these values to be already finely tuned and we have verified that this couple of values gives indeed better results than many other choices (not shown). We believe that these default values of the SGM hyperparameters enable a fair comparison to MSGM. For some test cases, found another SGM scheduling that works better and we use it instead. All scheduling hyperparameters are provided in the tables summarizing test cases Section[M]
For small time s, the time remapping is linear: g²(s)
-+0 3.3m and 2(s) --+ fm whereas for
large time, (2) TT 3 and using the Taylor expansion around T, yielding
=(s) x() + (1) - T) o(s -T),
we find that
2(s)
3m +
2 1,30
(Bas - 2mm) )
(L.9)
2
BM + Bmn + Bxr -T))
+ o (s- - -T),
(L.10)
2
BM
Bm T.
(L.11)
61

---

Published as a conference paper at ICLR 2026
As such, SGM forward and backward SDEs become:
dæ, (8) .s + V2g(s)dB.
(L.
dx; (T t) + v2g(T - t) [0( tt, + odB.)
(L.13)
where ag (T -t, x1) approximates 22g V log pr The backward SDE can be integrated
with the Stochastic Munge-KKutta Algorithm 2 where
ft, x.) = (T - t) + v2g(T t) -t, x1))
and G(t, x.) = V2g(T - t). (L.14)
L.2.2 SCHEDULING FOR MSGM
We propose a similar scheduling for MSGM. Scheduled forward and backward SDEs write:
dæ. g(>) ((s))(z))) dB.:
(L.
dæ. () ao 1 ())d odBB),
(L.
where ag (T t, 21) approximates g(T G( VV log poggr-cl (x)) Numerically, following Algo- rithm 2 we can integrate the forward SDE with
fs x.) and G(s, w.) y G())G((z.))
(L.
and the backward SDE with
fs(t, = g( (a) -T,,, x)) and G(t, ã.) g(T -((-)))(3))
(L.18)
L.2.3 THEORETICAL RESULTS
We verify that theoretical results remain under this time scheduling. The Fokker-Planck equation is
0s Ps
V
(L.
which can be rewritten
0s D2
V1 ( V
(L.20)
1/2 = Pz(s)-
(L.21)
Besides, for 8 S' (T) for BM > Bm Taylor expansion at s/ Fx T yields
-BmT + /B2.7 2- + 4T Bmm
(L.22)
Bas 3m
2 (s' -su) (1+ s6
(1+)0/
(1)
T T
(L.23)
3M
T.
(L.24)
Therefore, from the convergence of p (already proofed) we have the convergence of v2 p2-1(e9)
The rate of convergence is still exponential:
Ips/ Pooll (2
llp ||p-- Poo 122
(L.25)
1|p2 100) - -Pao 120llE exp(-a d)= (t()))
(L.26)
IlPo Poe 11 ((())) exp (-a d)T)
(L.27)
aGG d) BM + 3m
sf Pro-am (1) )
|+
BM/2
4
62

---

Published as a conference paper at ICLR 2026
Besides the ELBO remains valid: Po( > & (x)
== E log Po (opPo(T)lee [ ]]
(L.28)
E
2 a +]||+ (a)))(, 5)) 2o ds.
The above results guaranty the applicability of a classical scheduling to MSGM. An analyze with a scheduling more adapted to MSGM has still to be done and is left to future work. For example, the work of Strasman et al. 2024) may be adapted to the multiplicative structure our SDE.
L.3 LOss EVALUATION
Following the existing code hhttps://gite
com/ CW Crom/CC-Haaagg///aeeee www.aaanngnn.lwoooeeeeooeeeeeeee
Huang
et al. (021)), p), we sample final integration time s of the forward SDEs uniformly on [F., with
with & small. According Theorem 3.4.1 we consider the following SSM loss:
CSSM (0) E., 20
E.- 21100 E., Rad(a) (s, s. g.gG., ag.. U2),
(L.29)
with
(s x, GG, ae.. v) = (0 ))I (v G))(())(()))((t 8))
(L.30)
where E is the averaged over the generated samples. Each training sample zo of batch is chosen randomly among the train set. For each of them, we sample one time 8, one solution Øs, and one Rad(d).
slicing direction
For SGM, we take G = 2 in the above expressions and following Song et al. 2021 the solution v. of the SGM scheduled forward SDE equation L.12
a. (xp(-==(())
To + /1 expp(=((s))
(L.31)
where =(s) = So (s')d's is given by equation L.9 at is chosen randomly among the train set and 2.0 (O0,, Ia).
Unfortunately, to evaluate the MSGM loss, we cannot apply the same methodology, since, for d > 2 we are not aware of an analytic expression with nor without scheduling (equation L. 1.13and equation 3.1 respectively). We integrate that SDE numerically with the stochastic Runge-Kure Kutta method with renormalization (see Section L.I.I and Section 1.2... Through this integration, we have to compute the solution Z., for many time steps 8% == kTyNT E [0, T]. Instead of sampling a random continuous time ~ ((tt,,T)) we choose random discrete time as follow
solution the MSGM forward SDE, neither
8~ sT((, T)) with T(I., T) (s = N,.
(1, NT), > 4e).
(L.32)
The numerical integration of the forward SDE implies larger computational cost compared to SGM. Therefore, as explained in Section M.3 for fair comparisons between SGM and MSGM, the number of ADAMS iterations will be smaller.
For two-dimensional test cases, we could have used the analytic example of Section to integrate the forward MSGM SDE. However, we prefer to propose and analyze an algorithm that is not tied to the dimension 2. So, we perform all our numerical experiments with the same algorithm whatever the dimension. SGM forward equation is integrated analytically, whereas the MSGM is integrated numerically.
L.4 NEURAL NETWORK ARCHITECTURE
I.4.1 SPHERICAL DECOMPOSITION AN INPUT LAYER
In line with our spherical decomposition equation 3.6 we add fixed input layer the network used in MSGM:
ag a s) āo (at/ be/||=/|1e log || s). with ||a]|e = || +
(L..3)
63

---

Published as a conference paper at ICLR 2026
The geometrical interpretation of Section H.3 also suggests that form. For SGM, if not stated otherwise, we use a default architecture:
a 8) āø( s).
(L.34)
L.4.2 NETWORK ARCHITECTURE FOR LOW- DIMENSIINOON TEST CASES (MLP)
Following the existing code
com https///gg..taaa....ooooootttaaa
Huang
et al. 20211), we parameterize the vector field äg IR RR - IR (d = d or d + 1) with 4-layer
ML. conditioned on an index R. by concatenation. Let H 128 be the hidden width. For input
a E Rd, we form ho
[ E Rd+1 and compute
Wi B
h swish W ho b1);
h2 = swish W2h + b2).
W₂ R
h3 swish Wah2 + ba).
W3 E RH
Wah3 ba.
W4 E BaxH
with (swish(=)); 2, (z)) and the logistic sigmoid. We set āg(æ, E IR No residual connections, normalization, or dropout are used. Table summarizees the hyperparameters this default architecture.
Table 1: MLP architecture hyperparameters.
Hyperparameter
Value
Input dimension
d d or d
Index dimension
Hidden width
128
Depth
3 hidden layers
Activation
Swish (æ: xo (I))
Output dimension
Output layer
Linear
Residual connections
None
Normalization Dropout
None
L.4.3 NETWORK ARCHITECTURE FOR NIGH-DIMENSIONAL TEST CASES (UNET FOR 32 x 32 VORTICITY FIELDS)
For high-dimensional experiments of Section M.7.2 we model the score field áz(x, t) using a 2D UNet operating on images x' of size H x W representing vorticity snapshots (Hf = W = 16 or 32). Some part of our algorithm was built for vectors rather than images. So depending on the portion of the algorithm, images x' E IR Hxx are reshaped into vectors Rd with d = HW or vectors are reshaped as one-channel images x' E BlxHxW
Optional spherical premodule. When enabled, apply the spherical decomposition of Sec- tion L.4.I
(x"" log [|=]]e]] = NormaaizeeLogRRaiis(())))
x
I|=
The normalized field x is passed to the UNet, while log || is embedded through small MLP and added to the temporal embedding, giving conditioning mechanism analogous to the time embedding of diffusion models.
UNet backbone. The core architecture follows the DDPM UNet of Dhariwal & Nichol (002
fully convolutional encoder-decodee with skip connections, residual blocks, and optional attention at intermediate resolutions. We one input channel and one output channel (vorticity). Let Co denote the base width. The feature width at resolution level & is Comy where Co is the base channel width and mg is the channel multiplier.
64

---

Published as a conference paper at ICLR 2026
The UNet receives (a: t) (and optionally log ||=]|e) and computes:
(æ( t) = UNeto(reshape(()) Emb(t) + Embiog (log [=]l))),
followed by flattening back to dimension d needed.
Table 2: UNet architecture hyperparameters.
Component
Setting
Input/ / output channels
1
Input resolution
H W E (16,32)
Base channels width Co
32
Channel multipliers
(1, 2, 4)
Residual blocks per stage
2
Attention resolutions
8 x and x 4 (for 16 16 input)
Activation
SiLU sinusoidal + MLP
Time embedding
Log-norm conditioning
optional MLP added to time embedding
Dropout
Upsampling downsampling
convolutional
Output
1-channel vorticity field
This UNet is used as a drop-in replacement for the small-dimensional MLP of Section L.4.2 enabling
MSGM/SGM to scale to image-like vorticity fields to dimension d = 1024.
M DETAILS ABOUT OUR NUMERICAL EXPERIMENTS
We will show that - for comparable training time - MSGM can generate distribution of better quality than SGM when data distribution tails are heavy or close being heavy. For distributions with lighter tails such as Gaussian ones, SGM and MSGM produce similar results, except for small number of backward time steps where SGM can become unstable. MSGM more robust this aspect.
Our code can be found here: https: /tle:///giii gighuit.hhhgggo www.w//wwww..////eeeooooo sde www.://wwwwdddwwweooooooooooo Lighttttree 3rdSub and the preprocessed vorticity data we used Section M.7 can found here: ht g vitubb.. cwm/ ccom///resseeeuu/sscVweo- are freely available at
Original experimental data (Georgeault & Heitz, 2026) wttp:////......//1//////0.
through the multidisciplinary
repository Recherche Data Gouv hattps::/eetrreett..cee
data gouv. fr).
M.1 TEST CASES
We will illustrate MSGM and compare it to SGM through different test cases. We first consider four examples sampled from known distributions: the Swiss roll, a multidimensional Gaussian distribution, and the multidimensional Cauchy distribution with and without correlations. Then, we will address the experimental fluid dynamics data. For each test case, a table summarizes the nominal parameters and 7). All are performed on CPU. In addition, we additional cover a high-dimensional application with imagine processing, see section Section M.7 with GPU A40 NVL. with 48 Go of VRAM.
used in the experiments (see tables 4 5 6
M.2 DATA PREPROCESSING
The data set and distribution are centered before processing. For SGM, data sets are renormalized, component by component, by their estimated standard deviations. This preconditioning cantly reduce the number of conditioning of the covariance of the data set, and therefore facilitate the SGM (Guth al. 20022.. Generated data are then re-scaled for plots and other post- post-processings. For MSGM, it is not necessary and may even be counterproductive for conservative dynamical systems. In fact. it changes the definition of energy || | 12. The modified energy has physical meaning. It may have a very different distribution, possibly much less relevant for the data structure. So, we do not renormalize the data set before training MSGM.
signifi-
65

---

Published as a conference paper at ICLR 2026
Table Swiss roll test case: parameters of the nominal numerical experiments.
Parameter
SGM
MSGM
Dimension
2
Number of used training data points (M)
20
256
220
x 256
Number of test data points
104
104
Reference number of ADAMS steps
220
220
Number of ADAMS steps (Niter)
220
220
CPU time ADAMS steps (in ms)
4
3
Batch size
256
256
Number of time steps (forward) Ni
16
Number of time steps (backward) N7.
16
16
8min Bmax ts
0.1
0.1
20
20
10
10 10-3
Learning rate
10
10 10-3
Neural network architecture
default
spherical (equation L.33
MMD (MMD(train)= 0.9 x 10
1.9 x 10
0.9 x 10-2 2
M.3 COMPARISON STRATEGY
We will perform different qualitative visual comparisons with pairplots and quantitative assessment with Maximum Mean Discrepancy (MMD) (Gretton et al. 2012). Given two ensembles X (x())
(TR) & and V ()))) (R)V samples of random variables X and Y respectively, define MMD(r, y) as:
MMD (X
((x(0),2(01)) f-=1
Σ
-2k(x(0), ,1(0) + k(y y())),
(M.1)
N2
k(u,z) = exp(-1/ux 21)2).
(M.2)
If Xist is the test set and geen generated ensemble, MMD(X X is metric of the precision of our generated ensemble and hence our Al generative algorithm. A small MMD means close distribu- tions. However, MMD is relative metric. So we compare MMD Xtest | SGM MMD (Ttest MSGGM)
and MMD(X" (Xtest X where X and MMGM generated from SGM and MSGM respectively, and Ximin is the train set. MMD(X Ximin) provides a reference MMD, encoding in particular possible distribution shifts between the train and the test sets.
The numerical integration of the MSGM forward SDE is additional significant computational cost compared to SGM, and hence a slower training procedure. This cost scales linearly Ny due to the "for" loop in time. Empirically, appears to scale as ( = vd N2/24 (not shown), probably due to the vectorized d x d x d tensor products involved in each integration time step. In most of the numerical experiments below, N, 2 and thus Vd. The SGM iteration steps ( times faster than the MSGM iteration steps. Consequently, the number of iterations for the SGM is max(1, (SJ) times larger than the number of iterations for the MSGM. As such, we can compare the results of SGM and MSGM at similar training cost. By convention, we take the number of iterations for SGM as a reference and refer ? provide the values for the reference number of iterations, the true number of iterations, and the execution time per ADAMS step.
the reference number of iterations. Summary tables
and
M.4 SWISS ROLL
We first illustrate our method with the Swiss roll distribution. It is simple two-dimensional distribution:
py .ecaaageseesss...eeee ecs. vuw.as htt.ccc..///aaasssstt///oeee///
Its curved shape makes it difficult to grasp by linear Gaussian approaches. Both MSGM and SGM mimic the Swiss roll distribution well, as illustrated by the pairplot g However, the diffusion distribution ps differs from Figure 1 Figure 10 In particular, latent distributions are completely
66

---

Published as a conference paper at ICLR 2026
test gen.
test
gen.
4
2
4
4
-3 3 -4 4
-3
3
-4
4
Figure 8: Pair plots of generated data (orange dots) compared ground truth data (blue dots) with the SGM (left) and MSGM (right) for Swissroll data. On the diagonal, log-histogram ground truth data (continuous blue line) and logarithm of the pdf KDE estimation of generated data (orange line) are superimposed.
[+
[+2
I+
I*9
1+5
12
1*:4
1-
Figure 9: Evolution of the solution log-pdf log(p, of SGM forward SDE (top) and backward SDE (bottom) for Swiss roll data.
different. Figure I illustrates the convergence of the SGM and MSGM approaches as a function of the reference number of ADAMS iterations and a function number of time steps for integrating the backward SDE. The precision of each sampler is quantified through MMD and the confidence intervals of MMD are estimated from the samples of 10 MMD.
M.5 ANISOTROPIC GAUSSIAN DISTRIBUTION
For a complete numerical analysis, we compare SGM and MSGM on correlated Gaussian data 20 (O, AA"), with fixed matrix, A, initialized with i.i.d. coefficients A,,,~((( 1). For 32 time steps backward, the pairplots in Figures 12 13, and [1epeeenttt similar generative skills, but for or 16 time steps backward, only MSGM gives good results. For time steps backwards, MSGM still provides good distribution, whereas the SGM backward SDE completely diverges. Figures 15
Figure 10: Evolution of the solution log-pdf log(ps (21, (2), of MSGM forward SDE (top) and backward SDE (bottom) for Swiss roll data.
67

---

Published as a conference paper at ICLR 2026
10-1
10-1
SGM
MMD
MSGM
train data
10"²
10-2
28
212
216
20
2
2
2
2
2
effective number iterations
nb timesteps backward SDE
Figure 11: Convergence of MMD (mean and 80% confidence interval) for Swiss roll distribution as function of reference number of ADAMS iterations (left) and as function of number of time steps for integrating the backward SDE (right).
Table 4: Gaussian test case: parameters of the nominal numerical experiments.
Parameter
SGM
MSGM
Dimension
16
Number of used training data points (M)
1048576 x 256
262144 x 256
Number of test data points
10
10a
Reference number of ADAMS steps
220
220
Number of ADAMS steps (Niter)
2
1048576
262144
CPU time ADAMS steps (in ms)
23
Batch size
256
256
Number of time steps (forward) N/.
16
Number of time steps (backward) N₁
16
16
8min Siax t₂
0.1
0.1
20
20
10
10 10-3
Learning rate
10 33
10 3
Neural network architecture
default
spherical (equation L.33
MMD (MMD(train)=1.5 10-2)
10-2
2.5 10-2
16 and 17 also highlight this and show that the converged dynamics of the pdf ps differs between SGM and MSGM. Figure also confirms that MSGM converges faster with the number of time steps, and is generally more stable.
M.6 MULTIVARIATE CAUCHY DISTRIBUTION
Cauchy distributions are worst-case heavy-tail distributions in the sense that they do not have finite moments. Still, they appear in applications of hydrology, e.g. annual maximum one-day rainfalls and river discharges. Consequently, we analyze the expressivity of MSGM this extreme case. Note that due to the absence of finite moments, convergence in common metrics such Wasserstein-p total variation not well defined.
M.6.1 VECTOR OF INDEPENDENT CAUCHY VARIABLES
We first illustrate our method with a vector of independent Cauchy variables: &p = ICa with aCa defined by equation 6.2 with scale parameter confirm that SGM does not reproduce fat tails unlike MSGM. Moreover, SGM misaligns the far data points that have the coordinate I1 < -3 An explanation of the superior skills is the similarity between the data distribution and the latent distribution in MSGM: property not shared by SGM, as illustrated in Figures 21 and 22
1/50. As expected, Figure 19 and Figure 20
68

---

Published as a conference paper at ICLR 2026
test gen.
gen.
15
15
10
10
10
10
20
20
20
20
15
-10
10
-20
20
-10
10
-15
-10
Figure 12: Pair plots of generated data (orange dots) compared to ground truth data (blue dots) with the SGM (left) and MSGM (right) with 8 time steps backward for vector of 4 correlated Gaussian variables, among 16 correlated Gaussian variables used for training. On the diagonal, log-histogram of ground truth data (continuous blue line), and logarithm of the pdf KDE estimation of generated data (orange line) are superimposed.
test gen.
gen.
15
15
10
10
10
10
20
20
20
20
-10
15
-10
10
-20
20
-10
10
155
Figure 13: Pair plots of generated data (orange dots) compared to ground truth data (blue dots) with the SGM (left) and MSGM (right) with 16 time steps backward for vector of 4 correlated Gaussian variables, among 16 correlated Gaussian variables used for training. On the diagonal, log-histogram of ground truth data (continuous blue line), and logarithm of the pdf KDE estimation of generated data (orange line) are superimposed.
69

---

Published as a conference paper at ICLR 2026
test gen.
test gen.
15
15
15
15
10
10
10
10
20
20
1
-10
10
-15
15
-10
10
-20
20
-10
10
-15
15
10
10
-20
20
Figure 14: Pair plots of generated data (orange dots) compared to ground truth data (blue dots) with the SGM (left) and MSGM (right) with 32 time steps backward for a vector of 4 correlated Gaussian variables, among 16 correlated Gaussian variables used for training. On the diagonal, log-histogram of ground truth data (continuous blue line), and logarithm of the pdf KDE estimation of generated data (orange line) are superimposed.
F-6
2-3
F-4
1-5
F-6
/-:0
F-12
F-:4
16
Figure 15: Evolution of the solution log-pdf log(p. (x .22) of SGM forward SDE (top) and backward SDE (bottom, with time steps) for Gaussian data.
Figure 16: Evolution the solution log-pdf (og( x22) of SGM forward SDE (top) and backward SDE (bottom, with 32 time steps) for Gaussian data.
70

---

Published as a conference paper at ICLR 2026
Figure 17: Evolution of the solution log-pdf log(ps x22) of MSGM forward SDE (top) and backward SDE (bottom, with 16 time steps) for Gaussian data.
10
SGM
MSGM
train data
10-2
25
nb timesteps in backward SDE
10-1
10-1
MMD
MMD
10
10-2
2
212
236
20
28
212
216
200
effective number of iterations
effective number iterations
6x-2-
10-1
4x 2
MMD
3x-2
1x
10-2
10-2
4
16
dimension
dimension
Figure 18: Convergence of MMD (mean and 80% confidence interval) for the Gaussian data as
function of number of time steps for integrating the backward SDE Nš. (top), as a function of reference number of ADAMS iterations (middle) for Nà 8 (left) and N 32 (right), and function of dimension (bottom) for N₁. 8 (left) and Nē. 32 (right).
71

---

Published as a conference paper at ICLR 2026
test gen.
gen.
10
10
5
-5
10
10
-10
-5
5
10
-10
-5
5 10
Figure 19: Pair plots of generated data (orange dots) compared to ground truth data (blue dots) with the SGM (left) and MSGM (right) for vector of two independent Cauchy variables. On the diagonal, log-histogram of ground truth data (continuous blue line), theoretical log-pdf (dashed blue line), and logarithm of the pdf KDE estimation of generated data (orange line) superimposed.
10
10
XX|X > R)
S(R) PX|
10
10
S(R) II
2
lest
lest
gen.
gen.
10
10
10"
10 R
10'
10"
10
10'
Figure 20: Survival function of generated data (orange line) compared to ground truth data (blue line) with the SGM (left) and MSGM (right) for a vector of two independent Cauchy variables.
1+3
14
1-24
1+32
(+4)
(+45
1*35
Figure 21: Evolution of the solution log-pdf log(ps (2), of MSGM forward SDE (top) and backward SDE (bottom) for a vector two independent Cauchy variables, with fast scheduling: Bm, 0.1, 8M 0.4 and our neural network architecture based on spherical decomposition
equation L.33
72

---

Number of time steps (forward) N/. Published as a conference paper at ICLR 2026
Table Vector of independent Cauchy variables: parameters the nominal numerical experiments.
Parameter
SGM
MSGM
Dimension d
Number of used training data points (M)
20
256
209715 x 256
Number of test data points
10
109
Reference number of ADAMS steps
220
2
Number of ADAMS steps (Niter)
220
1048576
209715
CPU time / ADAMS steps (in ms)
27
Batch size
256
256
64
Number of time steps (backward) N4
128
128
Buin Binax
0.1
0.1
20
0.4
tE
10 10-3
10
Learning rate
10
10 10-3
Neural network architecture
default
spherical (equation L.33
MMD (MM. x 10 10-3)
7.5 x 10
3.3 10
[+:0
Figure 22: Evolution of the solution log-pdf (pg x2)) of SGM forward SDE (top) and backward SDE (bottom) for vector of independent Cauchy variables, with default scheduling: f3 = 0.1, 3M 20 and default neural network architecture.
Figure 19 compares MSGM with fast scheduling (3n 0.1, 8M 0.4) and neural network architecture based on spherical decomposition equation L.33 with SGM with default scheduling (8m. 0.1, BM 20) and default neural network architecture. For a fair comparison of MSGM, we complement our numerical analysis with Figures 23 we test SGM with fast and default schedulings, and with both cithercall-ecomppsitte
based and default network architectures. This fast scheduling seems not adapted to SGM, making the sample generation highly inaccurate in the pairplot of Figure[24 contrast, the network architecture with spherical decomposition does improve the SGM sampling procedure, especially for distribution tails. However, even with this architecture, SGM remains less efficient than MSGM. First, the estimated tail is less clean. Secondly, the samples generated far are not properly aligned with the test samples, especially for x2 -3. Third, outside the 21 and Ø2 axes, SGM generates too few samples close to the origin (say points æ with ||>:|| > and ||0:|| < 2).
For Cauchy distributions, we still compare MMD values. However, is well defined mathe- matically and is hardly relevant numerically. Indeed, the Gaussian kernel structure of the MMD is probably not adapted to samples that are so far from each other. In our experiments, we used 10a samples to compute an approximate MMD. Other quantities of interest can also be utilized, such as the survival function t + |[ > R), illustrated in Figure 20 As expected MSGM clearly outperforms SGM this metric. Indeed by construction
learning method is robust in terms of the radial distribution ||>||| obtained directly from the data and not after the noising process. This is valid since the norm distribution does not change in time due to equation 3.7
73

---

Published as a conference paper at ICLR 2026
test gen.
10
5
-5
10
-10
-5
5
10
-10
-5
5
10
test
test gen.
gen.
10
10
5
-5
10
10
-10
-5
5
10
-10
-5
5
10
-10
-
10 -10
-
10
Figure 23: Generated data (orange lines and dots) compared to ground truth data (blue lines and dots) with the MSGM (top left corner) and the SGM (bottom) for two-dimensional Cauchy distribution. SGM plots correspond to default scheduling: 3m, 0.1, Aa = 20. Left plots correspond to our neural network architecture based on spherical decomposition equation [.33 whereas the right plot correspond to default neural network architecture. On the diagonal, log histoiroggam of ground truth data (continuous blue line), theoretical log-pá agg (dashed blue line), and logarithm of the pdf KDE estimation of generated data (orange line) are superimposed.
74

---

Published as a conference paper at ICLR 2026
test gen.
gen.
10
5
-5
DW
1
10
10
-10
10
-10
10
Figure 24: Generated data (orange lines and dots) compared to ground truth data (blue lines and dots) with the SGM for two-dimensional Cauchy distribution. Plots correspond to a fast scheduling: 3. = 0.1, 3M 0.4. The left plot corresponds to our neural network architecture based on spherical decomposition equation L.33 whereas right plot corresponds to default neural network architecture. On the diagonal, log- histogram ground truth data (continuous blue line), theoretical log- t-f (dashed blue line), and logarithm of the pdf KDE estimation of generated data (orange line) superimposed.
1-14
J-24
1-30
F-4
[-
Figure 25: Evolution of the solution log- oog-pdf log(ps .2)) of SGM forward SDE (top) and backward
SDE (bottom) for vector of two independent Cauchy variables, with default scheduling: 3n 0.1,
BM = 20 and our neural network architecture based on spherical decomposition equation L.33
Figure 26: Evolution of the solution log-pdf log(ps (T x2)) of SGM forward SDE (top) and backward SDE (bottom) for vector of two independent Cauchy variables, with fast scheduling: 3. 0.1. 9M = 0.4 and default neural network architecture.
75

---

Published as a conference paper at ICLR 2026
Figure 27: Evolution of the solution log-pdf log (P x2)) of SGM forward SDE (top) and backward SDE (bottom) for vector of two independent Cauchy variables, with fast scheduling: f., 0.1, AM = 0.4 and our neural network architecture based on spherical decomposition equation L.33
Table 6: Vector of correlated Cauchy variables: parameters of the nominal numerical experiments.
Parameter
SGM
MSGM
Dimension
Number of used training data points (M)
220
256
220
256
Number of test data points
10a
103
Reference number of ADAMS steps
220
220
Number ADAMS steps (Niter)
220
220
CPU time / ADAMS steps (in ms)
45
Batch size
256
256
Number time steps (forward) N/
128
Number time steps (backward) No.
128
128
8min Bmax
0.1
0.01
t
10 10-4
10 10-4
Learning rate
10
10 10-3
Neural network architecture
default
spherical (equation L.33
MMD 3M3 x 10-3)
11.2 x 10-3
5.2 10
M.6.2
VECTOR OF CORRELATED CAUCHY VARIABLES
To address dimensionality issues, we consider the correlated Cauchy variables already presented in Section 6.I. In terms of survival function, MSGM expected more accurate than SGM (see
Figure 28
10
10
S(R) = xX| > R)
S(R) n [[| > R)
10
10
10
10
test
gen
gen.
10
10
10" -2
10
10'
10
10"
10"
10
10'
10
R
R
Figure 28: Survival function of generated data (orange line) compared to ground truth data (blue line) with the SGM (left) and MSGM (right) for vector of correlated Cauchy variables.
76

---

Published as a conference paper at ICLR 2026
test gen.
test gen.
10
10
10
10
10
10
10
10
10
10
-10
10
-10
10
-10
10
-10
10
-10
10
-10
10
-10
10
-10
10
test
gen.
gen.
10
10
10
10
10
10
10
10
10
-10
10
-10
10
-10
10
Figure 29: Generated data (orange lines and dots) compared to ground truth data (blue lines and dots) with the MSGM (top) with dense tensor G (top left corner), with sparse local tensor G (top right corner), and the SGM (bottom) for a vector of 4 correlated Cauchy variables. SGM correspond default scheduling: 8m, 0.1, A3 20. Left and top plots correspond network architecture based on spherical decomposition equation L.33 whereas the right bottom plot corresponds to default neural network architecture. On the diagonal, log-histogram of ground truth data (continuous blue line), and logarithm of the pdf KDE estimation of generated data (orange line) are superimposed.
neural
As for independent Cauchy variables, we present complementary numerical experiments with different scheduling and different neural network architectures Figures 29 and 30 Figures 35 unveil the corresponding diffusion dynamics from = 0 to ' and from / 0 to / T. Again, the neural network architecture based on spherical decomposition significantly improves the SGM generative skills but MSGM remains more efficient sampler. Not all the branches of the star-like pdf are well sampled and, outside the branches, the regions near the origin is not well sampled.
One wonder if the poorer results of SGM would improve for larger number of ADAMS iterations. To answer this question, we run longer experiments with 224 = 16777216 ADAMS iterations. Figures 36 37| and 4a show that MSGM slightly improves with increasing number of iterations, whereas SGM diverges. For a fair comparison, the MMD convergence Figure 4a is expressed in terms of effective number of ADAMS iterations, i.e. proportionally reduce
77

---

Published as a conference paper at ICLR 2026
test gen.
gen.
10
10
10
10
10
10
10
10
10
-10
10
-10
10
-10
10
Figure 30: Generated data (orange lines and dots) compared to ground truth data (blue lines and dots) with the SGM for vector of 4 correlated Cauchy variables with fast scheduling: 8m 0.01, 3M The left plot corresponds to our neural network architecture based on spherical decomposition equation 33 whereas the right plot corresponds to default neural network architecture. On the diagonal, log- histogram of ground truth data (continuous blue line), and logarithm of the pdf KDE estimation of generated data (orange line) are superimposed.
1-B
1-16
1-30
2-0
7-44
--80
Figure 31: Evolution of the solution log-pdf log(p,( #3)) of MSGM forward SDE (top) and backward SDE (bottom) for vector of 4 correlated Cauchy variables, with fast scheduling: 3. 0.01, BM and our neural network architecture based on spherical decomposition equation L.33
Figure 32: Evolution of the solution log-pdf log(ps (r)) of SGM forward SDE (top) and backward SDE (bottom) for of 4 correlated Cauchy variables, with default scheduling: B. 0.1. 9M = 20 and default neural network architecture.
78

---

Published as a conference paper at ICLR 2026
Figure 33: Evolution of the solution log-pdf ( of SGM forward SDE (top) and backward
SDE (bottom) for vector of 4 correlated Cauchy variables, with default scheduling: Am 0.1,
3M = 20 and our neural network architecture based on spherical decomposition equation L.33
06
1-:4
1-04
1-00 /-95
1-323
Figure 34: Evolution solution log-pdf (og( .3)) of SGM forward SDE (top) and backward SDE (bottom) for vector of 4 correlated Cauchy variables, with fast scheduling: 8... = 0.01, 8M 1 and default neural network architecture.
1-:4
J-41
1-112
1-
Figure 35: Evolution of the solution log-pdf log (ag(p .x)) of SGM forward SDE (top) and backward
SDE (bottom) for a vector of 4 correlated Cauchy variables, with fast scheduling: 3,, 0.01, 8 = and our neural network architecture based on spherical decomposition equation L33
79

---

Published as a conference paper at ICLR 2026
10
10
Figure 36: Pair plots of generated data (orange dots) compared to ground truth data (blue dots) with the SGM (left) and MSGM (right) with 2 = 16777216 ADAMS iterations for vector of 4 correlated Cauchy variables. On the diagonal, log-histogram of ground truth data (continuous blue line), and logarithm of the pdf KDE estimation of generated data (orange line) superimposed.
F+B
1-16
1-30
2-4
7-64
C-80
F-323
Figure 37: Evolution of the solution log- log(p. Fz)) of SGM forward SDE (top) and backward
SDE (bottom, 221 16777216 ADAMS iterations and 2 512 time steps) for vector of 4 correlated Cauchy variables, with default scheduling: 3m = 0.1, BM = and default neural network architecture.
the number of ADAMS iterations for MSGM in order to make the CPU training time of SGM and MSGM similar (see Section M.3 for details). For SGM with very large number of iterations (221_ = 16777216), use larger number of time steps (2 512) for the backward SDE prevent all samples generated by SGM to diverge.
M.7 VORTICITY FIELD FROM PARTICLE IMAGE VELOCIMETRY MEASUREMENTS
Particle Image Velocimetry (PIV) an experimental technique to measure velocity fields in fluids by tracking the displacement of tracer particles between consecutive images illuminated with lasers (Adrian & Westerweel 2011). We used two-dimensional, two-component (2D2C) PIV data of Figure 38 which provide both in-plane velocity components. Here PIV is not time-resolved, i.e. each velocity image is uncorrelated to the next. The flow observed is benchmark configuration a wake flow at Reynolds number Re 3900 created by circular cylinder embedded in mean stream (Parnaudeau et al. 2008). We compute the ww--dimensional curl the velocity. Named vorticity, presented in Figure 39
80

---

Published as a conference paper at ICLR 2026
10.0
10.0
60
50
50
V(mmn)
0.0
0.
-2.
-2.
-00
-00
60
-72
400
-13
[[mm]
Figure 38: 2D2C PIV velocity field: velocity component along x (left) and velocity component along y (right).
vorticity (1as
vorticity (1/s)
[mml]
-2
-00
-10
-20
40
60
[mml]
Figure 39: The full two-dimensional vorticity (left) and a zoom (right) of PIV field
81

---

Number of time steps (forward) N/. Published as a conference paper at ICLR 2026
Table Low-dimensional vorticity test case: parameters of the nominal numerical experiments.
Parameter
SGM
MSGM
Dimension d
16
16
Number of used training data points (M)
2
1024
1024
Number of test data points
6476
6476
Reference number of ADAMS steps
22
2
Number of ADAMS steps (Niter)
2
1048576
262144
CPU time ADAMS steps (in ms)
32
Batch size
256
256
16
Number of time steps (backward) N1
8
9min Pinxx
0.025
0.025
te
10
2.5 x 10-5
Learning rate
10 3
10 10-3
Neural network architecture
default
spherical (equation L.33
MMD (MMD(train)=0.9>> x 10-3)
1.5 10-2
1.3 10-2
vorticity
5.5
15.
50
-11.0
-10.0
k[mm]
Figure 40: Spatial cropping and spatial subsampling of vorticity field to obtain data sample at low dimension 16.
M.7.1 LOW-DIMENSIONAL TEST CASE: VORTICITY EVALUATED SEVERAL SPATIAL POINTS
To reduce the dimension of the data, we severely crop the vorticity images and subsample them spatially, keeping only 4 x 4 pixels by images as illustrated by Figure 40 Once reshaped as a vector, each small image represents data point of dimension 16. we choose a dimension & 16, we just keep the first d coefficients of the vector. For this experimental dataset, we investigate the influence of the amount data available for learning. Our default experiments will train the models with 210 1024 data points only.
As seen previously in Section 6.2 MSGMM is more robust in low-data mode and better represents rare events, also confirmed by the survival function Figure 41 We explain it by latent distribution
close to the data distribution as illustrated in Figure 42
For fair numerical comparison, also test SGM with and without our neural network architecture
based on spherical decomposition equation L.33 in Figures 43 44 and 45 This architecture improves the quality of the generated samples. However, tails are still underestimated and some regions of the space remain clearly badly sampled. In contrast, MSGM samples well the data distribution both with dense and with sparse tensor, G.
82

---

Published as a conference paper at ICLR 2026
10
10
S(R) P[] >> R)
S(R) = P[]| > R)
10
10
-2
10
test
test
gen.
10
10
10
10'
10"
10
Figure 41: Survival function of generated data (orange line) compared ground truth data (blue line) with the SGM (left) and MSGM (right) trained on 1024 16-dimensional data points representing PIV-based vorticity fields.
Figure 42: Evolution of the solution log-pdf ( of MSGM forward SDE (top) and backward SDE (bottom) for the vorticity images distribution, with nominal scheduling (8,, 0.025. AM = 5.0) and our neural network architecture based on spherical decomposition equation L.33
To complete the numerical study, we evaluated the MMD between generated samples and test samples for different values of the reference number of ADAMS iterations. different number of time steps to integrate the backward SDE, different dimension d, and different numbers of training data. The convergence plots are visible rare events. We can observe some tendencies, but definite conclusions may not be obtained from those convergence plots. For a very small training set (26 = 64 data points), both SGM and MSGM fail and MMDs are similarly large. The biggest MMD gap between SGM and MSGM appears to be the intermediate region: 2 = 1024 data points. As expected, this gap seems to increase with dimension, even though this tendency is not fully clear for the plot. For small numbers of ADAMS iterations or small numbers time steps, MSGM seems much better than SGM. This is expected since the MSGM latent space is already close the data distribution. Without enough ADAMS iterations, neither the MSGM nor the SGM samples accurately mimic the data distribution, and in any case, it better to let the optimization procedure for long enough time.
[igures/46landpp Again, MMD may not be the best tool for studying
M.7.2 HIGH-DIMENSIONAL TEST CASE: VORTICITY IMAGE PROCESSING
To demonstrate that MSGM can address high-dimensional problems, we propose here image generator based on the sparse local tensor of Equation (K.8) and the Unet detailed in Section L.4.3 From the original high-resolution PIV-based vorticity images of Figure 39 crop, subsample at resolution 64 64, smooth and subsample again images them spatially, keeping 32 x 32 pixels by
images as illustrated by Figure 48 Once reshaped vector, each small image represents data point of dimension 1024.
Figures 49 and 50 present generated images with MSGM and SGM respectively. Table summmrizes the parameters of our numerical experiment. The numerical evaluation of image generation skills of MSGM is beyond the scope of this paper and postpone this study to future work.
83

---

Published as a conference paper at ICLR 2026
test gen.
test gen.
3
3
4
3
3
-4
33
-4
3
3
test
test gen.
gen.
3
4
3
3
3
3
-4
-3
4
-3
Figure 43: Pair plots of generated data (orange dots) compared to ground truth data (blue dots) with the MSGM with dense tensor G (top left), sparse local tensor G (top right), and SGM (bottom) trained on 1024 16-dimensional data points representing PIV-based vorticity fields. Left and top plots correspond our neural network architecture based on spherical decomposition equation L.33 whereas the right bottom plot correspond to default neural network architecture. On the diagonal log-histogran of ground truth data (blue line) and logarithm of the pdf KDE estimation of generated data (orange line) are superimposed.
2-3
F-4
1-5
F-0
F-:4
Figure 44: Evolution of the solution log-pdf log (x) xx)) of SGM forward SDE (top) and backward SDE (bottom) for the vorticity images distribution, with nominal scheduling (S., 0.025, 3a 5.0) and default neural network architecture.
84

---

Published as a conference paper at ICLR 2026
/+33
/+32
Figure 45: Evolution of the solution log-pdf (og(p.( xz)) of SGM forward SDE (top) and backward SDE (bottom) for the vorticity images distribution, with nominal scheduling (3.. 0.025, ,a 5.0) and our neural network architecture based on spherical decomposition equation L.33
10-1
SGM
MMD
MMD
MSGM
train data
10-2
10-2
26
210
26
26
210
216
Training set size
Training set size
10-1
10-1
MMD
MMD
10-2
10
26
26
210
266
Training set size
Training set size
Figure 46: Convergence MMD (mean and 80% confidence interval) for the vorticity images distribution function of number of training data for (from left to right and from top to bottom) dimension 2, 4, 8, and 16.
85

---

Published as a conference paper at ICLR 2026
2x10-
SGM
MMD
MSGM
10
train data
10-2
2
2
26
2
23
effective number of iterations
nb timesteps in backward SDE
3 10
x 10
MMO
10"²
dimension
Figure 47: Convergence of MMD (mean and 80% confidence interval) for the vorticity images distribution as a function of reference number of ADAMS iterations (top left), as a function of number of time steps for integrating the backward SDE (top right), and as function of dimension (bottom).
-55
Figure 48: Spatial cropping and subsampling (left), spatial smoothing (middle), and spatial sub- sampling again (right) of vorticity field to obtain a data sample at lower but still high dimension 1024.
86

---

Published as a conference paper at ICLR 2026
24
30
Neriody
Figure 49: 32 32 image generation from MSGM with forward (top) and backward diffusion (bottom) at time (from left to right) T t = 0, 0.25, and T = 1.
wotk 25N0
Auruuyy CTN
keruuty Cle)
Figure 32 x 32 image generation from SGM with forward (top) and backward diffusion (bottom) at time (from left to right) T - / 0, 0.25, and 7 = 1. The apparent heteroskedasticity in the diffusion due to the data normalization (pixel-wise variance larger on top and bottom boundaries).
87

---

Published as a conference paper at ICLR 2026
Table 8: High-dimensional vorticity test case: parameters of the nominal numerical experiments.
Parameter
SGM
MSGM
Dimension
1024
1024
Number of used training data points (M)
5000
5000
Number of test data points
2500
2500
Reference number of ADAMS steps
10a
10a
Number of ADAMS steps (Niter GPU time ADAMS steps (in ms)
10
10a
410
590
Batch size
128
128
Number of time steps (forward) N/
128
Number of time steps (backward) N
2048
2048
Mmin Smax ta
0.8
0.8
160
160
8 x 10-3
8 x 10
Learning rate
10-4
10-4
Neural network architecture
default
spherical (equation L.33
MMD (MMD(train)==1.44 10-3)
2.4 10-3
3.2 10-3
N SUMMARIZED COMPARISON OF MSGM AND SGM
This section is devoted to a brief comparison of these two concepts of generative modeling both from theoretical and empirical point of views.
Each strategy follows its noising process, leading to different invariant distributions, i.c. Gaussian for SGM and rotational invariant for MSGM. Both latent spaces are tractable, allowing for fast initial sample generation for the reverse process. As particular added on, the latent distribution of MSGM allows for finite KL. divergence when compared to heavy-tail distribution, e.g., as discussed and motivated by Section E.6 From the convergence speed, both dynamics allow for exponential convergence to the invariant distribution, assuming the rank condition A2 is satisfied for G. We will conclude this section with comparison discussion beyond the heavy tail case.
N.1 THEORETICAL ASPECTS
The latent space of MSGM is data aware, which ensures smaller KL-divergence of target distribution and latent distribution compared to classical SGM, see Section and Proposition E.3.1 The method allows for inductive bias based on physics in the design of G. For example, in the context of transport noise, making the noising/denoising process more physically relevant. This topic is part of future work by the authors and is briefly discussed in Section 7 Moreover, the conservation of norm in the denoising/backwarr process of MSGM serves stabilization tool, both for training and for sampling stage. In particular, samples cannot diverge.
At first glance, MSGM offers drawbacks compared to SGM. First, we have rely SSM and cannot apply DSM since we do not have access to an analytic score solution of the noising process. Second, we have to rely on numerical integration in the training because of no available analytic solution; see also the empirical discussion N.2 below.
When it comes to scalability, as d oc, the current theoretical analysis is not yet complete. The current analysis is built on the (strong) rank-condition which can be verified in the case of dense tensors; see Section]] This is limit terms of scalability due the d scaling of G. Here, the sparse tensors discussed in Section K will serve as solution when it comes to scalability. However, in this context the rank condition has to be relaxed and new analysis is required as outlined in Section K.I.I
N.2 EMPIRICAL ASPECTS
SGM offers exact integration of the noising process, while MSGM relies on numerical integration. Although this first glance looks like drawback in praxis, for most of our test cases, only a few
88

---

Published as a conference paper at ICLR 2026
forward steps were needed in the training process, making the training traceable and comparable to SGM training based on exact integration, while offering the same quality. For a more detailed discussion, refer the fair comparison discussion in Section M.3 As current experiments suggest, MSGM requires less data in training. From approximation theory, learning the score reduces to training a support that is the hyper-sphere in TRd with conditioning variable log || IR.
In particular, the effective domain for learning neuronal remains bounded in d. It may affect the stability of the approximation using such an approximation class, Finally, the stabilization due the conservation of norm avoids divergence instabilities of SSM solvers for MSGM, when compared to well known instabilities of SSM solvers for SGM.
89


Reviews:
[{"review_id": "GADaNVAmqt", "reviewer": "Reviewer_Bt3f", "summary": "This paper introduces a new framework for diffusion models based on the backward prcess corresponding to a forward process which incorporates a geometric brownian motion in its SDE. This geometric term is the result of multiplying a matrix to dB_t. They further assume this term is given by the product of a vector to a three dimensional tensor for which all the sliced matrices are skew symmetric. This skew symmetric property leads to a nice norm-invariant property for the forward and backward process. The authors claim that this norm invariance property enables the diffusion model to sample from rare events in the case where the target distribution is heavy tailed. In this case, the norm invariance implies the latent distributions are also heavy tailed, so one does not end of with a gaussian variable at the end. On the other hand, starting from a gaussian, one cannot transform a gaussian random variable back to the target using the backward process; e.g. the authors show the kl divergence of gaussian with a heavy tailed distribution such as Cauchy is infinite. Hence, this diffusion could cover cases in which the data distribution is heavy tailed.", "questions": "1) is there a benefit of in adding a bias term to the forward process (eg f(x) = -x similar to the OU process) in the multiplicative framework?\n\n2) is there an argument one can make about the effectiveness of your method, either theoretically or empirically, compare to traditional diffusions beyond the heavy tail case?\n\n3) it would be nice to show the derivation of the loss and the backward sde for the new multiplicative diffusion.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The authors propose a new framework for diffusions which can be useful for heavy tailed targets, and cover more natural data distributions in some cases. The study the fokker plank equation for this process and drive the properties of the SDE such as norm-invariance.\n\nThe writing and proofs are written in a clear way and without typos, unfortunately I didn't have time to check the details of most of the proofs.", "weaknesses": "Beyond the heavy tail case, they don't really study the effectiveness of their setup compared to traditional diffusion and normalizing flows, hence when would be beneficial to use their framework. their experiments in Figure 3 is vague and I don't understand the catch from the plot on the effectivess of their method compared to gassian diffusion. the amount of experiments is also very limited, and they further don't study this usefulness version traditional non-multiplicative frameworks theoretically.\n\nothers:\nline 206 possibility of rare events is not clear at this point of the paper. maybe explain more what you mean", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "368JMKv403", "reviewer": "Reviewer_QLpk", "summary": "The authors propose a \"Multiplicative Score-based Generative Model\" (MSGM), a new diffusion framework using skew-symmetric multiplicative noise instead of additive Gaussian noise. This construction preserves the $L_2$-norm of the data throughout the forward process. The resulting latent distribution is non-Gaussian; its norm distribution is identical to the data's norm distribution, and its direction is uniform on the d-sphere. The authors derive the reverse dynamics and an ELBO-equivalent score-matching objective. They demonstrate superior performance on low-dimensional, heavy-tailed data (Cauchy and fluid vorticity) compared to standard SGMs.", "questions": "* The $O(d^3)$ complexity of the dense tensor $G$ is the primary barrier to practical use. You mention sparse tensors as future work. Have you investigated what degree of sparsity is permissible while still satisfying the rank condition (A2) and ensuring the empirical convergence to a uniform spherical distribution?\n*  The numerical integration of the forward SDE (Algorithm 1, line 5) is a major bottleneck. How sensitive is the final model quality to the number of steps ($N_T$) used in this integration? Using a small $N_T$ would be faster but would also increase the error in the $\\vec{x}_s$ samples used to approximate the score matching loss.\n\n*  In the reverse process, the score for the direction is conditioned on the norm $||x_0||$. This implies the direction and magnitude are re-coupled during generation. Can you provide more intuition on how the learned score  successfully models this complex conditional distribution to reverse the \"whitening\" of the forward process?\n\n* The paper's premise is that preserving the data norm is a significant advantage over SGMs, whose latent norm converges to a $\\chi_d$ distribution. However, in the high-dimensional regime ($d \\to \\infty$), the norm of a standard Gaussian vector $\\vec{x} \\sim \\mathcal{N}(0, I_d)$ strongly concentrates around $\\sqrt{d}$. If the input data is similarly normalized (a standard practice), the norm is effectively stable in the SGM forward process, just converging to a different constant. Does the exact preservation of the norm distribution offer a fundamental advantage over the concentration of the norm, especially given the $O(d^3)$ computational cost?\n\n* Relatedly, your argument for MSGM's superiority rests on modeling heavy-tailed data. In high dimensions, heavy-tailed distributions are often characterized by anisotropic, sparse, or low-dimensional structures, not just by an extreme radial (norm) component. Are you sure that this 1D norm-preservation mechanism is the correct inductive bias for high-dimensional heavy tails, or is the advantage you observe in low-d experiments an artifact that will not scale?", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "* **Novel Process:** The paper introduces a genuinely new forward diffusion process based on conservative, multiplicative noise. This is a clear theoretical departure from the standard additive noise paradigm.\n* **Theoretical Grounding:** The model is well-supported by theory. The authors provide the Fokker-Planck equationm prove exponential convergence to the latent distribution, and formally link the Sliced Score Matching loss to an ELBO , justifying the training procedure.\n* **Tractable Non-Gaussian Latent:** The key result—that the data norm is preserved—leads to a tractable, \"data-aware\" latent space. The 1D norm distribution can be modeled simply (e.g., with eCDF), and the direction is uniform. This avoids the prior mismatch (e.g., infinite KL divergence) that plagues standard SGMs when applied to heavy-tailed data.\n* **Empirical Niche:** The model shows strong performance in a specific, known-weakness area for SGMs: heavy-tailed distributions. The results on Cauchy data, while low-dimensional, are compelling and demonstrate a clear advantage in modeling both tail decay and directionality.", "weaknesses": "* **Prohibitive Scalability:** The model, as presented, is fundamentally unscalable. The multiplicative noise operator $G$ is a third-order tensor with $O(d^3)$ coefficients. The authors' experiments are restricted to trivial dimensions (e.g., $d=4, d=16$). The authors concede this $d^3$ cost is \"prohibitive\" for any high-dimensional application, such as image processing. This limitation confines the method to being a theoretical curiosity, not a practical tool.\n* **Computational Bottleneck:** The forward process \n\n$d\\vec{x}_{s}=G(\\vec{x}_{s})\\circ d\\vec{B}_{s} $ \n\nhas no analytic solution for $d>2$. This forces the use of numerical SDE integration (SRK4) during *every single training step* to obtain $\\vec{x}_s$. This adds a massive computational burden and introduces discretization error directly into the loss calculation, a problem SGMs with analytic forward processes do not have.\n* **Training Objective:** The lack of an analytic forward score compels the use of Sliced Score Matching (SSM). The authors admit this is a \"less stable approach\" than Denoising Score Matching (DSM). This instability is visible in the experiments (Figure 32), where the standard SGM (also trained with SSM) diverges completely, suggesting the training dynamic is brittle. MSGM's success here may be less about the superiority of its process and more about its \"closer\" prior simply making the unstable SSM objective *manageable*.\n* **Arbitrary Noise Operator:** The tensor $G$ is constructed by sampling random skew-symmetric matrices. This choice is arbitrary, lacks any physical or data-driven justification, and serves only to satisfy the theoretical assumptions (A1, A2) \"almost surely\". The authors' own \"future work\" suggestion to use sparse, physics-based tensors essentially acknowledges the non-viability of the operator used in the paper.\n* **Limited Scope:** The experiments are hyper-focused on a niche (low-d, heavy-tails) where the model is *designed* to win. There is no evidence it offers any advantage—and significant disadvantages in cost and complexity—for standard, high-dimensional benchmarks where SGMs are state-of-the-art.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "XvLn8wnn93", "reviewer": "Reviewer_MUAR", "summary": "The paper proposes to replace the additive Ornstein-Uhlenbeck (OU) diffusion process used in score-based generative models (SGM) with a *multiplicative* noising process $d \\overrightarrow{\\mathbf{x}}_s = \\mathbf{G}(\\overrightarrow{\\mathbf{x}}_s) \\circ d \\overrightarrow{\\mathbf{B}}_s$ where $\\mathbf{G}: \\mathbb{R} \\to \\mathbb{R}^{d \\times d}$ is linear with $\\operatorname{rank}(\\mathbf{G}(\\mathbf{x})) = d-1$ and each slice $\\mathbf{G}^k$ is skew-symmetric, which is inspired by transport noises in fluid dynamics.\nThe authors provide the corresponding reverse SDE and probability-flow ODE and show that the forward SDE preserves the distribution of data norms, while the distribution of the directions converge exponentially fast to the uniform distribution on the sphere. Hence, the latent distribution is non-Gaussian, which may be more appropriate for heavy-tailed and anisotropic data, while still being easy to sample.\nSince an analytic solution for the conditional score is unavailable, popular denoising score matching losses cannot be used, and the authors resort to sliced score matching to train multiplicative SGMs, which is shown to be equivalent to ELBO maximization.\nEmpirical results on low-dimensional data (correlated Cauchy, fluid dynamics data) suggest that multiplicative SGMs fit the data distribution's tail more accurately than vanilla SGMs.", "questions": "+ It seems that in all experiments, $\\mathbf{G}$ is randomly constructed (Eq. 6.1) and all $\\mathbf{G}^k$ are independent. Would it make sense to \"couple\" them in a way that makes the forward process more benign? \n+ Could $\\mathbf{G}$ be learned from data?\n+ Since the network learns $a_\\theta(\\mathbf{x},t)\\approx \\mathbf{G}(\\mathbf{x})^\\top\\nabla\\log p_{T-t}(\\mathbf{x})$, it seems that the network is tightly coupled to the choice of $\\mathbf{G}$ in the forward process (which is different in regular SGMs). Under what conditions can the score be recovered from $a_\\theta$?", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "+ **Novel diffusion process in generative modeling**\n\t+ Modifying the diffusion process which underlies modern SGMs is an important line of research, as it may lead to more meaningful latent spaces, and can incorporate prior knowledge into the corruption process (e.g. physical priors in $\\mathbf{G}$). to yield data-dependent latent distributions adapt to the   and the standard forward process  \n+ **Data-dependent latent distribution**\n\t+ Typical OU diffusion processes yields a latent distribution independent of the data, while multiplicative SGMs (MSGMs) are equipped with a latent distribution which depends on the data distribution, which may be beneficial in low-data environments, or in distributions where rare events dominate.\n\t+ While the latent distribution is non-Gaussian in MSGMs, it is still easy to sample by e.g. using the 1D empirical CDF of data norms.\n+ **Theoretical results**\n\t+ The authors provide an extensive theoretical analysis of the proposed technique, which is an important contribution. Future work on multiplicative diffusion processes may benefit from these results.", "weaknesses": "+ **Model training & Scalability**\n\t+ Since the conditional score does not admit a known analytic solution (in finite time), denoising score matching (DSM)---the de-facto standard approach to train modern diffusion models---cannot be used to train MSGM, and techniques like sliced score matching must be used, which are less stable in practice.\n\t+ The forward SDE must be numerically integrated since no analytical solution is known for $d > 2$. This is computationally expensive, which hinders scalability to larger models/problems. \n\t+ A dense tensor $\\mathbf{G}$ requires $d^3$ coefficients, which prohibits the use of this general method in high-dimensional problems (e.g. images, video, physical problems).\n\t\n+ **Empirical results**\n\t+ The empirical results in the main text are purely qualitative. While there are quantitative evaluations in the appendix, I suggest that important results are moved to the main text.\n\t+ Some of the experimental results are unconvincing: For example, I cannot see any improvement over SGM in Figure 4, and Figure 7 (right) seems to suggest that SGM converges faster to the data distribution (in terms of MMD).\n\t+ Due to the computational burden, only low-dimensional experiments are conducted. However, as SGMs thrive especially in complex, high-dimensional domains, providing experiments in these settings would be a valuable addition to this work.\n\t\n+ **Experimental setup**\n\t+ The main text is missing important details regarding the experimental setup, especially w.r.t. the baseline SGM: Was it trained with DSM or also with SSM?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "CPqDCOyapR", "reviewer": "Reviewer_Cw2o", "summary": "This paper proposes multiplicative score-based generative models, replacing the usual additive Gaussian noise with skew-symmetric multiplicative noise that preserves data norms and yields conservative dynamics. The author prove properties of the forward process and show that norm is invariant, leading to a tractable, non-Gaussian latent that is data-aware. The proposed method is tested on heavy-tailed correlated Cauchy vectors and experimental fluid vorticity data.", "questions": "- How sensitive are results to the choice and structure of the tensor $G$? Could physics-motivated sparse or low-rank $G$, as discussed, preserve the theory while reducing cost?\n- While the multiplicative noise idea is interesting, the derivation of the reverse process and training objective is solid, and the experiments are clear and well-motivated, I suggest the authors include missing related literature [1,2] that has already addressed similar generalizations of diffusion models or the so-called denoising Markov models from a higher level. Especially, Section 5.1 of [2] designs a denoising diffusion model with geometrical Brownian motion, which looks very much a prototype of the multiplicative noise studied in this manuscript.\n- It seems that the authors claim that using multiplicative noise could lead to possible steady norm through the forward/backward process. I wonder if the same objective could be achieved by Riemannian manifold-related techniques, e.g., in [3].\n\n[1] Benton, Joe, et al. \"From denoising diffusions to denoising markov models.\" Journal of the Royal Statistical Society Series B: Statistical Methodology 86.2 (2024): 286-301.\n[2] Ren, Yinuo, Grant M. Rotskoff, and Lexing Ying. \"A unified approach to analysis and design of denoising markov models.\" arXiv preprint arXiv:2504.01938 (2025).\n[3] De Bortoli, Valentin, et al. \"Riemannian score-based generative modelling.\" Advances in neural information processing systems 35 (2022): 2406-2422.", "limitations": "", "rating": 4, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "- Using skew-symmetric multiplicative noise preserves data norms, yielding a principled generalization of Gaussian latents.\n- Reverse SDE/ODE is derived and the equivalence to ELBO is proved for the proposed method.\n- Experiments show meaningful improvements.", "weaknesses": "- The dense third-order tensor $G$ has $O(d^3)$ parameterse and should be prohibitive for high-dimensional applications.\n- The reported CPU time is also an order of magnitude higher for MSGM than SGM, which may limit practicality.\n- While the focus on heavy-tailed distributions is interesting, adding stronger baselines or more metrics showing the proposed method's improvement in this setting would better support the superiority claims.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)