We highlight the fact that the assumption A4 is weaker than the assumption A2 as stated by Proposition K.1.2. It is actually a strictly weaker assumption since the tensors defined in Appendix K.2.1 and Appendix K.2.2 respect assumption A4 but not assumption A2.

**Proposition K.1.2.** Let the assumption A1 holds. Then, the assumption A2 implies the assumption A4.

*Proof.* If the assumption A2 holds, then, Theorem 3.3.1 implies that $\vec{x}_s \xrightarrow{\mathcal{L}} \vec{x}_\infty = \|\vec{x}_\infty\|\vec{x}_\infty^n$. The asymptotic latent direction, $\vec{x}_\infty^n$, is independent of the initial condition $\vec{x}_0$ and has zero mean. Therefore,

$$ \mathbb{E}[\vec{x}_s | \vec{x}_0] \xrightarrow{s \to \infty} \mathbb{E}[\vec{x}_\infty | \vec{x}_0] = \mathbb{E}[\|\vec{x}_\infty\| \vec{x}_\infty^n | \vec{x}_0] = \|\vec{x}_0\| \mathbb{E}[\vec{x}_\infty^n | \vec{x}_0] = \|\vec{x}_0\| \mathbb{E}[\vec{x}_\infty^n] = 0, \quad (K.3) $$

and by Proposition K.1.1, the assumption A4 holds. □

## K.2 Sparse tensors

Here we propose several possible choices of sparse tensors.

First, we will consider a simple low-rank tensor in Appendix K.2.1 and show that it makes the latent distribution untractable. Then, we will introduce a sparse local tensor in Appendix K.2.2, which is adapted to MSGM and leads to good generative skills in practice. Finally, we propose a sparse nonlocal tensor in Appendix K.2.3 that involves more Brownian motions but meets the original assumptions A1 and A2 of our paper.

### K.2.1 Low-rank tensor

A simple choice of tensor with $d^2 = O(d^2)$ non-zero coefficients would be to take $d$ times the same dense random skew-symmetric matrix $G^1$ i.e.

$$ G_{i,j}^k = G_{i,j}^1 = \frac{1}{2}(M_{i,j}^1 - M_{j,i}^1), \quad (K.4) $$

$$ M_{i,j}^1 \stackrel{iid}{\sim} \mathcal{N}(0, 1). \quad (K.5) $$

Appendix I.2 provides an analytic solution for the forward diffusion in this case. Such a solution would be a 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 K.2.1 below shows that there is a rank deficiency, probably inducing the existence of non-rotation-invariant latent distribution, $p_\infty$, preventing MSGM sampling in practice. Indeed, numerically illustrated in dimension $d=4$, the analytic solution of Appendix I.2 shows a latent distribution intractable in practice. The latent distribution is not rotation-invariant and does depends on the initial direction distribution, $p_0^n$. It seems to be a direct consequence of the rank deficiency.

We conclude that low-rank tensors as in equation K.4 is not a suitable choice for MSGM. Moreover, it suggests that assumptions A1 and A4 as in Appendix K.1.2 are not sufficient for MSGM.

**Proposition K.2.1.** If $G$ is defined from equation K.4 and equation K.5, then, for any $x \in \mathbb{R}^d$, $\text{rank}(G(x)) \le 1$. Assumption A1 is verified, assumptions A2 and A3 are not for $d \ge 3$, and assumption A4 is verified almost surely if and only if the dimension $d$ if even. Moreover, we have $S = \frac{d}{2}G^1(G^1)^\top$ and $\mathbb{E}S = \frac{d(d-1)}{4}I_d$.

*Proof.* The tensor defined by equation K.4 and equation K.5 obviously matches the skew-symmetric condition A1.

For odd dimension $d$, $G^1$ – 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: \mathbb{R}^{d \times d} \to \mathbb{R}$, $M \mapsto \det(\frac{1}{2}(M - M^\top))$, which is non-zero since there exists invertible skew-symmetric matrices. As a non-zero polynomial, the set $\{M \in \mathbb{R}^{d,d} | \det(M - M^\top)\}$ 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

$$ \mathbb{P}(\det(G^1) = 0) = 0, \quad (K.6) $$
and so $G^1$ is invertible with full rank with probability 1. Thus

$$S = \frac{d}{2}G^1(G^1)^{\top}. \qquad (K.7)$$

is positive definite. Therefore, A4 is verified for even dimension $d$.

However, for any $d \ge 3$, neither conditions A2 nor condition A3 is satisfied. Indeed, for any $x \in \mathbb{R}^d \setminus \{0\}$, $\text{rank}(G(x)) = \text{rank}[G^1x, \dots, G^1x] \le 1$. This is expected since the diffusion involves a single Brownian motion (see Appendix I.2).

Since the entries in $M^1$ are independent standard normal Gaussian, we have $\mathbb{V}(G_{i,j}^1) = \frac{1}{4}(\mathbb{V}(M_{i,j}^1) + \mathbb{V}(M_{j,i}^1)) = \frac{1}{2}$. Then, $[G^1(G^1)^{\top}]_{ik} = \sum_{j=1}^d G_{ij}^1 G_{kj}^1$. Hence for $i = k$

$$\mathbb{E}[[(G^1(G^1)^{\top})]_{ik}] = \sum_{j=1}^{d} \mathbb{E}[(G_{ij}^1)^2] = \sum_{j \ne i} \frac{1}{2} = \frac{d-1}{2},$$

since $G_{ii}^1 = 0$. For $i \ne k$, $G_{ij}^1$ and $G_{kj}^1$ involve independent entries of $M^1$, leading to $\mathbb{E}[G_{ij}^1 G_{kj}^1] = 0$. As a consequence

$$\mathbb{E}[S] = \frac{d}{2}\mathbb{E}[G^1(G^1)^{\top}] = \frac{d(d-1)}{4}I_d \in \mathbb{R}^{d,d}. \quad \square$$

### K.2.2 Local sparse tensor

Let us define a tensor with only $2d = O(d)$ non-zero coefficients.

$$G_{i,j}^k = \begin{cases} 1 & \text{if } i = j - 1[d] = k \\ -1 & \text{if } i - 1[d] = j = k \\ 0 & \text{otherwise.} \end{cases}, \quad 1 \le i, j, k \le d, \qquad (K.8)$$

with $[d]$ stands for modulo $d$. It is built from a subset of the canonical basis of skew-symmetric matrices, keeping only $d$ matrices with most non-zero values close to the diagonal. It ensures a strong sparsity and a local structure for $x \to G^k x$.

The skew-symmetry assumption A1 is obviously fulfilled from the definition K.8. However, the strict rank condition assumption A2 is not in general. Fortunately, the assumptions A3 and A4 still hold. In particular, the Itô term matrix simplifies as shown by the following proposition.

We implemented this version of sparse tensor. For small dimension applications in Appendix M.6.2 and Appendix 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 tensor can complicate or even prevent MSGM applications. There, we obtained satisfactory results with local sparse tensor (see Figure 49 in Appendix M.7.2).

**Proposition K.2.2.** If $G$ is defined from equation K.8, then, for any $x \in (\mathbb{R} \setminus \{0\})^d$, $\text{rank}(G(x)) = d - 1$. Moreover, we have $S = I_d$ and the assumptions A1, A3, and A4 are verified.

*Proof.* For any $x \in (\mathbb{R} \setminus \{0\})^d$, we have

$$G(x) = [G^1x, \dots, G^d x] = \begin{pmatrix} x_2 & 0 & \cdots & 0 & -x_d \\ -x_1 & x_3 & \cdots & 0 & 0 \\ 0 & -x_2 & \cdots & \vdots & \vdots \\ \vdots & \vdots & \ddots & x_d & 0 \\ 0 & 0 & \cdots & -x_{d-1} & x_1 \end{pmatrix}. \qquad (K.9)$$

To simplify notations, all the indices in this proof will be defined modulo $d$. For instance, $x_{i+1}$ for $i = d$ stands for $x_1$.
For any $\mathbf{y} \in \mathbb{R}^d$, $\mathbf{G}(\mathbf{x})^T \mathbf{y} = 0$ ($\mathbf{x} \in \mathbb{R}^d$) if and only if, for all $i \le d$, $x_{i+1}y_i - x_iy_{i+1} = 0$ and $y_{i+1} = \frac{x_{i+1}}{x_i} y_i$. Finally,

$$y_i = \prod_{j=1}^{i-1} \frac{x_{i+1}}{x_i} y_1 = \frac{x_i}{x_1} y_1. \qquad (K.10)$$

Therefore, $\mathbf{y} \in \mathbb{R}^d$. Reciprocally, we can verify that $\mathbb{R}^d \subset \mathrm{Ker}(\mathbf{G}(\mathbf{x}))$. We conclude that $\mathrm{rank}(\mathbf{G}(\mathbf{x})) = d - \dim(\mathrm{Ker}(\mathbf{G}(\mathbf{x}))) = d - 1$.

To evaluate the matrix $\mathbf{S}$, we note that

$$\mathbf{G}^k = e_k e_{k+1}^\top - e_{k+1} e_k^\top, \qquad (K.11)$$

with $(e_k)_k$ the canonical basis of $\mathbb{R}^d$. Then,

$$\mathbf{S} = -\frac{1}{2} \sum_{k=1}^{d} (\mathbf{G}^k)^2, \qquad (K.12)$$

$$= -\frac{1}{2} \sum_{k=1}^{d} (e_k e_{k+1}^{\top} - e_{k+1} e_k^{\top})^2, \qquad (K.13)$$

$$= -\frac{1}{2} \sum_{k=1}^{d} (0 - e_k e_k^{\top} - e_{k+1} e_{k+1}^{\top} + 0), \qquad (K.14)$$

$$= \frac{1}{2}(\mathbf{I}_d + \mathbf{I}_d), \qquad (K.15)$$

$$= \mathbf{I}_d. \qquad (K.16)$$

□

### K.2.3 Non-local sparse tensor

We also propose another tensor with $d(d-1) = O(d^2)$ non-zero coefficients.

$$\mathbf{G}_{i,j}^{k,k'} = \begin{cases} 1 & \text{if } i-k'[d] = j=k, \quad 1 \le i,j,k \le d, \\ -1 & \text{if } i = j-k'[d] = k, \quad 1 \le k' \le \lceil \frac{d-1}{2} \rceil, \\ 0 & \text{otherwise.} \end{cases} \qquad (K.17)$$

where $\lceil \frac{d-1}{2} \rceil$ is the least integer greater than or equal to $\frac{d-1}{2}$. It is the canonical basis for skew-symmetric matrices. It ensures a relative sparsity and encodes a non-local structure for $\mathbf{x} \to \mathbf{G}^{k,k'} \mathbf{x}$.

Here, the sparse tensor $\mathbf{G}$ is of size $d \times d \times d(d-1)/2$ instead of $d \times d \times d$. Our theoretical framework differs slightly. The forward diffusion involves $d(d-1)/2$ one-dimensional Brownian motions. Consequently, the neural network, $\mathbf{a}_\theta(\mathbf{x}, s)$, approximating the scaled score, $\mathbf{G}(\mathbf{x})\nabla \log p_s(\mathbf{x})$, has $d(d-1)/2$ coefficients. The size of the neural network parameters $\theta$ can increase and may complicate the training procedure. An alternative could be to work with a neural network, $\mathbf{s}_\theta(\mathbf{x}, s)$, which approximates the true score, $\nabla \log p_s(\mathbf{x})$, having $d$ coefficients only.

This choice of tensor meets all the assumptions, including A1 and A2 as proofed below. However, because of the additional implementation complexity mentioned above, we postpone its numerical evaluation to MSGM for future work.

**Proposition K.2.3.** If $\mathbf{G}$ is defined from equation K.17, then, for any $\mathbf{x} \in \mathbb{R}^d \setminus \{0\}$, $\mathrm{rank}(\mathbf{G}(\mathbf{x})) = d - 1$. Moreover, we have $\mathbf{S} = \mathbf{I}_d$ and the assumptions A1, A2, A3, and A4 are verified.

*Proof.* For any $\mathbf{x} \in \mathbb{R}^d \setminus \{0\}$, we have

$$\mathbf{G}(\mathbf{x}) = [\mathbf{G}^{1,1}\mathbf{x}, \dots, \mathbf{G}^{d,\lceil \frac{d-1}{2} \rceil}\mathbf{x}] \qquad (K.18)$$

We already know that $\mathrm{Im}(\mathbf{G}(\mathbf{x})) \subset \mathbf{x}^{\perp}$ since $\mathbf{x}^{\top}\mathbf{G}(\mathbf{x}) = 0$. Now we assume that $\mathbf{y} \in \mathbf{x}^{\perp}$, and we define

$$Q = \frac{1}{\|\mathbf{x}\|^2} (\mathbf{y}\mathbf{x}^{\top} - x\mathbf{y}^{\top}). \qquad (K.19)$$
Applying on $x$, we get:

$$Qx = \frac{1}{\|x\|^2} (y \|x\|^2 - xy \cdot x) = y. \qquad (K.20)$$

Besides, $(G^{k,k'})_{k,k'}$ is the canonical basis of skew-symmetric matrices and $Q$ is skew-symmetric so there exists $\alpha \in \mathbb{R}^{\frac{d(d-1)}{2}}$ such that

$$Q = \sum_{k'=1}^{\left\lfloor \frac{d-1}{2} \right\rfloor} \sum_{k=1}^{d} \alpha_{k,k'} G^{k,k'} \qquad (K.21)$$

and thus

$$y = Qx = \sum_{k'=1}^{\left\lfloor \frac{d-1}{2} \right\rfloor} \sum_{k=1}^{d} \alpha_{k,k'} G^{k,k'} x = G(x) \alpha \in \operatorname{Im}(G(x)). \qquad (K.22)$$

We conclude that $\operatorname{Im}(G(x)) = x^\perp$ and $\operatorname{rank}(G(x)) = d-1$.

To evaluate the matrix $S$, we note that

$$G^{k,k'} = e_k e_{k+k'}^{\top} - e_{k+k'} e_k^{\top}, \qquad (K.23)$$

with $(e_k)_k$ the canonical basis of $\mathbb{R}^d$, and defining again all the indices modulo $d$.

$$S = -\frac{1}{2} \sum_{k'=1}^{\left\lfloor \frac{d-1}{2} \right\rfloor} \sum_{k=1}^{d} (G^{k,k'})^2, \qquad (K.24)$$

$$= -\frac{1}{2} \sum_{k'=1}^{\left\lfloor \frac{d-1}{2} \right\rfloor} \sum_{k=1}^{d} (e_k e_{k+k'}^{\top} - e_{k+k'} e_k^{\top})^2, \qquad (K.25)$$

$$= -\frac{1}{2} \sum_{k'=1}^{\left\lfloor \frac{d-1}{2} \right\rfloor} \sum_{k=1}^{d} (0 - e_k e_k^{\top} - e_{k+k'} e_{k+k'}^{\top} + 0), \qquad (K.26)$$

$$= \frac{1}{2}(I_d + I_d), \qquad (K.27)$$

$$= I_d. \qquad (K.28)$$

□

### K.3 Discussion about local and non local structure

The random tensor of Appendix J and the large sparse tensor of Appendix K.2.3 may be interpreted as non-local since $x \to G^k x$ changes coefficients $x_i$ of $x$ which are not sorted next to each other in $x$. For large dimension $d$, we believe that this can accelerate the convergence in comparison with local tensors, like the sparse tensor of Appendix K.2.2 or a 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 can spread directly in the whole state space 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 our conjecture. However, diffusion models in large dimension strongly rely on the powerful skills of convolutional neural networks (CNN), which have – up to attention layers – an intrinsic local structures. Accordingly, it may be difficult for a CNN to learn how to denoise a non-local noising process. More theoretical and experimental works would be needed to confirm this intuition. This is out 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):

$$d\mathbf{x}_t = \mathbf{f}_S(t, \mathbf{x}_t) dt + \tilde{\mathbf{G}}(t, \mathbf{x}_t) \circ dB_t, \quad (L.1)$$

where $\mathbf{f}_S: \mathbb{R} \times \mathbb{R}^d \to \mathbb{R}^d$ is the drift, $\tilde{\mathbf{G}}: \mathbb{R} \times \mathbb{R}^d \to \mathbb{R}^{d \times m}$ is the diffusion term, and $B_t$ is an $d$-dimensional Wiener process.

The following Runge-Kutta (RK) method (Kloeden et al., 1992) approximates the solution $\mathbf{x}_{n+1} \approx \mathbf{x}(t_{n+1})$ over the interval $[t_n, t_{n+1}]$, with time step $\Delta t = t_{n+1} - t_n$ and Wiener increment $\Delta B_n = B_{t_{n+1}} - B_{t_n}$:

$$K_1 = \mathbf{f}_S(t_n, \mathbf{x}_n) \Delta t + \tilde{\mathbf{G}}(t_n, \mathbf{x}_n) \Delta B_n, \quad (L.2)$$

$$K_2 = \mathbf{f}_S \left( t_n + \frac{\Delta t}{2}, \mathbf{x}_n + \frac{\mathbf{K}_1}{2} \right) \Delta t + \tilde{\mathbf{G}} \left( t_n + \frac{\Delta t}{2}, \mathbf{x}_n + \frac{\mathbf{K}_1}{2} \right) \Delta B_n, \quad (L.3)$$

$$K_3 = \mathbf{f}_S \left( t_n + \frac{\Delta t}{2}, \mathbf{x}_n + \frac{\mathbf{K}_2}{2} \right) \Delta t + \tilde{\mathbf{G}} \left( t_n + \frac{\Delta t}{2}, \mathbf{x}_n + \frac{\mathbf{K}_2}{2} \right) \Delta B_n, \quad (L.4)$$

$$K_4 = \mathbf{f}_S(t_n + \Delta t, \mathbf{x}_n + \mathbf{K}_3) \Delta t + \tilde{\mathbf{G}}(t_n + \Delta t, \mathbf{x}_n + \mathbf{K}_3) \Delta B_n, \quad (L.5)$$

$$\mathbf{x}_{n+1} = \mathbf{x}_n + \frac{1}{6}(\mathbf{K}_1 + 2\mathbf{K}_2 + 2\mathbf{K}_3 + \mathbf{K}_4). \quad (L.6)$$

This method leverages the structure of Stratonovich SDEs and their differential geometry properties. It is particularly well-suited to our SDE equation 3.1 with skew-symmetric noise and no Stratonovich drift.

### L.1.2 Renormalisation

Both our forward SDE equation 3.1 and backward SDE equation 2.2 preserve the solution norm $||\mathbf{x}_t||$. However, even the above Runge Kutta discretization can break this symmetry. To enforce it numerically, we normalize after each time step.

The final integration scheme is 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 is 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, we implemented a 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.

### L.2.1 Usual scheduling

In continuous time (Song et al., 2021), a convenient way is to make a change of variable, replacing the time $s$ by

$$z(s) = \int_{0}^{s} g^{2}(s')ds'. \quad (L.7)$$

with

$$g^2(s) = \frac{1}{2}\beta(s) = \frac{1}{2}\left(\beta_m + (\beta_M - \beta_m)\frac{s}{T}\right), \quad (L.8)$$

and $\beta_M > \beta_m > 0$. We first describe the hyperparameters values chosen in our numerical experiments and then explain how scheduling affects SGM and MSGM theories.
<caption>**Algorithm 2:** SRK4 for conservative Stratonovich SDEs with renormalization.</caption>

**Input:** Integration time $T$, number of time step $N_T$, initial condition $x_0$, drift $f_S$, diffusion $\tilde{G}$
1: $\Delta t \leftarrow \frac{T}{N_T}$; \{Time step\}
2: **for** $n = 0$ to $N_T - 1$ **do**
3: &nbsp;&nbsp;&nbsp;&nbsp;$\Delta B_n \sim \mathcal{N}(0, \Delta t I_d)$ \{Wiener increment\}
4: &nbsp;&nbsp;&nbsp;&nbsp;$t_n \leftarrow n\Delta t$
5: &nbsp;&nbsp;&nbsp;&nbsp;$K_1 \leftarrow f_S(t_n, x_n) \Delta t + \tilde{G}(t_n, x_n) \Delta B_n$
6: &nbsp;&nbsp;&nbsp;&nbsp;$K_2 \leftarrow f_S(t_n + \frac{\Delta t}{2}, x_n + \frac{K_1}{2}) \Delta t + \tilde{G}(t_n + \frac{\Delta t}{2}, x_n + \frac{K_1}{2}) \Delta B_n$
7: &nbsp;&nbsp;&nbsp;&nbsp;$K_3 \leftarrow f_S(t_n + \frac{\Delta t}{2}, x_n + \frac{K_2}{2}) \Delta t + \tilde{G}(t_n + \frac{\Delta t}{2}, x_n + \frac{K_2}{2}) \Delta B_n$
8: &nbsp;&nbsp;&nbsp;&nbsp;$K_4 \leftarrow f_S(t_n + \Delta t_n, x_n + K_3) \Delta t + \tilde{G}(t_n + \Delta t_n, x_n + K_3) \Delta B_n$
9: &nbsp;&nbsp;&nbsp;&nbsp;$\tilde{x}_{n+1} \leftarrow x_n + \frac{1}{6}(K_1 + 2K_2 + 2K_3 + K_4)$ \{Classical RK4 blend\}
10: &nbsp;&nbsp;&nbsp;&nbsp;$x_{n+1} \leftarrow \frac{\|x_0\|}{\|\tilde{x}_{n+1}\|} \tilde{x}_{n+1}$ \{Optional step : Enforce $\|x_{n+1}\| = \|x_0\|$\}
11: **end for**
12: **return** $x_{N_T}$ \{Approximation of $x_T$\}

Since we built our code from an existing one (https://github.com/CW-Huang/sdeflow-light,Huang et al. (2021)), by default we choose the values provided there for SGM scheduling: $\beta_m = 0.1$ and $\beta_M = 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, we found another SGM scheduling that works better and we use it instead. All scheduling hyperparameters are provided in the tables summarizing test cases in Appendix M.

For small time $s$, the time remapping is linear : $g^2(s) \sim_{s \to 0} \frac{1}{2}\beta_m$ and $z(s) \sim_{s \to 0} \frac{1}{2}\beta_m s$ whereas for large time, $g^2(s) \sim_{s \to T} \frac{1}{2}\beta_M$ and using the Taylor expansion around $T$, yielding

$$z(s) = z(T) + z'(t)(s - T) + o(s - T),$$

we find that

$$z(s) = \frac{1}{2} \left( \beta_m + \frac{1}{2}(\beta_M - \beta_m) \frac{s}{T} \right) s, \qquad (L.9)$$

$$= \frac{1}{2} \left( \frac{\beta_M + \beta_m}{2} T + \beta_M(s-T) \right) + o_T(s-T), \qquad (L.10)$$

$$\underset{s \to T}{\stackrel{\beta_M + \beta_m}{\rightleftharpoons}} \frac{T}{4}. \qquad (L.11)$$

As such, SGM forward and backward SDEs become:

$$d\vec{x}_s = -g^2(s)\vec{x}_s ds + \sqrt{2g(s)}d\vec{B}_s, \qquad (L.12)$$

$$d\vec{x}_t = g^2(T-t)\vec{x}_t dt + \sqrt{2g(T-t)}\left(a_\theta(T-t, \vec{x}_t)dt + \circ d\vec{B}_t\right), \qquad (L.13)$$

where $a_\theta(T-t, \vec{x}_t)$ approximates $\sqrt{2g(T-t)}\nabla \log p_{T-t}(\vec{x}_t)$. The backward SDE can be integrated with the Stochastic Runge-Kutta Algorithm 2 where

$$f_S(t, \vec{x}_t) = g^2(T-t)\vec{x}_t + \sqrt{2g(T-t)}a_\theta(T-t, \vec{x}_t) \quad \text{and} \quad \tilde{G}(t, \vec{x}_t) = \sqrt{2g(T-t)}. \qquad (L.14)$$

### L.2.2 Scheduling for MSGM

We propose a similar scheduling for MSGM. Scheduled forward and backward SDEs write:

$$d\vec{x}_s = g(s)\mathbf{G}(\vec{x}_s) \circ d\vec{B}_s, \qquad (L.15)$$

$$d\vec{x}_t = g(T-t)\mathbf{G}(\vec{x}_t) \left( a_\theta(T-t, \vec{x}_t)dt + \circ d\vec{B}_t \right), \qquad (L.16)$$
where $\mathbf{a}_{\theta}(T-t, \vec{\mathbf{x}}_t)$ approximates $g(T-t)\mathbf{G}(\vec{\mathbf{x}}_t)^{\mathrm{T}}\nabla \log p_{T-t}(\vec{\mathbf{x}}_t)$. Numerically, following Algorithm 2 we can integrate the forward SDE with

$$f_S(s, \vec{\mathbf{x}}_s) = 0 \quad \text{and} \quad \tilde{\mathbf{G}}(s, \vec{\mathbf{x}}_s) = g(s)\mathbf{G}(\vec{\mathbf{x}}_s), \qquad (L.17)$$

and the backward SDE with

$$f_S(t, \vec{\mathbf{x}}_t) = g(T-t)\mathbf{G}(\vec{\mathbf{x}}_t)\mathbf{a}_{\theta}(T-t, \vec{\mathbf{x}}_t) \quad \text{and} \quad \tilde{\mathbf{G}}(t, \vec{\mathbf{x}}_t) = g(T-t)\mathbf{G}(\vec{\mathbf{x}}_t). \qquad (L.18)$$

### L.2.3 Theoretical results

We can verify that our theoretical results remain under this time scheduling. The new Fokker-Planck equation is

$$\frac{\partial}{\partial s} p_s = \nabla_{\perp} \cdot (\frac{1}{2} g^2(s) \mathbf{\Sigma}(\mathbf{x}) \nabla_{\perp} p_s(\mathbf{x})). \qquad (L.19)$$

which can be rewritten as

$$\frac{\partial}{\partial s} p_s^g = \nabla_{\perp} \cdot (\frac{1}{2} \mathbf{\Sigma}(\mathbf{x}) \nabla_{\perp} p_s^g(\mathbf{x})). \qquad (L.20)$$

$$p_s^g = p_{z(s)}. \qquad (L.21)$$

Besides, for $0 \le s' \le z(T)$ for $\beta_M > \beta_m$ Taylor expansion at $s'_0 = \frac{\beta_M + \beta_m}{4} T$ yields

$$z^{-1}(s') = \frac{-\beta_M T + \sqrt{\beta_M^2 T^2 + 4T(\beta_M - \beta_m)z}}{\beta_M - \beta_m} \qquad (L.22)$$

$$= T + \frac{2}{\beta_M} (s' - s_0) \left( 1 + \underset{s' \to s'_0}{o}(1) \right), \qquad (L.23)$$

$$\underset{s' \to s'_0}{\longrightarrow} T. \qquad (L.24)$$

Therefore, from the convergence of $p_s^g$ (already proofed) we have the convergence of $p_{s'} = p_{z^{-1}(s')}^g$. The rate of convergence is still exponential:

$$\|p_{s'} - p_{\infty}\|_{L^2(\mathbb{R}^d)}^2 = \|p_{z^{-1}(s')}^g - p_{\infty}\|_{L^2(\mathbb{R}^d)}^2, \qquad (L.25)$$

$$\le \|p_{z^{-1}(0)}^g - p_{\infty}\|_{L^2(\mathbb{R}^d)}^2 \exp(-\alpha(\mathbf{G}, d)z^{-1}(s')), \qquad (L.26)$$

$$= \|p_0 - p_{\infty}\|_{L^2(\mathbb{R}^d)}^2 \exp(-\alpha(\mathbf{G}, d)T) \qquad (L.27)$$

$$\left( 1 - \frac{\alpha(\mathbf{G}, d)}{\beta_M/2} \left( s' - \frac{\beta_M + \beta_m}{4} T \right) \left( 1 + \underset{s' \to \frac{\beta_M + \beta_m}{4} T}{o}(1) \right) \right).$$

Besides the ELBO remains valid:

$$p_0(\mathbf{x}) \ge \mathcal{E}_{\infty}(\mathbf{x}) := \mathbb{E} \left[ \log p_0(\vec{\mathbf{x}}_T) | \vec{\mathbf{x}}_0 = \mathbf{x} \right] \qquad (L.28)$$

$$- \int_0^T \mathbb{E}_{\vec{\mathbf{x}}_s} \left[ \frac{1}{2} \| \mathbf{a}_{\theta}(\vec{\mathbf{x}}_s, s) \|^2 + \nabla \cdot (g(s) \mathbf{G}(\vec{\mathbf{x}}_s) \mathbf{a}_{\theta}(\vec{\mathbf{x}}_s, s)) \right] \Big|_{\vec{\mathbf{x}}_0 = \mathbf{x}} ds.$$

## L.3 Loss evaluation

Following the existing code (https://github.com/CW-Huang/sdeflow-light,Huang et al. (2021)), we sample final integration time $s$ of the forward SDEs uniformly on $[t_{\epsilon}, T]$ with $T = 1$ with $t_{\epsilon}$ small. According to Theorem 3.4.1, we consider the following SSM loss:

$$\mathcal{L}_{\text{SSM}}(\boldsymbol{\theta}) = \hat{\mathbb{E}}_{\vec{\mathbf{x}}_0} \hat{\mathbb{E}}_{s \sim \mathcal{U}[t_{\epsilon}, T]} \hat{\mathbb{E}}_{\vec{\mathbf{x}}_s | \vec{\mathbf{x}}_0} \hat{\mathbb{E}}_{v_s \sim \text{Rad}(d)} \mathcal{L}_{\text{SSM}}(s, \vec{\mathbf{x}}_s, g\mathbf{G}, \mathbf{a}_{\theta_n}, v_s), \qquad (L.29)$$

with

$$\mathcal{L}_{\text{SSM}}(s, \mathbf{x}, g\mathbf{G}, \mathbf{a}_{\theta_n}, \mathbf{v}) = \frac{1}{2} \|\mathbf{a}_{\theta}(\mathbf{x}, s)\|^2 + (\mathbf{v} \cdot \nabla)(g(s)\mathbf{G}(\mathbf{x})\mathbf{a}_{\theta}(\mathbf{x}, s)) \cdot \mathbf{v}, \qquad (L.30)$$
where $\hat{\mathbb{E}}$ is the averaged over the generated samples. Each training sample $\vec{x}_0$ of a batch is chosen randomly among the train set. For each of them, we sample one time $s$, one solution $\vec{x}_s$, and one slicing direction $v_s \sim \text{Rad}(d)$.

For SGM, we take $\mathbf{G} = \sqrt{2}$ in the above expressions and following Song et al. (2021), the solution $\vec{x}_s$ of the SGM scheduled forward SDE equation L.12 is

$$ \vec{x}_s = \exp(-\frac{1}{2}z(s))\vec{x}_0 + \sqrt{1 - \exp(-z(s))}\vec{x}_{\infty}, \qquad (L.31) $$

where $z(s) := \int_0^s g^2(s')ds'$ is given by equation L.9, $\vec{x}_0$ is chosen randomly among the train set and $\vec{x}_{\infty} \sim \mathcal{N}(0, I_d)$.

Unfortunately, to evaluate the MSGM loss, we cannot apply the same methodology, since, for $d > 2$ we are not aware of an analytic expression for the solution of the MSGM forward SDE, neither with nor without scheduling (equation L.15 and equation 3.1 respectively). We integrate that SDE numerically with the stochastic Runge-Kutta method with renormalization (see Appendix L.1.1 and Appendix L.1.2). Through this integration, we have to compute the solution $\vec{x}_{s_k}$ for many time steps $s_k := kT/N_T \in [0, T]$. Instead of sampling a random continuous time $s \sim \mathcal{U}([t_\epsilon, T])$, we choose a random discrete time as follow

$$ s \sim \mathcal{U}(I(t_\epsilon, T)) \quad \text{with} \quad I(t_\epsilon, T) = \{s_k|s_k = k\frac{T}{N_T}, k \in \{1, \dots, N_T\}, s_k \ge t_\epsilon\}. \qquad (L.32) $$

The numerical integration of the forward SDE implies a larger computational cost compared to SGM. Therefore, as explained in Appendix 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 Appendix I 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

### L.4.1 Spherical decomposition as an input layer

In line with our spherical decomposition equation 3.6, we add a fixed input layer to the network used in MSGM:

$$ a_{\theta}(\mathbf{x}, s) = \tilde{a}_{\theta}(\mathbf{x}/\|\mathbf{x}\|_{\epsilon}, \log \|\mathbf{x}\|_{\epsilon}, s), \quad \text{with} \quad \|\mathbf{x}\|_{\epsilon} := \|\mathbf{x}\| + \epsilon. \qquad (L.33) $$

The geometrical interpretation of Appendix H.3 also suggests that form.

For SGM, if not stated otherwise, we use a default architecture:

$$ a_{\theta}(\mathbf{x}, s) = \tilde{a}_{\theta}(\mathbf{x}, s). \qquad (L.34) $$

### L.4.2 Network Architecture for low-dimensional test cases (MLP)

Following the existing code (https://github.com/CW-Huang/sdeflow-light,Huang et al. (2021)), we parameterize the vector field $\tilde{a}_{\theta}: \mathbb{R}^d \times \mathbb{R} \to \mathbb{R}^d$ ($\tilde{d} = d$ or $d+1$) with a 4-layer MLP conditioned on an index $t \in \mathbb{R}$ by concatenation. Let $H = 128$ be the hidden width. For input $\mathbf{x} \in \mathbb{R}^d$, we form $\mathbf{h}_0 = [\mathbf{x}; t] \in \mathbb{R}^{d+1}$ and compute

$$ \begin{aligned} \mathbf{h}_1 &= \text{swish}(\mathbf{W}_1 \mathbf{h}_0 + \mathbf{b}_1), & \mathbf{W}_1 &\in \mathbb{R}^{H \times (\tilde{d}+1)}, \\ \mathbf{h}_2 &= \text{swish}(\mathbf{W}_2 \mathbf{h}_1 + \mathbf{b}_2), & \mathbf{W}_2 &\in \mathbb{R}^{H \times H}, \\ \mathbf{h}_3 &= \text{swish}(\mathbf{W}_3 \mathbf{h}_2 + \mathbf{b}_3), & \mathbf{W}_3 &\in \mathbb{R}^{H \times H}, \\ \mathbf{y} &= \mathbf{W}_4 \mathbf{h}_3 + \mathbf{b}_4, & \mathbf{W}_4 &\in \mathbb{R}^{d \times H}, \end{aligned} $$

with $(\text{swish}(z))_i = z_i \sigma(z_i)$ and $\sigma$ the logistic sigmoid. We set $\tilde{a}_{\theta}(\mathbf{x}, t) = \mathbf{y} \in \mathbb{R}^d$. No residual connections, normalization, or dropout are used. Table 1 summarizes the hyperparameters of this default architecture.
<table><caption>Table 1: MLP architecture hyperparameters.</caption><thead><tr><td>Hyperparameter</td><td>Value</td></tr></thead><tbody><tr><td>Input dimension</td><td>d̃ = d or d + 1</td></tr><tr><td>Index dimension</td><td>1</td></tr><tr><td>Hidden width</td><td>128</td></tr><tr><td>Depth</td><td>3 hidden layers</td></tr><tr><td>Activation</td><td>Swish (x ↦ xσ(x))</td></tr><tr><td>Output dimension</td><td>d</td></tr><tr><td>Output layer</td><td>Linear</td></tr><tr><td>Residual connections</td><td>None</td></tr><tr><td>Normalization / Dropout</td><td>None</td></tr></tbody></table>

### L.4.3 Network Architecture for high-dimensional test cases (UNet for 32 × 32 vorticity fields)

For high-dimensional experiments of Appendix M.7.2, we model the score field $\tilde{a}_{\theta}(\mathbf{x}, t)$ using a 2D UNet operating on images $\mathbf{x}'$ of size $H \times W$ representing vorticity snapshots ($H = 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 $\mathbf{x}' \in \mathbb{R}^{1 \times H \times W}$ are reshaped into vectors $\mathbf{x} \in \mathbb{R}^d$ with $d = HW$ or vectors are reshaped as a one-channel images $\mathbf{x}' \in \mathbb{R}^{1 \times H \times W}$.

**Optional spherical premodule.** When enabled, we apply the spherical decomposition of Appendix L.4.1:

$$ (\mathbf{x}_{\epsilon}^{n}, \log \|\mathbf{x}\|_{\epsilon}) = \text{NormalizeLogRadius}(x), \quad \mathbf{x}_{\epsilon}^{n} = \frac{x}{\|\mathbf{x}\|_{\epsilon}}. $$

The normalized field $\mathbf{x}_{\epsilon}^{n}$ is passed to the UNet, while $\log \|\mathbf{x}\|_{\epsilon}$ is embedded through a small MLP and added to the temporal embedding, giving a conditioning mechanism analogous to the time embedding of diffusion models.

**UNet backbone.** The core architecture follows the DDPM UNet of Dhariwal & Nichol (2021): a fully convolutional encoder-decoder with skip connections, residual blocks, and optional attention at intermediate resolutions. We use one input channel and one output channel (vorticity). Let $C_0$ denote the base width. The feature width at resolution level $k$ is $C_0 m_k$ where $C_0$ is the base channel width and $m_k$ is the channel multiplier.

The UNet receives $(\mathbf{x}, t)$ (and optionally $\log \|\mathbf{x}\|_{\epsilon}$) and computes:

$$ \tilde{a}_{\theta}(\mathbf{x}, t) = \text{UNet}_{\theta}(\text{reshape}(x), \text{Emb}(t) + \text{Emb}_{\log}(\log \|\mathbf{x}\|_{\epsilon})), $$

followed by flattening back to dimension $d$ if needed.

<table><caption>Table 2: UNet architecture hyperparameters.</caption><thead><tr><td>Component</td><td>Setting</td></tr></thead><tbody><tr><td>Input / output channels</td><td>1</td></tr><tr><td>Input resolution</td><td>$H = W \in \{16, 32\}$</td></tr><tr><td>Base channels width $C_0$</td><td>32</td></tr><tr><td>Channel multipliers</td><td>(1, 2, 4)</td></tr><tr><td>Residual blocks per stage</td><td>2</td></tr><tr><td>Attention resolutions</td><td>8 × 8 and 4 × 4 (for 16 × 16 input)</td></tr><tr><td>Activation</td><td>SiLU</td></tr><tr><td>Time embedding</td><td>sinusoidal + MLP</td></tr><tr><td>Log-norm conditioning</td><td>optional MLP added to time embedding</td></tr><tr><td>Dropout</td><td>0</td></tr><tr><td>Upsampling / downsampling</td><td>convolutional</td></tr><tr><td>Output</td><td>1-channel vorticity field</td></tr></tbody></table>

This UNet is used as a drop-in replacement for the small-dimensional MLP of Appendix L.4.2, enabling MSGM/SGM to scale to image-like vorticity fields up 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 to being heavy. For distributions with lighter tails such as Gaussian ones, SGM and MSGM produce similar results, except for a small number of backward time steps where SGM can become unstable. MSGM is more robust in this aspect.

Our code can be found here: https://anonymous.4open.science/r/MSGM-submission-6E1E and the preprocessed vorticity data we used in Appendix M.7 can be found here: https://anonymous.4open.science/r/MSGM-data-6E64.

## 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 used in the experiments (see tables 3, 4, 5, 6, and 7). All are performed on CPU. In addition, we additional cover a high-dimensional application with imagine processing, see section Appendix M.7 with a GPU A40 NVL with 48 Go of VRAM.

## 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 can significantly reduce the number of conditioning of the covariance of the data set, and therefore facilitate the SGM (Guth et al., 2022). Generated data are then re-scaled for plots and other 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 $||\vec{x}_0||^2$. The modified energy has no 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.

## 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^{(i)}) \in (\mathbb{R})^N$ and $Y = (y^{(i)}) \in (\mathbb{R})^N$ samples of random variables $X$ and $Y$ respectively, we define MMD$(x, y)$ as:

$$ \text{MMD}^2(\mathcal{X}, \mathcal{Y}) = \frac{1}{N^2} \sum_{i,j=1}^{N} \left( k(x^{(i)}, x^{(j)}) - 2k(x^{(i)}, y^{(j)}) + k(y^{(i)}, y^{(j)}) \right), \quad (\text{M.1}) $$

$$ k(u, v) = \exp(-\|u - v\|^2). \quad (\text{M.2}) $$

If $\mathcal{X}^{\text{test}}$ is the test set and $\mathcal{X}^{\text{gen}}$ our generated ensemble, MMD($\mathcal{X}^{\text{test}}, \mathcal{X}^{\text{gen}}$) is a metric of the precision of our generated ensemble and hence our AI generative algorithm. A small MMD means close distributions. However, MMD is a relative metric. So we compare MMD($\mathcal{X}^{\text{test}}, \mathcal{X}^{\text{gen}}_{\text{SGM}}$), MMD($\mathcal{X}^{\text{test}}, \mathcal{X}^{\text{gen}}_{\text{MSGM}}$) and MMD($\mathcal{X}^{\text{test}}, \mathcal{X}^{\text{train}}$) where $\mathcal{X}^{\text{gen}}_{\text{SGM}}$ and $\mathcal{X}^{\text{gen}}_{\text{MSGM}}$ are generated from SGM and MSGM respectively, and $\mathcal{X}^{\text{train}}$ is the train set. MMD($\mathcal{X}^{\text{test}}, \mathcal{X}^{\text{train}}$) 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 an additional significant computational cost compared to SGM, and hence a slower training procedure. This cost scales linearly in $N_T$ due to the "for" loop in time. Empirically, it appears to scale as $\zeta = \sqrt{d} N_T / 2^4$ (not shown), probably due to the vectorized $d \times d \times d$ tensor products involved in each integration time step. In most of the numerical experiments below, $N_T = 2^4$ and thus $\zeta = \sqrt{d}$. The SGM iteration steps are $\zeta$ times faster than the MSGM iteration steps. Consequently, the number of iterations for the SGM is $\max(1, \lfloor \zeta \rfloor)$ times larger than the number of iterations for the MSGM. As such, we can compare the results of SGM and MSGM at a similar training cost. By convention, we take the number of iterations for SGM
<table><caption>Table 3: Swiss roll test case: parameters of the nominal numerical experiments.</caption><thead><tr><th>Parameter</th><th>SGM</th><th>MSGM</th></tr></thead><tbody><tr><td>Dimension <i>d</i></td><td>2</td><td>2</td></tr><tr><td>Number of used training data points (<i>M</i>)</td><td>2<sup>20</sup> × 256</td><td>2<sup>20</sup> × 256</td></tr><tr><td>Number of test data points</td><td>10<sup>4</sup></td><td>10<sup>4</sup></td></tr><tr><td>Reference number of ADAMS steps</td><td>2<sup>20</sup></td><td>2<sup>20</sup></td></tr><tr><td>Number of ADAMS steps (<i>N</i><sub>iter</sub>)</td><td>2<sup>20</sup></td><td>2<sup>20</sup></td></tr><tr><td>CPU time / ADAMS steps (in ms)</td><td>4</td><td>3</td></tr><tr><td>Batch size</td><td>256</td><td>256</td></tr><tr><td>Number of time steps (forward) <i>N</i><sub>T</sub><sup>f</sup></td><td>1</td><td>16</td></tr><tr><td>Number of time steps (backward) <i>N</i><sub>T</sub><sup>b</sup></td><td>16</td><td>16</td></tr><tr><td>β<sub>min</sub></td><td>0.1</td><td>0.1</td></tr><tr><td>β<sub>max</sub></td><td>20</td><td>20</td></tr><tr><td><i>t</i><sub>ε</sub></td><td>10<sup>-3</sup></td><td>10<sup>-3</sup></td></tr><tr><td>Learning rate</td><td>10<sup>-3</sup></td><td>10<sup>-3</sup></td></tr><tr><td>Neural network architecture</td><td>default</td><td>spherical (equation L.33)</td></tr><tr><td>MMD (MMD(train)= 0.9 × 10<sup>-2</sup>)</td><td>1.9 × 10<sup>-2</sup></td><td>0.9 × 10<sup>-2</sup></td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>Figure 8: Pair plots of generated data (orange dots) compared to ground truth data (blue dots) with the SGM (left) and MSGM (right) for Swissroll data. 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.</figcaption></figure>

as a reference and refer to it as the reference number of iterations. Summary tables 3, 4, 5, 6, and 7) provide the values for the reference number of iterations, the true number of iterations, and the execution time per ADAMS step.

## M.4 Swiss roll

We first illustrate our method with the Swiss roll distribution. It is a simple two-dimensional distribution: https://homepages.ecs.vuw.ac.nz/~marslast/Code/Ch6/lle.py. 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 8. However, the diffusion distribution $p_s$ differs from Figure 9 to Figure 10. In particular, latent distributions are completely different.

Figure 11 illustrates the convergence of the SGM and MSGM approaches as a function of the reference number of ADAMS iterations and as a function of 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.
<figure><img src="image_13.png" /><figcaption>Figure 9: Evolution of the solution log-pdf $\log(p_s(x_1, x_2))$ of SGM forward SDE (top) and backward SDE (bottom) for Swiss roll data.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 10: Evolution of the solution log-pdf $\log(p_s(x_1, x_2))$ of MSGM forward SDE (top) and backward SDE (bottom) for Swiss roll data.</figcaption></figure>

## M.5 Anisotropic Gaussian distribution

For a complete numerical analysis, we compare SGM and MSGM on correlated Gaussian data $x_0 \sim N(0, AA^T)$, with a fixed matrix, $A$, initialized with i.i.d. coefficients $A_{i,j} \sim N(0, 1)$. For 32 time steps backward, the pairplots in Figures 12, 13, and 14 present similar generative skills, but for 8 or 16 time steps backward, only MSGM gives good results. For 8 time steps backwards, MSGM still provides a good distribution, whereas the SGM backward SDE completely diverges. Figures 15, 16, and 17 also highlight this and show that the converged dynamics of the pdf $p_s$ differs between SGM and MSGM. Figure 18 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 in this extreme case. Note that due to the absence of finite moments, convergence in common metrics such as Wasserstein-$p$ or total variation is not well defined.

<figure><img src="image_15.png" /><figcaption>Figure 11: Convergence of MMD (mean and 80% confidence interval) for Swiss roll distribution as a function of reference number of ADAMS iterations (left) and as a function of number of time steps for integrating the backward SDE (right).</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>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 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.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>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 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.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>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.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 15: Evolution of the solution log-pdf log($p_s(x_1, x_2)$) of SGM forward SDE (top) and backward SDE (bottom, with 8 time steps) for Gaussian data.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 16: Evolution of the solution log-pdf log($p_s(x_1, x_2)$) of SGM forward SDE (top) and backward SDE (bottom, with 32 time steps) for Gaussian data.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Figure 17: Evolution of the solution log-pdf log($p_s(x_1, x_2)$) of MSGM forward SDE (top) and backward SDE (bottom, with 16 time steps) for Gaussian data.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Figure 18: Convergence of MMD (mean and 80% confidence interval) for the Gaussian data as a function of number of time steps for integrating the backward SDE $N_T^b$ (top), as a function of reference number of ADAMS iterations (middle) for $N_T^b = 8$ (left) and $N_T^b = 32$ (right), and as a function of dimension (bottom) for $N_T^b = 8$ (left) and $N_T^b = 32$ (right).</figcaption></figure>
<table><caption>Table 4: Gaussian test case: parameters of the nominal numerical experiments.</caption><thead><tr><th>Parameter</th><th>SGM</th><th>MSGM</th></tr></thead><tbody><tr><td>Dimension <i>d</i></td><td>16</td><td>16</td></tr><tr><td>Number of used training data points (<i>M</i>)</td><td>1048576 × 256</td><td>262144 × 256</td></tr><tr><td>Number of test data points</td><td>10<sup>4</sup></td><td>10<sup>4</sup></td></tr><tr><td>Reference number of ADAMS steps</td><td>2<sup>20</sup></td><td>2<sup>20</sup></td></tr><tr><td>Number of ADAMS steps (<i>N</i><sub>iter</sub>)</td><td>2<sup>20</sup> = 1048576</td><td>262144</td></tr><tr><td>CPU time / ADAMS steps (in ms)</td><td>3</td><td>23</td></tr><tr><td>Batch size</td><td>256</td><td>256</td></tr><tr><td>Number of time steps (forward) <i>N</i><sub>T</sub><sup>f</sup></td><td>1</td><td>16</td></tr><tr><td>Number of time steps (backward) <i>N</i><sub>T</sub><sup>b</sup></td><td>16</td><td>16</td></tr><tr><td>β<sub>min</sub></td><td>0.1</td><td>0.1</td></tr><tr><td>β<sub>max</sub></td><td>20</td><td>20</td></tr><tr><td><i>t</i><sub>ε</sub></td><td>10<sup>-3</sup></td><td>10<sup>-3</sup></td></tr><tr><td>Learning rate</td><td>10<sup>-3</sup></td><td>10<sup>-3</sup></td></tr><tr><td>Neural network architecture</td><td>default</td><td>spherical (equation L.33)</td></tr><tr><td>MMD (MMD(train)=1.5 × 10<sup>-2</sup>)</td><td>11 × 10<sup>-2</sup></td><td>2.5 × 10<sup>-2</sup></td></tr></tbody></table>

### M.6.1 Vector of independent Cauchy variables

<table><caption>Table 5: Vector of independent Cauchy variables: parameters of the nominal numerical experiments.</caption><thead><tr><th>Parameter</th><th>SGM</th><th>MSGM</th></tr></thead><tbody><tr><td>Dimension <i>d</i></td><td>2</td><td>2</td></tr><tr><td>Number of used training data points (<i>M</i>)</td><td>2<sup>20</sup> × 256</td><td>209715 × 256</td></tr><tr><td>Number of test data points</td><td>10<sup>5</sup></td><td>10<sup>5</sup></td></tr><tr><td>Reference number of ADAMS steps</td><td>2<sup>20</sup></td><td>2<sup>20</sup></td></tr><tr><td>Number of ADAMS steps (<i>N</i><sub>iter</sub>)</td><td>2<sup>20</sup> = 1048576</td><td>209715</td></tr><tr><td>CPU time / ADAMS steps (in ms)</td><td>3</td><td>27</td></tr><tr><td>Batch size</td><td>256</td><td>256</td></tr><tr><td>Number of time steps (forward) <i>N</i><sub>T</sub><sup>f</sup></td><td>1</td><td>64</td></tr><tr><td>Number of time steps (backward) <i>N</i><sub>T</sub><sup>b</sup></td><td>128</td><td>128</td></tr><tr><td>β<sub>min</sub></td><td>0.1</td><td>0.1</td></tr><tr><td>β<sub>max</sub></td><td>20</td><td>0.4</td></tr><tr><td><i>t</i><sub>ε</sub></td><td>10<sup>-3</sup></td><td>10<sup>-3</sup></td></tr><tr><td>Learning rate</td><td>10<sup>-3</sup></td><td>10<sup>-3</sup></td></tr><tr><td>Neural network architecture</td><td>default</td><td>spherical (equation L.33)</td></tr><tr><td>MMD (MMD(train)=2.8 × 10<sup>-3</sup>)</td><td>7.5 × 10<sup>-3</sup></td><td>3.3 × 10<sup>-3</sup></td></tr></tbody></table>

We first illustrate our method with a vector of independent Cauchy variables: $x_0 = x_{Ca}$ with $x_{Ca}$ defined by equation 6.2 with scale parameter $\gamma = 1/50$. As expected, Figure 19 and Figure 20 confirm that SGM does not reproduce fat tails unlike MSGM. Moreover, SGM misaligns the far data points that have the coordinate $x_1 < -3$. An explanation of the superior skills is the similarity between the data distribution and the latent distribution in MSGM: a property not shared by SGM, as illustrated in Figures 21 and 22.

Figure 19 compares MSGM with fast scheduling ($\beta_m = 0.1, \beta_M = 0.4$) and a neural network architecture based on spherical decomposition equation L.33 with SGM with default scheduling ($\beta_m = 0.1, \beta_M = 20$) and default neural network architecture. For a fair comparison of MSGM, we complement our numerical analysis with Figures 23-27: we test SGM with fast and default schedulings, and with both spherical-decomposition-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. In contrast, the network architecture with spherical decomposition does improve the SGM sampling procedure, especially for distribution tails. However, even with this architecture,
<figure><img src="image_23.png" /><figcaption>Figure 19: Pair plots of generated data (orange dots) compared to ground truth data (blue dots) with the SGM (left) and MSGM (right) for a 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) are superimposed.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>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.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Figure 21: Evolution of the solution log-pdf $\log(p_s(x_1, x_2))$ of MSGM forward SDE (top) and backward SDE (bottom) for a vector of two independent Cauchy variables, with fast scheduling: $\beta_m = 0.1$, $\beta_M = 0.4$ and our neural network architecture based on spherical decomposition equation L.33.</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>Figure 22: Evolution of the solution log-pdf log($p_s(x_1, x_2)$) of SGM forward SDE (top) and backward SDE (bottom) for a vector of two independent Cauchy variables, with default scheduling: $\beta_m = 0.1$, $\beta_M = 20$ and default neural network architecture.</figcaption></figure>

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 $x_2 < -3$. Third, outside the $x_1$ and $x_2$ axes, SGM generates too few samples close to the origin (say points $\boldsymbol{x}$ with $\|\boldsymbol{x}\|_2 > 2$ and $\|\boldsymbol{x}\|_1 < 2$).

For Cauchy distributions, we still compare MMD values. However, it is not well defined mathematically 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 $10^4$ samples to compute an approximate MMD. Other quantities of interest can also be utilized, such as the survival function $t \mapsto \mathbb{P}(\|\boldsymbol{x}\| > R)$, illustrated in Figure 20. As expected MSGM clearly outperforms SGM on this metric. Indeed by construction our learning method is robust in terms of the radial distribution $\|\boldsymbol{x}\|$ 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.

### M.6.2 Vector of correlated Cauchy variables

<table><caption>Table 6: Vector of correlated Cauchy variables: parameters of the nominal numerical experiments.</caption><thead><tr><th>Parameter</th><th>SGM</th><th>MSGM</th></tr></thead><tbody><tr><td>Dimension <i>d</i></td><td>4</td><td>4</td></tr><tr><td>Number of used training data points (<i>M</i>)</td><td>2<sup>20</sup> × 256</td><td>2<sup>20</sup> × 256</td></tr><tr><td>Number of test data points</td><td>10<sup>5</sup></td><td>10<sup>5</sup></td></tr><tr><td>Reference number of ADAMS steps</td><td>2<sup>20</sup></td><td>2<sup>20</sup></td></tr><tr><td>Number of ADAMS steps (<i>N</i><sub>iter</sub>)</td><td>2<sup>20</sup></td><td>2<sup>20</sup></td></tr><tr><td>CPU time / ADAMS steps (in ms)</td><td>3</td><td>45</td></tr><tr><td>Batch size</td><td>256</td><td>256</td></tr><tr><td>Number of time steps (forward) <i>N</i><sub>T</sub><sup><i>f</i></sup></td><td>1</td><td>128</td></tr><tr><td>Number of time steps (backward) <i>N</i><sub>T</sub><sup><i>b</i></sup></td><td>128</td><td>128</td></tr><tr><td>β<sub>min</sub></td><td>0.1</td><td>0.01</td></tr><tr><td>β<sub>max</sub></td><td>20</td><td>1</td></tr><tr><td><i>t</i><sub>ε</sub></td><td>10<sup>-4</sup></td><td>10<sup>-4</sup></td></tr><tr><td>Learning rate</td><td>10<sup>-3</sup></td><td>10<sup>-3</sup></td></tr><tr><td>Neural network architecture</td><td>default</td><td>spherical (equation L.33)</td></tr><tr><td>MMD (MMD(train)=3.5 × 10<sup>-3</sup>)</td><td>11.2 × 10<sup>-3</sup></td><td>5.2 × 10<sup>-3</sup></td></tr></tbody></table>

To address dimensionality issues, we consider the correlated Cauchy variables already presented in Section 6.1. In terms of survival function, MSGM is as expected more accurate than SGM (see Figure 28).

As for independent Cauchy variables, we present complementary numerical experiments with different scheduling and different neural network architectures in Figures 29 and 30. Figures 31-35 unveil the corresponding diffusion dynamics from $s = 0$ to $s = T$ and from $t = 0$ to $t = T$. Again, the neural network architecture based on spherical decomposition significantly improves the SGM generative
<figure><img src="image_27.png" /><figcaption>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 a default scheduling: $\beta_m = 0.1$, $\beta_M = 20$. Left plots correspond to our neural network architecture based on spherical decomposition equation L.33 whereas the right plot correspond to default neural network architecture. 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) are superimposed.</figcaption></figure>
<figure><img src="image_28.png" /><figcaption>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: $\beta_m = 0.1$, $\beta_M = 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 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) are superimposed.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>Figure 25: Evolution of the solution log-pdf log($p_s(x_1, x_2)$) of SGM forward SDE (top) and backward SDE (bottom) for a vector of two independent Cauchy variables, with default scheduling: $\beta_m = 0.1$, $\beta_M = 20$ and our neural network architecture based on spherical decomposition equation L.33.</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>Figure 26: Evolution of the solution log-pdf log($p_s(x_1, x_2)$) of SGM forward SDE (top) and backward SDE (bottom) for a vector of two independent Cauchy variables, with fast scheduling: $\beta_m = 0.1$, $\beta_M = 0.4$ and default neural network architecture.</figcaption></figure>
<figure><img src="image_31.png" /><figcaption>Figure 27: Evolution of the solution log-pdf log($p_s(x_1, x_2)$) of SGM forward SDE (top) and backward SDE (bottom) for a vector of two independent Cauchy variables, with fast scheduling: $\beta_m = 0.1$, $\beta_M = 0.4$ and our neural network architecture based on spherical decomposition equation L.33.</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>Figure 28: 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 4 correlated Cauchy variables.</figcaption></figure>

skills but MSGM remains a 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 can wonder if the poorer results of SGM would improve for a larger number of ADAMS iterations. To answer this question, we run longer experiments with $2^{24} = 16777216$ ADAMS iterations. Figures 36, 37, and 4a show that MSGM slightly improves with an 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. we proportionally reduce the number of ADAMS iterations for MSGM in order to make the CPU training time of SGM and MSGM similar (see Appendix M.3 for details). For SGM with very large number of iterations ($2^{24} = 16777216$), we use a larger number of time steps ($2^9 = 512$) for the backward SDE to prevent all samples generated by SGM to diverge.

## M.7 Vorticity field from Particle Image Velocimetry measurements

Particle Image Velocimetry (PIV) is 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 a benchmark configuration : a wake flow at Reynolds number $Re = 3900$ created by a circular cylinder embedded in a mean stream (Parnaureau et al., 2008). We compute the two-dimensional curl of the velocity. Named vorticity, it is presented in Figure 39.

### M.7.1 Low-dimensional test case: vorticity evaluated on several spatial points

To reduce the dimension $d$ of the data, we severely crop the vorticity images and subsample them spatially, keeping only 4 × 4 pixels by images as illustrated by Figure 40. Once reshaped as a vector, each small image represents a data point of dimension 16. If we choose a dimension $d \le 16$, we just keep the first $d$ coefficients of the vector. For this experimental dataset, we investigate the influence
<figure><img src="image_33.png" /><figcaption>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 to a default scheduling: $\beta_m = 0.1, \beta_M = 20$. Left and top plots correspond to our neural 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.</figcaption></figure>
<figure><img src="image_34.png" /><figcaption>Figure 30: Generated data (orange lines and dots) compared to ground truth data (blue lines and dots) with the SGM for a vector of 4 correlated Cauchy variables with a fast scheduling: $\beta_m = 0.01$, $\beta_M = 1$. The left plot corresponds to our neural network architecture based on spherical decomposition equation L.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.</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>Figure 31: Evolution of the solution log-pdf $\log(p_s(x_1, x_3))$ of MSGM forward SDE (top) and backward SDE (bottom) for a vector of 4 correlated Cauchy variables, with fast scheduling: $\beta_m = 0.01$, $\beta_M = 1$ and our neural network architecture based on spherical decomposition equation L.33.</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>Figure 32: Evolution of the solution log-pdf $\log(p_s(x_1, x_3))$ of SGM forward SDE (top) and backward SDE (bottom) for a vector of 4 correlated Cauchy variables, with default scheduling: $\beta_m = 0.1$, $\beta_M = 20$ and default neural network architecture.</figcaption></figure>
<figure><img src="image_37.png" /><figcaption>Figure 33: Evolution of the solution log-pdf log($p_s(x_1, x_3)$) of SGM forward SDE (top) and backward SDE (bottom) for a vector of 4 correlated Cauchy variables, with default scheduling: $\beta_m = 0.1$, $\beta_M = 20$ and our neural network architecture based on spherical decomposition equation L.33.</figcaption></figure>

<figure><img src="image_38.png" /><figcaption>Figure 34: Evolution of the solution log-pdf log($p_s(x_1, x_3)$) of SGM forward SDE (top) and backward SDE (bottom) for a vector of 4 correlated Cauchy variables, with fast scheduling: $\beta_m = 0.01$, $\beta_M = 1$ and default neural network architecture.</figcaption></figure>

of the amount of data available for learning. Our default experiments will train the models with $2^{10} = 1024$ data points only.

As seen previously in Section 6.2, MSGM is more robust in low-data mode and better represents rare events, as also confirmed by the survival function Figure 41. We explain it by a latent distribution close to the data distribution as illustrated in Figure 42.

For a fair numerical comparison, we 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 fit well the data distribution both with dense and with sparse tensor, **G**.

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 in Figures 46 and 47. Again, MMD may not be the best tool for studying rare events. We can observe some tendencies, but definite conclusions may not be obtained from those convergence plots. For a very small training set ($2^6 = 64$ data points), both SGM and MSGM

<figure><img src="image_39.png" /><figcaption>Figure 35: Evolution of the solution log-pdf log($p_s(x_1, x_3)$) of SGM forward SDE (top) and backward SDE (bottom) for a vector of 4 correlated Cauchy variables, with fast scheduling: $\beta_m = 0.01$, $\beta_M = 1$ and our neural network architecture based on spherical decomposition equation L.33.</figcaption></figure>
<figure><img src="image_40.png" /><figcaption>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^{24} = 16777216$ ADAMS iterations for a 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) are superimposed.</figcaption></figure>

<figure><img src="image_41.png" /><figcaption>Figure 37: Evolution of the solution log-pdf $\log(p_s(x_1, x_3))$ of SGM forward SDE (top) and backward SDE (bottom, $2^{24} = 16777216$ ADAMS iterations and $2^9 = 512$ time steps) for a vector of 4 correlated Cauchy variables, with default scheduling: $\beta_m = 0.1$, $\beta_M = 20$ and default neural network architecture.</figcaption></figure>

<figure><img src="image_42.png" /><figcaption>Figure 38: 2D2C PIV velocity field: velocity component along *x* (left) and velocity component along *y* (right).</figcaption></figure>
<figure><img src="image_43.png" /><figcaption>Figure 39: The full two-dimensional vorticity (left) and a zoom (right) of a PIV field</figcaption></figure>

<figure><img src="image_44.png" /><figcaption>Figure 40: Spatial cropping and spatial subsampling of a vorticity field to obtain a data sample at low dimension $d = 16$.</figcaption></figure>

<figure><img src="image_45.png" /><figcaption>Figure 41: Survival function of generated data (orange line) compared to ground truth data (blue line) with the SGM (left) and MSGM (right) trained on 1024 16-dimensional data points representing PIV-based vorticity fields.</figcaption></figure>
<figure><img src="image_46.png" /><figcaption>Figure 42: Evolution of the solution log-pdf log($p_s(x_1, x_3)$) of MSGM forward SDE (top) and backward SDE (bottom) for the vorticity images distribution, with nominal scheduling ($\beta_m = 0.025$, $\beta_M = 5.0$) and our neural network architecture based on spherical decomposition equation L.33.</figcaption></figure>

<figure><img src="image_47.png" /><figcaption>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 to 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-histogram of ground truth data (blue line) and logarithm of the pdf KDE estimation of generated data (orange line) are superimposed.</figcaption></figure>
<table><caption>Table 7: Low-dimensional vorticity test case: parameters of the nominal numerical experiments.</caption><thead><tr><th>Parameter</th><th>SGM</th><th>MSGM</th></tr></thead><tbody><tr><td>Dimension d</td><td>16</td><td>16</td></tr><tr><td>Number of used training data points (M)</td><td>2<sup>10</sup> = 1024</td><td>1024</td></tr><tr><td>Number of test data points</td><td>6476</td><td>6476</td></tr><tr><td>Reference number of ADAMS steps</td><td>2<sup>20</sup></td><td>2<sup>20</sup></td></tr><tr><td>Number of ADAMS steps (N<sub>iter</sub>)</td><td>2<sup>20</sup> = 1048576</td><td>262144</td></tr><tr><td>CPU time / ADAMS steps (in ms)</td><td>4</td><td>32</td></tr><tr><td>Batch size</td><td>256</td><td>256</td></tr><tr><td>Number of time steps (forward) N<sub>T</sub><sup>f</sup></td><td>1</td><td>16</td></tr><tr><td>Number of time steps (backward) N<sub>T</sub><sup>b</sup></td><td>8</td><td>8</td></tr><tr><td>β<sub>min</sub></td><td>0.025</td><td>0.025</td></tr><tr><td>β<sub>max</sub></td><td>5</td><td>5</td></tr><tr><td>t<sub>ε</sub></td><td>10<sup>-4</sup></td><td>2.5 × 10<sup>-5</sup></td></tr><tr><td>Learning rate</td><td>10<sup>-3</sup></td><td>10<sup>-3</sup></td></tr><tr><td>Neural network architecture</td><td>default</td><td>spherical (equation L.33)</td></tr><tr><td>MMD (MMD(train)=0.9 × 10<sup>-2</sup>)</td><td>1.5 × 10<sup>-2</sup></td><td>1.3 × 10<sup>-2</sup></td></tr></tbody></table>

<figure><img src="image_48.png" /><figcaption>Figure 44: Evolution of the solution log-pdf log($p_s(x_1, x_3)$) of SGM forward SDE (top) and backward SDE (bottom) for the vorticity images distribution, with nominal scheduling ($\beta_m = 0.025, \beta_M = 5.0$) and default neural network architecture.</figcaption></figure>

fail and MMDs are similarly large. The biggest MMD gap between SGM and MSGM appears to be in the intermediate region: $2^{10} = 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 of time steps, MSGM seems much better than SGM. This is expected since the MSGM latent space is already close to the data distribution. Without enough ADAMS iterations, neither the MSGM nor the SGM samples accurately mimic the data distribution, and in any case, it is better to let the optimization procedure run for a long enough time.

<figure><img src="image_49.png" /><figcaption>Figure 45: Evolution of the solution log-pdf log($p_s(x_1, x_3)$) of SGM forward SDE (top) and backward SDE (bottom) for the vorticity images distribution, with nominal scheduling ($\beta_m = 0.025, \beta_M = 5.0$) and our neural network architecture based on spherical decomposition equation L.33.</figcaption></figure>
<figure><img src="image_50.png" /><figcaption>Figure 46: Convergence of MMD (mean and 80% confidence interval) for the vorticity images distribution as a function of number of training data for (from left to right and from top to bottom) dimension *d* = 2, 4, 8, and 16.</figcaption></figure>

<figure><img src="image_51.png" /><figcaption>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 a function of dimension (bottom).</figcaption></figure>
<table><caption>Table 8: High-dimensional vorticity test case: parameters of the nominal numerical experiments.</caption><thead><tr><th>Parameter</th><th>SGM</th><th>MSGM</th></tr></thead><tbody><tr><td>Dimension <i>d</i></td><td>1024</td><td>1024</td></tr><tr><td>Number of used training data points (<i>M</i>)</td><td>5000</td><td>5000</td></tr><tr><td>Number of test data points</td><td>2500</td><td>2500</td></tr><tr><td>Reference number of ADAMS steps</td><td>10<sup>5</sup></td><td>10<sup>5</sup></td></tr><tr><td>Number of ADAMS steps (<i>N</i><sub>iter</sub>)</td><td>10<sup>5</sup></td><td>10<sup>5</sup></td></tr><tr><td>GPU time / ADAMS steps (in ms)</td><td>410</td><td>590</td></tr><tr><td>Batch size</td><td>128</td><td>128</td></tr><tr><td>Number of time steps (forward) <i>N</i><sub>T</sub><sup>f</sup></td><td>1</td><td>128</td></tr><tr><td>Number of time steps (backward) <i>N</i><sub>T</sub><sup>b</sup></td><td>2048</td><td>2048</td></tr><tr><td>β<sub>min</sub></td><td>0.8</td><td>0.8</td></tr><tr><td>β<sub>max</sub></td><td>160</td><td>160</td></tr><tr><td><i>t</i><sub>e</sub></td><td>8 × 10<sup>-3</sup></td><td>8 × 10<sup>-3</sup></td></tr><tr><td>Learning rate</td><td>10<sup>-4</sup></td><td>10<sup>-4</sup></td></tr><tr><td>Neural network architecture</td><td>default</td><td>spherical (equation L.33)</td></tr><tr><td>MMD (MMD(train)=1.4 × 10<sup>-3</sup>)</td><td>2.4 × 10<sup>-3</sup></td><td>3.2 × 10<sup>-3</sup></td></tr></tbody></table>

<figure><img src="image_52.png" /><figcaption>Figure 48: Spatial cropping and subsampling (left), spatial smoothing (middle), and spatial subsampling again (right) of a vorticity field to obtain a data sample at lower but still high dimension $d = 1024$.</figcaption></figure>

### M.7.2 High-dimensional test case : vorticity image processing

To demonstrate that MSGM can address high-dimensional problems, we propose here an image generator based on the sparse local tensor of Appendix K.2.2 and the Unet detailed in Appendix L.4.3. From the original high-resolution PIV-based vorticity images of Figure 39, we crop, subsample at resolution 64 × 64, smooth and subsample again images them spatially, keeping 32 × 32 pixels by images as illustrated by Figure 48. Once reshaped as a vector, each small image represents a data point of dimension 1024.

Figures 49 and 50 present generated images with MSGM and SGM respectively. Table 8 summarizes the parameters of our numerical experiment. The numerical evaluation of image generation skills of MSGM is beyond the scope of this paper and we postpone this study to future work.

## 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 own noising process, leading to different invariant distributions, i.e. Gaussian for SGM and rotational invariant for MSGM. Both latent spaces are tractable, allowing for fast initial sample generation for the reverse process. As a 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 Appendix E.6. From the convergence speed, both dynamics allow for exponential
<figure><img src="image_53.png" /><figcaption>Figure 49: 32 × 32 image generation from MSGM with forward (top) and backward diffusion (bottom) at time (from left to right) $s = T - t = 0, 0.25$, and $T = 1$.</figcaption></figure>

<figure><img src="image_54.png" /><figcaption>Figure 50: 32 × 32 image generation from SGM with forward (top) and backward diffusion (bottom) at time (from left to right) $s = T - t = 0, 0.25$, and $T = 1$. The apparent heteroskedasticity in the diffusion is due to the data normalization (pixel-wise variance is larger on top and bottom boundaries).</figcaption></figure>
convergence to the invariant distribution, assuming the rank condition A2 is satisfied for $G$. We will
conclude this section with a 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 Appendix E and Proposition E.5.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/backward process of MSGM serves as a 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 to rely on 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 \to \infty$, 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 Appendix J. This is a limit in terms of scalability due to the $d^3$ scaling of $G$. Here, the sparse tensors discussed in Appendix K will serve as a 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 Appendix K.1.1.

## N.2 Empirical aspects

SGM offers exact integration of the noising process, while MSGM relies on numerical integration. Although this at first glance looks like a drawback in praxis, for most of our test cases, only a few 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, we refer to the *fair comparison* discussion in Appendix M.3. As our current experiments suggest, MSGM requires less data in training. From approximation theory, learning the score reduces to training on a support that is the hyper-sphere in $\mathbb{R}^d$, with a conditioning variable $\log \|x\| \in \mathbb{R}$. In particular, the effective domain for learning a neuronal remains bounded in $d$. It may affect the stability of the approximation using such an approximation class. Finally, the stabilization due to the conservation of norm avoids divergence instabilities of SSM solvers for MSGM, when compared to well known instabilities of SSM solvers for SGM.