SymDiff: Equivariant Diffusion via Stochastic Symmetrisation

## Abstract

Abstract We propose SymDiff , a novel method for constructing equivariant diffusion models using the recently introduced framework of stochastic symmetrisation. SymDiff resembles a learned data augmentation that is deployed at sampling time, and is lightweight, computationally efficient, and easy to implement on top of arbitrary off-the-shelf models.
Notably, in contrast to previous work, SymDiff typically does not require any neural network components that are intrinsically equivariant, avoiding the need for complex parameterizations and the use of higher-order geometric features.
Instead, our method can leverage highly scalable modern architectures as drop-in
replacements for these more constrained alternatives.
We show that this additional flexibility yields significant empirical benefit on E ⁡ ( 3 ) E 3 \operatorname{E}(3) roman_E ( 3 ) -equivariant molecular generation.
To the best of our knowledge, this is the first application of symmetrisation to generative modelling, suggesting its potential in this domain more generally.

## 1 Introduction

For geometrically structured data such as $N$-body systems of molecules or proteins, it is often of interest to obtain a diffusion that is *equivariant* with respect to some group actions (Xu et al., ; Hoogeboom et al., ; Watson et al., ; Yim et al., ).
By accounting for the large number of permutation and rotational symmetries within such systems, equivariant diffusion models aim to improve the data efficiency (Batzner et al., ) and generalisation (Elesedy & Zaidi, ) of the generative model.
In previous work this has been achieved by the use of intrinsically equivariant neural networks, which constrain each of their linear and nonlinear layers individually to be equivariant, so that the overall model is also equivariant as a result (Bronstein et al., ).

There is a growing literature on the practical limitations of intrinsically equivariant neural networks (Wang et al., ; Canez et al., ; Abramson et al., ; Pertigkiozoglou et al., ).
It has been observed that these models can suffer from degraded training dynamics due to their imposition of architectural constraints, as well as increased computational cost and implementation complexity (Duval et al., ).
To address these issues, there has been recent interest in *symmetrisation* techniques for obtaining equivariance instead (Murphy et al., ; Puny et al., ; Duval et al., ; Kim et al., ; Kaba et al., ; Mondal et al., ; Dym et al., ; Gelberg et al., ).
These approaches offer a mechanism for constructing equivariant neural networks using subcomponents that are not equivariant, which previous work has shown leads to better performing models (Yarotsky, ; Kim et al., ).
However, the advantages of symmetrisation-based equivariance have not yet been explored for generative modelling.
This is possibly in part because previous work has solely focused on deterministic equivariance, rather than the more complex condition of *stochastic* equivariance that is required in the context of generative models.

In this paper, we introduce SymDiff, a novel methodology for obtaining equivariant diffusion models through symmetrisation, rather than intrinsic equivariance.
We build on the recent framework of *stochastic symmetrisation* developed by Cornish () using the theory of Markov categories (Fritz, ).
Unlike previous work on symmetrisation, which operates on deterministic functions, stochastic symmetrisation can be applied to *Markov kernels* directly in distribution space.
We show that this leads naturally to a more flexible approach for constructing equivariant diffusion models than is possible using intrinsic architectures.
We apply this concretely to obtain an $\operatorname{E}(3)$-equivariant diffusion architecture for modelling $N$-body systems, where $\operatorname{E}(3)$ denotes the Euclidean group.
For this task, we formulate a novel reverse process that is allowed to be non-Gaussian, for which we derive a tractable optimisation objective.
Our model is stochastically $\operatorname{E}(3)$-equivariant overall without needing any intrinically $\operatorname{E}(3)$-equivariant neural networks as subcomponents.
We also sketch how to extend SymDiff to score and flow-based generative models (Song et al., ; Lipman et al., ).

We implemented SymDiff for de novo molecular generation, using it as a drop-in replacement for the $\operatorname{E}(3)$-equivariant diffusion of Hoogeboom et al. (), which relies on intrinsically equivariant neural networks.
In contrast, our model is able to leverage highly scalable off-the-shelf architectures such as Diffusion Transformers (Peebles & Xie, ) for all of its subcomponents.
We demonstrate this leads to significantly improved empirical performance for both the QM9 and GEOM-Drugs datasets.

## 2 Background

### 2.1 Equivariant Markov Kernels

##### Markov kernels

At a high level, a *Markov kernel* $k:\mathcal{X}\to\mathcal{Y}$ may be thought of as a *conditional distribution* or *stochastic map* that, when given an input $\mathbf{x}\in\mathcal{X}$, produces a random output in $\mathcal{Y}$ with distribution $k(d\mathbf{y}|\mathbf{x})$.
For example, given a function $f:\mathcal{X}\times\mathcal{E}\to\mathcal{Y}$ and a random element $\eta$ of $\mathcal{E}$, there is a Markov kernel $k:\mathcal{X}\to\mathcal{Y}$ for which $k(d\mathbf{y}|\mathbf{x})$ is the distribution of $f(\mathbf{x},\eta)$.(^1^11For “nice” choices of $\mathcal{Y}$, the converse also holds by *noise outsourcing* (Kallenberg, , Lemma 3.22).)
As a special case, every deterministic function $f:\mathcal{X}\to\mathcal{Y}$ may be thought of as a Markov kernel $\mathcal{X}\to\mathcal{Y}$ also.
When $k(d\mathbf{y}|\mathbf{x})$ has a *density* (or *likelihood*), we will denote this by $k(\mathbf{y}|\mathbf{x})$, although we note that we can still reason about Markov kernels even when they do not admit a likelihood in this sense.

##### Stochastic equivariance

Let $\mathcal{G}$ be a group acting on spaces $\mathcal{X}$ and $\mathcal{Y}$.
We will denote this using “dot” notation, so that the action on $\mathcal{X}$ is a function $(g,\mathbf{x})\mapsto g\cdot\mathbf{x}$.
Recall that a function $f:\mathcal{X}\to\mathcal{Y}$ is then *equivariant* if

$$ $f(g\cdot\mathbf{x})=g\cdot f(\mathbf{x})\qquad\text{for all $\mathbf{x}\in \mathcal{X}$ and $g\in\mathcal{G}$}.$ $$

Here $f$ is purely deterministic, and so this concept must be generalised in order to encompass models whose outputs are stochastic (Bloem-Reddy & Teh, ).
To this end, Cornish () uses a notion defined for Markov kernels: a Markov kernel $k:\mathcal{X}\to\mathcal{Y}$ is *stochastically equivariant* if

$$ $k(d\mathbf{y}|g\cdot\mathbf{x})=(g\cdot k)(d\mathbf{y}|\mathbf{x})\qquad\text{ for all $\mathbf{x}\in\mathcal{X}$ and $g\in\mathcal{G}$},$ (1) $$

where the right-hand side denotes the distribution of $g\cdot\mathbf{y}$ when $\mathbf{y}\sim k(d\mathbf{y}|\mathbf{x})$, or in other words the *pushforward* of $k(d\mathbf{y}|\mathbf{x})$ under $g$.
When $k$ is obtained from $f$ and $\eta$ as above, equation 1 holds iff

$$ $f(g\cdot\mathbf{x},\eta)\overset{\mathrm{d}}{=}g\cdot f(\mathbf{x},\eta)\qquad \text{for all $\mathbf{x}\in\mathcal{X}$ and $g\in\mathcal{G}$},$ $$

where $\overset{\mathrm{d}}{=}$ denotes equality in distribution.(^2^22Note this is more general than the *almost sure* condition from equation (17) of Bloem-Reddy & Teh ().)
If $\eta$ is constant, this says $f$ is deterministically equivariant in the usual sense.
Likewise, when $k$ has a conditional density $k(\mathbf{y}|\mathbf{x})$, equation 1 holds if

$$ $k(g\cdot\mathbf{y}|g\cdot\mathbf{x})=k(\mathbf{y}|\mathbf{x})\qquad\text{for all $\mathbf{x}\in\mathcal{X}$, $\mathbf{y}\in\mathcal{Y}$, and $g\in\mathcal{ G}$,}$ $$

provided the action of $\mathcal{G}$ on $\mathcal{Y}$ has unit Jacobian (Cornish, , Proposition 3.1), as will be the case for all the actions we consider.
This latter condition recovers the usual formulation of stochastic equivariance considered in the diffusion literature by e.g. Xu et al. (); Hoogeboom et al. ().

### 2.2 Equivariant diffusion models

##### Denoising diffusion models

Diffusion models construct a generative model $p_{\theta}(\mathbf{z}_{0})$ of data $p_{\text{data}}(\mathbf{z}_{0})$, living in a space $\mathcal{Z}$, by learning to reverse an iterative forward noising process $\mathbf{z}_{t}$.
Following the notation of Kingma et al. (), the distribution of $\mathbf{z}_{t}$ is defined by $q(\mathbf{z}_{t}|\mathbf{z}_{0})=\mathcal{N}(\mathbf{z}_{t};\alpha_{t}\mathbf{
z}_{0},\sigma_{t}^{2}\mathbf{I})$ for some noise schedule $\alpha_{t},\sigma_{t}>0$ such that the signal-to-noise ration $\text{SNR}(t)\coloneqq\alpha_{t}^{2}/\sigma_{t}^{2}$ is strictly monotonically decreasing.
The joint distribution of the forward and reverse processes respectively then have the forms:

$$ $\displaystyle q(\mathbf{z}_{0:T})=q(\mathbf{z}_{0})\prod_{t=1}^{T}q(\mathbf{z} _{t}|\mathbf{z}_{t-1})\qquad p_{\theta}(\mathbf{z}_{0:T})=p(\mathbf{z}_{T}) \prod_{t=1}^{T}p_{\theta}(\mathbf{z}_{t-1}|\mathbf{z}_{t})$ (2) $$

with $q(\mathbf{z}_{0})\coloneqq p_{\text{data}}(\mathbf{z}_{0})$, and $q(\mathbf{z}_{t}|\mathbf{z}_{t-1})=\mathcal{N}(\mathbf{z}_{t};\alpha_{t|t-1}
\mathbf{z}_{t-1},\sigma_{t|t-1}^{2}\mathbf{I})$, with constants defined as $\alpha_{t|t-1}\coloneqq\alpha_{t}/\alpha_{t-1}$ and $\sigma_{t|t-1}^{2}\coloneqq\sigma_{t}^{2}-\alpha_{t|t-1}^{2}\sigma_{t-1}^{2}$. We take $p(\mathbf{z}_{T})$ to also be Gaussian.
The reverse process is then trained to maximise the ELBO of $\log p_{\theta}(\mathbf{z}_{0})$, which can be obtained as follows (Sohl-Dickstein et al., ):

$$ $\log p_{\theta}(\mathbf{z}_{0})\geq\mathbb{E}_{q(\mathbf{z}_{1}|\mathbf{z}_{0} )}[\log p_{\theta}(\mathbf{z}_{0}|\mathbf{z}_{1})]-D_{\mathrm{KL}}(q(\mathbf{z }_{T}|\mathbf{z}_{0})||p(\mathbf{z}_{T}))\\ -\sum_{t=2}^{T}\mathbb{E}_{q(\mathbf{z}_{t}|\mathbf{z}_{0})}[D_{\mathrm{KL}}(q (\mathbf{z}_{t-1}|\mathbf{z}_{t},\mathbf{z}_{0})||p_{\theta}(\mathbf{z}_{t-1}| \mathbf{z}_{t}))].$ (3) $$

This objective can be efficiently optimised when the reverse process is Gaussian.
This is due to the fact that the posterior distributions $q(\mathbf{z}_{t-1}|\mathbf{z}_{t},\mathbf{z}_{0})=\mathcal{N}(\mathbf{z}_{t-1}
;\mu_{q}(\mathbf{z}_{t},\mathbf{z}_{0}),\sigma_{q}^{2}(t)\mathbf{I})$ are Gaussian, and that the KL divergence between Gaussians can be expressed in closed form.
To match these posteriors, the reverse process is then typically defined in terms of a neural network $\mu_{\theta}:\mathcal{Z}\to\mathcal{Z}$ as

$$ $\displaystyle p_{\theta}(\mathbf{z}_{t-1}|\mathbf{z}_{t})\coloneqq\mathcal{N}( \mathbf{z}_{t-1};\mu_{\theta}(\mathbf{z}_{t}),\sigma_{q}^{2}(t)\mathbf{I}).$ (4) $$

##### Invariant and equivariant diffusion

It is often desirable for the model $p_{\theta}(\mathbf{z}_{0})$ to be *invariant* with respect to the action of a group $\mathcal{G}$.
Invariance here means that

$$ $p_{\theta}(g\cdot\mathbf{z}_{0})=p_{\theta}(\mathbf{z}_{0})\qquad\text{for all $g\in\mathcal{G}$.}$ $$

Intuitively, this says that the density $p_{\theta}(\mathbf{z}_{0})$ is constant on the orbits of this action.
For the model in equation 2, invariance follows if $p_{\theta}(\mathbf{z}_{T})$ is itself invariant, and if each $p_{\theta}(\mathbf{z}_{t-1}|\mathbf{z}_{t})$ is stochastically $\mathcal{G}$-equivariant (Xu et al., ).
All previous work we are aware of has approached this problem by obtaining a deterministically equivariant $\mu_{\theta}$, which then implies that the reverse process is stochastically equivariant (Le et al., ).

### 2.3 N-body systems and E ⁢ ( 3 ) 𝐸 3 E(3) italic_E ( 3 ) -equivariant diffusion

Equivariant diffusion models are often applied to model $N$-body systems such as molecules and proteins (Xu et al., ; Hoogeboom et al., ; Yim et al., ).
This is motivated by the large number of symmetries present in these system.
For example, intuitively speaking, neither the coordinate system nor the ordering of bodies in the system should matter for sampling.
We describe the standard components of such models below.

##### N 𝑁 N italic_N -body data

In 3-dimensions, the state of an $N$-body system can be encoded as a pair $\mathbf{z}=[\mathbf{x},\mathbf{h}]\in\mathbb{R}^{N\times(3+d)}$, where $\mathbf{x}\equiv(\mathbf{x}^{(1)},\ldots,\mathbf{x}^{(N)})\in\mathbb{R}^{N
\times 3}$ describes a set of $N$ points in 3D space, and $\mathbf{h}\equiv(\mathbf{h}^{(1)},\ldots,\mathbf{h}^{(N)})\in\mathbb{R}^{N
\times d}$ describes a set of $N$ feature vectors of dimension $d$.
Each feature vector $\mathbf{h}^{(i)}$ is associated with the point $\mathbf{x}^{(i)}$.
For example, Hoogeboom et al. () encodes molecules in this way, where each $\mathbf{x}^{(i)}$ denotes the location of an atom, and $\mathbf{h}^{(i)}$ some corresponding properties such as atom type (represented as continuous quantities).

##### Center of mass free space

Intuitively speaking, for many applications, the location of an $N$-body system in space should not matter.
For this reason, instead of defining a diffusion on the full space of $N$-body systems directly, previous work (Garcia Satorras et al., ; Xu et al., ; Hoogeboom et al., ) has set $\mathcal{Z}\coloneqq\mathcal{U}$, where $\mathcal{U}$ is the “center of mass (CoM) free” linear subspace of $\mathbb{R}^{N\times(3+d)}$ consisting of $[\mathbf{x},\mathbf{h}]$ such that $\bar{\mathbf{x}}\coloneqq\frac{1}{N}\sum_{i=1}^{N}\mathbf{x}^{(i)}=0$.
In this way, samples from their model are always guaranteed to be centered at the origin.

##### CoM-free diffusions

To construct their forward and reverse processes to now live entirely on $\mathcal{U}$ instead of $\mathbb{R}^{N\times(3+d)}$, Xu et al. () defines the *projected Gaussian distribution* $\mathcal{N}_{\mathcal{U}}(\mu,\sigma^{2}\mathbf{I})$, for $\mu\in\mathcal{U}$ and $\sigma^{2}>0$, as the distribution of $\mathbf{z}$ obtained via the following process:

$$ $\epsilon\sim\mathcal{N}(0,\mathbf{I})\qquad\mathbf{z}\coloneqq\mu+\sigma\,\pi( \epsilon),$ $$

where $\pi:\mathbb{R}^{N\times(3+d)}\to\mathcal{U}$ centers its input at the origin, so that $\pi([\mathbf{x},\mathbf{h}])\coloneqq[\mathbf{x}-(\bar{\mathbf{x}},\ldots,\bar
{\mathbf{x}}),\mathbf{h}]$.
By construction, the projected Gaussian distribution is then supported on the linear subspace $\mathcal{U}$.
Xu et al. () shows that this distribution has a density with Gaussian form $\mathcal{N}_{\mathcal{U}}(\mathbf{z};\mu,\sigma^{2}\mathbf{I})\propto\mathcal{
N}(\mathbf{z};\mu,\sigma^{2}\mathbf{I})$ defined for points $\mathbf{z}\in\mathcal{U}$ living in this subspace.
The forward and reverse processes $q(\mathbf{z}_{t-1}|\mathbf{z}_{t})$ and $p_{\theta}(\mathbf{z}_{t-1}|\mathbf{z}_{t})$ in equation 2 can then be defined exactly as before, but now as projected Gaussians.
Since these are still Gaussian (albeit restricted to a linear subspace), the KL terms in equation 3 remain tractable, which permits optimising the ELBO in the usual way.
This does require that $\mu_{\theta}:\mathcal{U}\to\mathcal{U}$ is now constrained to take values in the subspace $\mathcal{U}$.
Prior work has enforced this by taking $\mu_{\theta}\coloneqq\pi\circ\varphi_{\theta}$, where $\varphi_{\theta}:\mathcal{U}\to\mathbb{R}^{N\times(3+d)}$ is a neural network.

##### Invariance and equivariance

Intuitively, the ordering of the $N$ points and the orientation of the overall system in 3D space should not matter.
To formalise this, let $S_{N}$ denote the symmetric group of permutations of the integers $\{1,\ldots,N\}$, and $\operatorname{O}(3)$ denote the group of orthogonal $3\times 3$ matrices.
Their product $S_{N}\times\operatorname{O}(3)$ acts on $N$-body systems by reordering and orthogonally transforming points as follows:

$$ $(\sigma,R)\cdot[\mathbf{x},\mathbf{h}]\coloneqq[R\mathbf{x}^{(\sigma(1))}, \ldots,R\mathbf{x}^{(\sigma(N))},\mathbf{h}^{(\sigma(1))},\ldots,\mathbf{h}^{( \sigma(N))}],$ (5) $$

where $\sigma\in S_{N}$ and $R\in\operatorname{O}(3)$.
Previous work (Xu et al., ; Hoogeboom et al., ) has then chosen the model $p_{\theta}(\mathbf{z}_{0})$ to be invariant to this action.
They enforce this via the approach described in Section 2.2, by ensuring that $\mu_{\theta}$ is deterministically ($S_{N}\times\operatorname{O}(3)$)-equivariant, which implies that the reverse process is stochastically equivariant also.

Following standard terminology in the literature (Garcia Satorras et al., ; Xu et al., ; Hoogeboom et al., ), we refer to a ($S_{N}\times\operatorname{O}(3)$)-diffusion defined on the CoM-free space $\mathcal{U}$ as an *$\operatorname{E}(3)$-equivariant diffusion*.

## 3 Equivariant Diffusion via Stochastic Symmetrisation

In this section, we introduce our methodology SymDiff, and apply this concretely to the problem of obtaining $\operatorname{E}(3)$-equivariant diffusion models for $N$-body systems.
We also discuss extensions to score and flow-based generative models (Song et al., ; Lipman et al., ) in Appendix D.

### 3.1 Stochastic symmetrisation

Recently, Cornish () gave a general theory of *symmetrisation* using the framework of *Markov categories* (Fritz, ).
This work provides a characterisation of all possible symmetrisation procedures, recovering previous examples in the deterministic setting as special cases (Murphy et al., ; Puny et al., ; Kaba et al., ; Kim et al., ).
The theory applies flexibly and compositionally to general groups and actions, including in the non-compact case, and extends beyond deterministic functions to provide a methodology for symmetrising *Markov kernels*, which had not previously been considered.

In this work, we will consider a special case of the symmetrisation framework of Cornish ().
Concretely, for groups $\mathcal{H}$ and $\mathcal{G}$, we will consider the problem of symmetrising an $\mathcal{H}$-equivariant Markov kernel so that it becomes *both* $\mathcal{H}$-equivariant *and* $\mathcal{G}$-equivariant.
Intuitively, this allows “upgrading” a “partially equivariant” Markov kernel to become “more equivariant”.
By taking $\mathcal{H}$ to be the trivial group consisting of just the identity element, this in particular allows converting arbitrary unconstrained Markov kernels to ones that are $\mathcal{G}$-equivariant.
However, as we show later, it is often more computationally efficient to symmetrise a Markov kernel that is already “partially equivariant” as opposed to fully unconstrained, which motivates studying the case of general $\mathcal{H}$ here.

##### General procedure

Suppose the product group $\mathcal{H}\times\mathcal{G}$ acts on both the spaces $\mathcal{X}$ and $\mathcal{Y}$.
This means we have an action $(h,g)\cdot\mathbf{x}$ defined for $(h,g)\in\mathcal{H}\times\mathcal{G}$ and $\mathbf{x}\in\mathcal{X}$, and similarly for $\mathcal{Y}$ (see equation 5 below for a concrete example involving $N$-body systems).
This implies that both $\mathcal{H}$ and $\mathcal{G}$ act on both $\mathcal{X}$ and $\mathcal{Y}$ individually: for example, $\mathcal{H}$ acts on $\mathcal{X}$ via $h\cdot\mathbf{x}\coloneqq(h,e_{\mathcal{G}})\cdot\mathbf{x}$, where $e_{\mathcal{G}}\in\mathcal{G}$ is the identity element. The following result, which is a special case of Example 6.3 of Cornish (), then gives a procedure for symmetrising $\mathcal{H}$-equivariant Markov kernels to become ($\mathcal{H}\times\mathcal{G}$)-equivariant ones.
For the proof, see Appendix A.1.

###### Proposition 1 .

Let $\gamma:\mathcal{X}\to\mathcal{G}$ be a Markov kernel that is $\mathcal{H}$-invariant (so that $\gamma(dg|h\cdot\mathbf{x})=\gamma(dg|\mathbf{x})$) and $\mathcal{G}$-equivariant, where $\mathcal{G}$ acts on itself by left multiplication, i.e. $g\cdot g^{\prime}\coloneqq gg^{\prime}$.
Then every $\mathcal{H}$-equivariant Markov kernel $k:\mathcal{X}\to\mathcal{Y}$ gives rise to an ($\mathcal{H}\times\mathcal{G}$)-equivariant Markov kernel $\operatorname{\mathsf{sym}}_{\gamma}(k):\mathcal{X}\to\mathcal{Y}$, where $\operatorname{\mathsf{sym}}_{\gamma}(k)(d\mathbf{y}|\mathbf{x})$ is the distribution of $\mathbf{y}$ sampled as follows:

$$ $\displaystyle g\sim\gamma(dg|\mathbf{x})\qquad\mathbf{y}_{0}\sim k(d\mathbf{y} |g^{-1}\cdot\mathbf{x})\qquad\mathbf{y}\coloneqq g\cdot\mathbf{y}_{0}.$ (6) $$

Cornish () shows that this is the only natural symmetrisation procedure that can be defined here without further assumptions on the groups, spaces, or Markov kernels involved.

###### Proposition 1 .

##### Recursive symmetrisation

Notice that the symmetrisation procedure defined by Proposition 1 requires $\gamma$ already to satisfy some equivariance constraints.
In effect this pushes back the problem of obtaining ($\mathcal{H}\times\mathcal{G}$)-equivariant Markov kernels to the choice of $\gamma$, which mirrors the situation in the deterministic setting also (Puny et al., ; Kim et al., ).
To work around this, Cornish () proposes a *recursive* strategy for obtaining $\gamma$.
Specifically, the idea is to set

$$ $\gamma\coloneqq\operatorname{\mathsf{sym}}_{\gamma_{0}}(\gamma_{1})$ (7) $$

where $\gamma_{0},\gamma_{1}:\mathcal{X}\to\mathcal{G}$ are Markov kernels, and $\gamma_{0}$ is $\mathcal{H}$-invariant and $\mathcal{G}$-equivariant, but where now $\gamma_{1}$ is only required to be $\mathcal{H}$-invariant, and may behave arbitrarily with respect to $\mathcal{G}$.
For compact $\mathcal{G}$, a suitable choice of $\gamma_{0}$ here is always available as $\gamma_{0}(dg|\mathbf{x})\coloneqq\lambda(dg)$ (Cornish, , Example 6.3), where $\lambda$ denotes the *Haar measure* on $\mathcal{G}$ (Kallenberg, ).(^3^33We could use this $\gamma_{0}$ directly in place of $\gamma$ in Proposition 1, although we found it empirically better to use the more expressive model form in equation 7 instead.)
We note that this use of the Haar measure exploits the stochastic nature of the procedure in Proposition 1 and would not be possible using deterministic symmetrisation methods here instead.

### 3.2 SymDiff : Symmetrised diffusion

We propose to use stochastic symmetrisation to obtain a diffusion process as in Section 2.2 whose reverse kernels are stochastically equivariant.
Specifically, suppose some product group $\mathcal{H}\times\mathcal{G}$ acts on our state space $\mathcal{Z}$.
For each timestep $t\in\{1,\ldots,T\}$, we will choose some $\mathcal{H}$-equivariant Markov kernel $k_{\theta}:\mathcal{Z}\to\mathcal{Z}$ that admits a conditional density $k_{\theta}(\mathbf{z}_{t-1}|\mathbf{z}_{t})$.
Similarly, we will choose some Markov kernel $\gamma_{\theta}:\mathcal{Z}\to\mathcal{G}$ that satisfies the conditions of Proposition 1 when $\mathcal{X}=\mathcal{Z}$.
With these components, we will define our equivariant reverse process to be

$$ $\displaystyle p_{\theta}(\mathbf{z}_{t-1}|\mathbf{z}_{t})\coloneqq \operatorname{\mathsf{sym}}_{\gamma_{\theta}}(k_{\theta})(\mathbf{z}_{t-1}| \mathbf{z}_{t}),$ (8) $$

which is guaranteed to be ($\mathcal{H}\times\mathcal{G}$)-equivariant by Proposition 1.
This defines a conditional *density*, not just a Markov kernel, as a consequence of the next result.
For the proof, see Appendix A.2.

###### Proposition 2 .

Under the setup of Proposition 1, let $k(d\mathbf{y}|g,\mathbf{x})$ be the distribution of $\mathbf{y}$ sampled via

$$ $\mathbf{y}_{0}\sim k(d\mathbf{y}|g^{-1}\cdot\mathbf{x})\qquad\mathbf{y} \coloneqq g\cdot\mathbf{y}_{0}.$ $$

If $k(d\mathbf{y}|\mathbf{x})$ has a density $k(\mathbf{y}|\mathbf{x})$, then so do $k(d\mathbf{y}|g,\mathbf{x})$ and $\operatorname{\mathsf{sym}}_{\gamma}(k)(d\mathbf{y}|\mathbf{x})$, where the latter is

$$ $\operatorname{\mathsf{sym}}_{\gamma}(k)(\mathbf{y}|\mathbf{x})=\mathbb{E}_{ \gamma(dg|\mathbf{x})}[k(\mathbf{y}|g,\mathbf{x})].$ $$

When the action of $\mathcal{G}$ on $\mathcal{Y}$ has unit Jacobian, we also have $k(\mathbf{y}|g,\mathbf{x})=k(g^{-1}\cdot\mathbf{y}|g^{-1}\cdot\mathbf{x})$.

###### Proposition 2 .

##### Training objective

We would like to learn the parameters $\theta$ using the ELBO from equation 3.
However, in general, we do not have access to the densities $p_{\theta}(\mathbf{z}_{t-1}|\mathbf{z}_{t})$ from equation 8 in closed form, since this requires computing the expectation as in Proposition 2.
As such, we cannot compute the ELBO directly.
However, since $\log$ is concave, Jensen’s inequality allows us to bound

$$ $\log p_{\theta}(\mathbf{z}_{t-1}|\mathbf{z}_{t})\geq\mathbb{E}_{\gamma_{\theta }(dg|\mathbf{z}_{t})}[\log k_{\theta}(\mathbf{z}_{t-1}|g,\mathbf{z}_{t})].$ $$

Since the ELBO in equation 3 depends linearly on $\log p_{\theta}(\mathbf{z}_{t-1}|\mathbf{z}_{t})$, this allows us to also bound

$$ $\log p_{\theta}(\mathbf{z}_{0})\geq\overbrace{\mathbb{E}_{q(\mathbf{z}_{1}| \mathbf{z}_{0}),\gamma_{\theta}(dg|\mathbf{z}_{1})}[\log k_{\theta}(\mathbf{z} _{0}|g,\mathbf{z}_{1})]}^{-\mathcal{L}_{0}}-D_{\mathrm{KL}}(q(\mathbf{z}_{T}| \mathbf{z}_{0})||p(\mathbf{z}_{T}))\\ -\sum_{t=2}^{T}\overbrace{\mathbb{E}_{q(\mathbf{z}_{t}|\mathbf{z}_{0}),\gamma_ {\theta}(dg|\mathbf{z}_{t})}[D_{\mathrm{KL}}(q(\mathbf{z}_{t-1}|\mathbf{z}_{t} ,\mathbf{z}_{0})||k_{\theta}(\mathbf{z}_{t-1}|g,\mathbf{z}_{t})]}^{\mathcal{L} _{t}},$ (9) $$

where the right-hand side is a tractable lower bound to the original ELBO in equation 3. We take this new bound as our objective used to train our SymDiff model.

##### Comparison with deterministic symmetrisation

An alternative approach to equation 8 would be to obtain $\mu_{\theta}$ in equation 4 via *deterministic* symmetrisation.
An advantage of our method is that we require only a *single* pass through $\gamma_{\theta}$ and $k_{\theta}$ at sampling time.
In contrast, many deterministic techniques (Puny et al., ; Kim et al., ) involve a Monte Carlo averaging step that requires multiple passes through the model instead, and then only return an approximation, which would introduce sampling bias.
The canonicalisation method of Kaba et al. () avoids this Monte Carlo step, but instead suffers from pathologies associated with its deterministic analogue of $\gamma_{\theta}$, which Dym et al. () show must become discontinuous at certain inputs.
Additionally, canonicalisation requires an intrinsically $\mathcal{G}$-equivariant architecture for its analogue of $\gamma_{\theta}$.
In contrast, whenever $\mathcal{G}$ is compact, our recursive strategy in equation 7 requires no intrinsically $\mathcal{G}$-equivariant neural network components anywhere in the model at all.

### 3.3 SymDiff for N 𝑁 N italic_N -body systems

We now apply SymDiff in the setting of $N$-body systems considered in Section 2.3.
Specifically, we take $\mathcal{Z}\coloneqq\mathcal{U}$, $\mathcal{H}\coloneqq S_{N}$, and $\mathcal{G}\coloneqq\operatorname{O}(3)$, and consider the action on $\mathcal{Z}$ defined in equation 5.
This means that we start with $k_{\theta}(\mathbf{z}_{t-1}|\mathbf{z}_{t})$ in equation 8 that is already equivariant with respect to reorderings of the $N$ bodies, and then symmetrise this to obtain an ($S_{N}\times\operatorname{O}(3)$)-equivariant reverse kernel overall.
We choose to symmetrise in this way because highly scalable $S_{N}$-equivariant kernels based on Transformer architectures can be readily constructed for this purpose (Vaswani et al., ; Lee et al., ; Peebles & Xie, ), whereas intrinsically $\operatorname{O}(3)$-equivariant neural networks have not shown the same degree of scalability to-date (Abramson et al., ).

##### Choice of unsymmetrised kernels

It remains now to choose $k_{\theta}(\mathbf{z}_{t-1}|\mathbf{z}_{t})$.
We now do so in a way that equation 9 will resemble the standard diffusion objective in Ho et al. (), allowing for the scalable training of SymDiff.
Specifically, we take

$$ $k_{\theta}(\mathbf{z}_{t-1}|\mathbf{z}_{t})\coloneqq\mathcal{N}_{\mathcal{U}}( \mathbf{z}_{t-1};\mu_{\theta}(\mathbf{z}_{t}),\sigma_{q}^{2}(t)\mathbf{I}),$ (10) $$

where $\mu_{\theta}:\mathcal{Z}\to\mathcal{Z}$ is an arbitrary $S_{N}$-equivariant neural network, which in turn means $k_{\theta}$ is stochastically $S_{N}$-equivariant.
We highlight that $\mu_{\theta}$ is otherwise unconstrained and can process $[\mathbf{x},\mathbf{h}]$ jointly.
In contrast, previous work using intrinsically ($S_{n}\times\operatorname{O}(3)$)-equivariant components (Satorras et al., ; Thölke & De Fabritiis, ; Hua et al., ) has required complex parameterisations for $\mu_{\theta}$ that handle the $\mathbf{x}$ and $\mathbf{h}$ inputs separately.

##### Form of KL terms

We now show how our model yields a closed-form expression for the $\mathcal{L}_{t}$ terms in equation 9.
Standard arguments show that each $q(\mathbf{z}_{t-1}|\mathbf{z}_{t},\mathbf{z}_{0})=\mathcal{N}_{\mathcal{U}}(
\mathbf{z}_{t};\mu_{q}(\mathbf{z}_{t},\mathbf{z}_{0}),\sigma_{q}^{2}(t)\mathbf
{I})$ is a (projected) Gaussian.
We claim that our model also gives a (projected) Gaussian

$$ $k_{\theta}(\mathbf{z}_{t-1}|R,\mathbf{z}_{t})=\mathcal{N}_{\mathcal{U}}( \mathbf{z}_{t-1};R\cdot\mu_{\theta}(R^{T}\cdot\mathbf{z}_{t}),\sigma_{q}^{2}(t )\mathbf{I}).$ (11) $$

Indeed, by the definition of this kernel in Proposition 2 and the definition of $\mathcal{N}_{\mathcal{U}}$ in Section 2.3, and since $R^{-1}=R^{T}$ for $R\in\operatorname{O}(3)$, for $\mathbf{z}_{t-1}\sim k_{\theta}(\mathbf{z}_{t-1}|R,\mathbf{z}_{t})$ we have

$$ $\displaystyle\mathbf{z}_{t-1}$ $\displaystyle\overset{\mathrm{d}}{=}R\cdot\left(\mu_{\theta}(R^{T}\cdot\mathbf {z}_{t})+\sigma_{q}(t)\,\epsilon\right)=R\cdot\mu_{\theta}(R^{T}\cdot\mathbf{z }_{t})+\sigma_{q}(t)R\cdot\epsilon,$ $$

where $\epsilon\sim\mathcal{N}_{\mathcal{U}}(0,\mathbf{I})$.
Since $R\cdot\epsilon\overset{\mathrm{d}}{=}\epsilon$ for $R\in\operatorname{O}(3)$, equation 11 now follows.
The same argument as Hoogeboom et al. () now yields the closed-form expression

$$ $\displaystyle\mathcal{L}_{t}=\mathbb{E}_{q(\mathbf{z}_{t}|\mathbf{z}_{0}), \gamma_{\theta}(dR|\mathbf{z}_{t})}\left[\frac{1}{2\sigma_{q}^{2}(t)}\norm{\mu _{q}(\mathbf{z}_{t},\mathbf{z}_{0})-R\cdot\mu_{\theta}(R^{T}\cdot\mathbf{z}_{t })}^{2}\right].$ (12) $$

We can obtain unbiased gradients of this quantity whenever $\gamma_{\theta}(dR|\mathbf{z}_{t})$ is *reparametrisable* (Kingma, ).
In other words, we should define $\gamma_{\theta}(dR|\mathbf{z}_{t})$ to be the distribution of $\varphi_{\theta}(\mathbf{z}_{t},\xi)$, where $\varphi_{\theta}$ is a deterministic neural network, and $\xi$ is some noise variable whose distribution does not depend on $\theta$. For a discussion of how we can handle the $\mathcal{L}_{0}$ term in equation 9 in our framework, we refer to Appendix C.1.

##### ϵ italic-ϵ \epsilon italic_ϵ -parameterisation

When $\mu_{\theta}$ is taken to have the $\epsilon$-form (Ho et al., ; Kingma et al., )

$$ $\mu_{\theta}(\mathbf{z}_{t})\coloneqq\frac{1}{\alpha_{t|t-1}}\mathbf{z}_{t}- \frac{\sigma_{t|t-1}^{2}}{\alpha_{t|t-1}\sigma_{t}}\epsilon_{\theta}(\mathbf{z }_{t}),$ (13) $$

for some neural network $\epsilon_{\theta}:\mathcal{Z}\to\mathcal{Z}$, the same argument given by Ho et al. () now allows us to rewrite equation 12 in a way that resembles the standard diffusion objective:

$$ $\displaystyle\mathcal{L}_{t}=\mathbb{E}_{q(\mathbf{z}_{0}),\epsilon\sim \mathcal{N}_{\mathcal{U}}(0,\mathbf{I}),\gamma_{\theta}(dR|\mathbf{z}_{t})} \left[\frac{1}{2}w(t)\norm{\epsilon-R\cdot\epsilon_{\theta}(R^{T}\cdot\mathbf{ z}_{t})}^{2}\right]$ (14) $$

where $\mathbf{z}_{t}=\alpha_{t}\mathbf{z}_{0}+\sigma_{t}\epsilon$, and $w(t)=(1-\text{SNR}(t-1)/\text{SNR}(t))$.
Recall we require $\mu_{\theta}$ to be $S_{N}$-equivariant, which straightforwardly follows here whenever $\epsilon_{\theta}$ is. In practice, we set $w(t)=1$ during training as is commonly done in the diffusion literature (Kingma & Gao, ).

##### Recursive choice of γ θ subscript 𝛾 𝜃 \gamma_{\theta} italic_γ start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT

Recall that to apply Proposition 1, we require $\gamma_{\theta}$ to be $S_{N}$-invariant and $\operatorname{O}(3)$-equivariant.
To obtain this, we apply the recursive procedure from equation 7, where $\gamma_{0}:\mathcal{Z}\to\operatorname{O}(3)$ is obtained using the Haar measure on $\operatorname{O}(3)$ as described there, and $\gamma_{1}(dR|\mathbf{z}_{t})$ is defined as the distribution of $f_{\theta}(\mathbf{z}_{t},\eta)$, where $f_{\theta}$ is some $S_{N}$-invariant neural network, and $\eta$ is sampled from some noise distribution $\nu(d\eta)$.
Since the Haar measure does not depend on $\theta$, the overall $\gamma_{\theta}$ obtained in this way remains reparametrisable.
We emphasise that $f_{\theta}$ is *not* required to be $\operatorname{O}(3)$-equivariant, thus allowing for highly flexible choices such as Set Transformers (Lee et al., ).
At sampling time, we use the procedure from Section 5 of Mezzadri () to sample from the Haar measure on $\operatorname{O}(3)$, which is a negligible overhead compared with the cost of evaluating $f_{\theta}$.

Figure: Algorithm 1 SymDiff training step

Figure: Algorithm 2 SymDiff sampling process

### 3.4 Data Augmentation is a Special Case of SymDiff

Data augmentation is a popular method for incorporating “soft” inductive biases within neural networks.
Consider again our setup from the previous section, but suppose we now used the unsymmetrised kernels $k_{\theta}(\mathbf{z}_{t-1}|\mathbf{z}_{t})$ from equation 10 in place of $p_{\theta}(\mathbf{z}_{t-1}|\mathbf{z}_{t})$ in our backwards process directly.
Suppose moreover that we trained our model using the standard diffusion objective, applying a uniform random orthogonal transformation to the input of our $\epsilon_{\theta}$ network before each forward pass.
This is equivalent to optimising the following objective:

$$ $\displaystyle\mathcal{L}_{t}^{\mathrm{aug}}=\mathbb{E}_{q(\mathbf{z}_{0}), \epsilon\sim\mathcal{N}_{\mathcal{U}}(0,\mathbf{I}),\lambda(dR)}\left[\frac{1} {2}w(t)\norm{\epsilon-\epsilon_{\theta}(\alpha_{t}R\cdot\mathbf{z}_{0}+\sigma_ {t}\epsilon)}^{2}\right],$ (15) $$

where $\lambda$ is the Haar measure on $\operatorname{O}(3)$. (More general choices of $\lambda$ could also be considered.)
We then have the following result, proven in Appendix A.3.

###### Proposition 3 .

When $\gamma_{\theta}(dR|\mathbf{z}_{t})=\lambda(dR)$ for all $\mathbf{z}_{t}\in\mathcal{Z}$, our SymDiff objective recovers the data augmentation objective exactly, so that $\mathcal{L}_{t}=\mathcal{L}_{t}^{\mathrm{aug}}$.

In this way, SymDiff may be understood as a strict generalisation of data augmentation in which a more flexible augmentation process $\gamma_{\theta}$ is *learned* during training.
Moreover, our $\gamma_{\theta}$ is then deployed at sampling time, which guarantees stochastic equivariance.
In contrast, data augmentation is usually only applied during training, and the model deployed at sampling time then becomes only approximately equivariant.

###### Proposition 3 .

## 4 Experiments

We evaluated our $E(3)$-equivariant SymDiff model as a drop-in replacement for the $E(3)$-equivariant diffusion (EDM) of Hoogeboom et al. () on both the QM9 and GEOM-Drugs datasets for molecular generation.
We implemented this by taking the code base of Hoogeboom et al. () and substituting our symmetrised reverse process for their one, which relies on intrinsically equivariant neural networks.
We made minimal other changes to their codebase and experimental setup otherwise.
We found SymDiff led to significantly improved performance on both tasks.
We additionally compared against several other more sophisticated models than Hoogeboom et al. () that model more complex molecular structure.
Despite the relative simplicity of our approach, we obtained results that were on par or superior to these additional methods also.
We give an overview of these results now and provide full details in Appendix C.

### 4.1 QM9

##### Dataset

QM9 (Ramakrishnan et al., ) is a common benchmark dataset used for evaluating molecular generation.
It consists of molecular properties and atom coordinates for 130k small molecules with up to 9 heavy atoms and a total of 29 atoms including hydrogen.
For our experiments, we trained our SymDiff method to generate molecules with 3D coordinates, atom types, and atom charges where we explicitly modeled hydrogen atoms.
We used the same train-val-test split of 100K-8K-13K as in Anderson et al. ().

##### Our model

We took our model to be $\operatorname{\mathsf{sym}}_{\gamma_{\theta}}(k_{\theta})(\mathbf{z}_{t-1}|
\mathbf{z}_{t})$ from equation 8 above.
We chose the “backbone” neural network $\epsilon_{\theta}$ of the component $k_{\theta}$ to be the smallest Diffusion Transformer (DiT) model considered by Peebles & Xie (), with 29M parameters.
Likewise, we chose the “backbone” neural network $f_{\theta}$ of the component $\gamma_{\theta}$ to be a DiT with 2.2M parameters.
In this way, our $k_{\theta}$ was much larger than $\gamma_{\theta}$, following a similar approach taken by earlier work on deterministic symmetrisation (Kim et al., ; Kaba et al., ).
By construction, DiT architectures are $S_{N}$-equivariant, and can be made $S_{N}$-invariant using a Set Transformer (Lee et al., ) approach, which allowed us to achieve the requirements for these components from Section 3.3.
However, DiTs are not required to be $\operatorname{O}(3)$-equivariant, unlike the backbone E-GNN components used by Hoogeboom et al. ().
Overall, our model had 31.2M parameters in total.

##### Metrics

To measure the quality of generated molecules, we follow standard practise (Hoogeboom et al., ; Garcia Satorras et al., ) and report atom stability, molecular stability, validity and uniqueness. We exclude results for the novelty metric for the same reasons as discussed in Vignac & Frossard () and refer the reader to these works for a more extensive discussion of these metrics. For all of our metrics, we used 10,000 samples and report the mean and standard deviation over three evaluation runs. To demonstrate the efficiency of our approach, we also report the number of iterations per second and time taken to generate one sample.

##### Baselines

We compared our method with EDM (Hoogeboom et al., ), whose framework is the most similar to ours.
We also considered several other baselines that baked in more sophisticated inductive biases related to molecular generation. These methods included END (Cornet et al., ), GeoLDM (Xu et al., ) and MUDiff (Hua et al., ).
For EDM, we trained their 5.3M parameter model using their experimental setup and for the remaining methods we state the results reported in the papers, noting that all the papers used the same train-val-test splits.
We conducted an ablation study to make sure that the comparison between our 31.2M parameter model and the 5.3M parameter EDM was fair (see Appendix C.2.3).
In addition, we trained an unsymmetrised reverse process with a 29M parameter DiT backbone (DiT), as well as the same model using data augmentation as described in Section 3.4 (DiT-Aug).
We also compared our SymDiff model to the case where we simply let $\gamma_{\theta}$ be obtained using the Haar measure on $\operatorname{O}(3)$, rather than learning this component (SymDiff-H).

##### Scalability

To demonstrate the scalability of our approach we also trained the SymDiff and EDM models for different sizes. For EDM we used sizes of 5.3M (EDM^+, original model size), 9.5M (EDM^++) and 12.4M (EDM^+++) parameters. To make the comparisons fair, we scaled the size of our original model (31.2M, SymDiff^+++) down to 23.5M (SymDiff^++) and 13.5M (SymDiff^+) parameters, where the ratio of the sizes to our largest model is the same as for the EDM models. As a further test of scalability, we also scaled our method to 117.8M (SymDiff^∗) parameters (DiT-B from Peebles & Xie ()). For additional details see Appendix C.2.1.

##### Results

From Table 1, we see that our SymDiff model comfortably beats its EDM counterpart on all metrics, bar uniqueness. Additionally, we see that our model is also competitive with the more recent, sophisticated baselines from the literature, outperforming all of them on validity.
From Table 3 we also see that our method scales significantly better than the EDM approach, both in terms of training and sampling time.
This is intuitive as the intrinsically equivariant E-GNN used by the EDM model uses message passing for training and inference, which is computationally very costly.
In contrast, our symmetrisation approach allows us to use computationally efficient DiT components that parallelise and scale much more effectively (Fei et al., ).

**Table 1: Test NLL, atom stability, molecular stability, validity and uniqueness on QM9 for 10,000 samples and 3 evaluation runs. We omit the results for NLL where not available.**
| Method | NLL $\downarrow$ | Atm. stability (%) $\uparrow$ | Mol. stability (%) $\uparrow$ | Val. (%) $\uparrow$ | Uniq. (%) $\uparrow$ |
| --- | --- | --- | --- | --- | --- |
| GeoLDM | – | 98.90 ${\scriptstyle\pm{0.10}}$ | 89.40 ${\scriptstyle\pm{0.50}}$ | 93.80 ${\scriptstyle\pm{0.40}}$ | 92.70 ${\scriptstyle\pm{0.50}}$ |
| MUDiff | -135.50${\scriptstyle\pm{2.10}}$ | 98.80${\scriptstyle\pm{0.20}}$ | 89.90${\scriptstyle\pm{1.10}}$ | 95.30${\scriptstyle\pm{1.50}}$ | 99.10${\scriptstyle\pm{0.50}}$ |
| END | – | 98.90 ${\scriptstyle\pm{0.00}}$ | 89.10 ${\scriptstyle\pm{0.10}}$ | 94.80 ${\scriptstyle\pm{0.10}}$ | 92.60 ${\scriptstyle\pm{0.20}}$ |
| EDM | -110.70${\scriptstyle\pm{1.50}}$ | 98.70${\scriptstyle\pm{0.10}}$ | 82.00 ${\scriptstyle\pm{0.40}}$ | 91.90${\scriptstyle\pm{0.50}}$ | 90.70${\scriptstyle\pm{0.60}}$ |
| SymDiff^∗ | -133.79${\scriptstyle\pm{1.33}}$ | 98.92${\scriptstyle\pm{0.03}}$ | 89.65${\scriptstyle\pm{0.10}}$ | 96.36${\scriptstyle\pm{0.27}}$ | 97.66${\scriptstyle\pm{0.22}}$ |
| SymDiff | -129.35 ${\scriptstyle\pm{1.07}}$ | 98.74${\scriptstyle\pm{0.03}}$ | 87.49 ${\scriptstyle\pm{0.23}}$ | 95.75 ${\scriptstyle\pm{0.10}}$ | 97.89 ${\scriptstyle\pm{0.26}}$ |
| SymDiff-H | -126.53${\scriptstyle\pm{0.90}}$ | 98.57${\scriptstyle\pm{0.07}}$ | 85.51${\scriptstyle\pm{0.18}}$ | 95.22${\scriptstyle\pm{0.18}}$ | 97.98${\scriptstyle\pm{0.09}}$ |
| DiT-Aug | -126.81${\scriptstyle\pm{1.69}}$ | 98.64${\scriptstyle\pm{0.03}}$ | 85.85${\scriptstyle\pm{0.24}}$ | 95.10${\scriptstyle\pm{0.17}}$ | 97.98${\scriptstyle\pm{0.08}}$ |
| DiT | -127.78${\scriptstyle\pm{2.49}}$ | 98.23${\scriptstyle\pm{0.04}}$ | 81.03${\scriptstyle\pm{0.25}}$ | 94.71${\scriptstyle\pm{0.31}}$ | 97.98${\scriptstyle\pm{0.12}}$ |
| Data |  | 99.00 | 95.20 | 97.8 | 100 |

We attribute the scalability of our method to the extra architectural flexibility provided by our approach to symmetrisation, enabling it to not only outperform its EDM counterpart but also perform competitively with more sophisticated alternatives (see Tables 1, 2).
In particular, we found that our largest model, SymDiff^∗, outperformed all our baselines on atom stability and validity, and is within variance for molecular stability.
By training the model for 500 more epochs, we further found that it surpassed all our baselines on all metrics except uniqueness (see Appendix C.2.4). We conjecture that similar performance improvements can be achieved by using our SymDiff approach as a drop-in replacement for the reverse kernels in more sophisticated methods.

**Table 2: Test NLL, atom stability, molecular stability, validity and uniqueness on QM9 for 10,000 samples and 3 evaluation runs. We omit the results for NLL where not available.**
| Method | NLL $\downarrow$ | Atm. stability (%) $\uparrow$ | Mol. stability (%) $\uparrow$ | Val. (%) $\uparrow$ | Uniq. (%) $\uparrow$ |
| --- | --- | --- | --- | --- | --- |
| EDM^+++ | -119.12${\scriptstyle\pm{1.41}}$ | 98.65${\scriptstyle\pm{0.11}}$ | 85.68${\scriptstyle\pm{0.83}}$ | 93.55${\scriptstyle\pm{0.57}}$ | 98.44${\scriptstyle\pm{0.18}}$ |
| EDM^++ | -110.97${\scriptstyle\pm{1.42}}$ | 98.59${\scriptstyle\pm{0.02}}$ | 84.63${\scriptstyle\pm{0.16}}$ | 93.06${\scriptstyle\pm{0.18}}$ | 98.37${\scriptstyle\pm{0.10}}$ |
| EDM^+ | -110.70${\scriptstyle\pm{1.50}}$ | 98.70${\scriptstyle\pm{0.10}}$ | 82.00 ${\scriptstyle\pm{0.40}}$ | 91.90${\scriptstyle\pm{0.50}}$ | 90.70${\scriptstyle\pm{0.60}}$ |
| SymDiff^+++ | -129.35${\scriptstyle\pm{1.07}}$ | 98.74${\scriptstyle\pm{0.03}}$ | 87.49${\scriptstyle\pm{0.23}}$ | 95.75 ${\scriptstyle\pm{0.10}}$ | 97.89 ${\scriptstyle\pm{0.26}}$ |
| SymDiff^++ | -125.40${\scriptstyle\pm{0.63}}$ | 98.47${\scriptstyle\pm{0.01}}$ | 83.51${\scriptstyle\pm{0.24}}$ | 95.16${\scriptstyle\pm{0.16}}$ | 98.16${\scriptstyle\pm{0.18}}$ |
| SymDiff^+ | -110.68 ${\scriptstyle\pm{2.55}}$ | 97.46${\scriptstyle\pm{0.07}}$ | 71.25 ${\scriptstyle\pm{0.50}}$ | 91.87 ${\scriptstyle\pm{0.27}}$ | 98.82 ${\scriptstyle\pm{0.08}}$ |
| Data |  | 99.00 | 95.20 | 97.8 | 100 |

**Table 3: Iterations per second, and sampling time for different sizes of SymDiff and EDM. Iterations per and sampling time were found by training both models with a batch size of 256.**
| Method | # Parameters | Iters./sec | Sampling time (s) |
| --- | --- | --- | --- |
| EDM^+++ | 12.4M | 2.84 | 0.56 |
| EDM^++ | 9.4M | 3.31 | 0.46 |
| EDM^+ | 5.4M | 5.52 | 0.27 |
| SymDiff^+++ | 31.2M | 14.32 | 0.09 |
| SymDiff^++ | 23.5M | 15.76 | 0.08 |
| SymDiff^+ | 13.5M | 18.98 | 0.07 |

From Table 1, we note that DiT-Aug performs notably better than the DiT model on all metrics, highlighting its strength as a baseline.
Moreover, we observe that our SymDiff model outperforms both SymDiff-H and DiT-Aug on all metrics apart from uniqueness.
This shows the benefit of our approach as a more effective way of incorporating $\operatorname{O}(3)$ equivariance than data augmentation.

### 4.2 GEOM-Drugs

##### Dataset, model and training

GEOM-Drugs (Axelrod & Gomez-Bombarelli, ) is a larger and more complicated dataset than QM9, containing 430,000 molecules with up to 181 atoms.
We processed the dataset in the same way as Hoogeboom et al. (), where we again model hydrogen explicitly.
For our model, we used the same hyperparameters as for QM9 and, given the size of the dataset, a smaller model size of 6.4M parameters with $f_{\theta}$ containing 1.1M parameters.
We trained our model for 40 epochs to match the same number of steps as used in Hoogeboom et al. ().

##### Metrics and baselines

We report the same metrics as for QM9 but exclude molecular stability and uniqueness for the same reasons discussed in Hoogeboom et al. (). We compared our method to the same baselines as for QM9, excluding the methods that did not evaluate on GEOM-Drugs. Again, all our baselines used the same train-test-val splits.
We state the reported results from the original papers rather than retraining these models ourselves.

##### Results

From Table 4 we again see that our approach comfortably outperforms its EDM counterpart. It is also again competitive with the more sophisticated baselines. We also observe desirable scaling behavior like with QM9 (see Appendix C.3.3); in particular we note that our smallest model took only 20h to train and outperformed the EDM’s atom stability by 2.85%, whereas the EDM model resulted in out-of-memory errors on our NVIDIA H100 80GB GPUs.

**Table 4: Test NLL, atom stability and validity on GEOM-Drugs for 10,000 samples and 3 evaluation runs. GeoLDM and EDM ran their results for just one evaluation run. We omit the results for NLL and validity where not available.**
| Method | NLL $\downarrow$ | Atm. stability (%) $\uparrow$ | Val. (%) $\uparrow$ |
| --- | --- | --- | --- |
| GeoLDM | – | 84.4 | 99.3 |
| END | – | 87.8${\scriptstyle\pm{0.99}}$ | 92.9${\scriptstyle\pm{0.3}}$ |
| EDM | -137.1 | 81.3 | – |
| SymDiff | -196.77${\scriptstyle\pm{0.03}}$ | 84.15${\scriptstyle\pm{0.07}}$ | 95.99${\scriptstyle\pm{0.25}}$ |
| Data |  | 86.50 | 99.9 |

We also repeated comparisons with data augmentation and found similar results to QM9 (see Appendix C.3.4).

## 5 Conclusion

We have introduced SymDiff: a novel, lightweight, and scalable framework for constructing equivariant diffusion models based on stochastic symmetrisation.
We applied this approach to $\operatorname{E}(3)$-equivariance for $N$-body data, obtaining an overall model that is stochastically equivariant but that does not rely on any intrinsically equivariant neural network subcomponents.
Our approach leads to significantly greater modelling flexibility, which allows leveraging powerful off-the-shelf architectures such as Transformers (Vaswani et al., ).
We showed empirically that this leads overall to improved performance on several relevant benchmarks.

#### Acknowledgments

The authors are grateful to Tom Rainforth, Emile Mathieu, Saifuddin Syed and Ahmed Elhag for helpful discussions. LZ and KA are supported by the EPSRC CDT in Modern Statistics and Statistical Machine Learning (EP/S023151/1).

## Appendix A Proofs

### A.1 Proof of Proposition 1

###### Proof.

This follows by appropriately instantiating Example 6.3 of Cornish (), whose definitions and notation we will import freely here.
We will also make use of *string diagrams* (Selinger, ), an introduction to which can be found in Section 2 of Cornish ().
Intuitively, these diagrams represent (possibly stochastic) computational processes that should be read up the page, with the inputs applied at the bottom, and outputs produced at the top.

First, we will set $\rho$ in Example 6.3 of Cornish () to be the trivial action of $H$ on $N$, so that the semidirect product $N\rtimes_{\rho}H$ becomes simply the product $N\times H$.
As per Definition 3.7 of Cornish (), the inversion operation for $N\times H$ can then be written as

$$ $=$ $(-)^{-1}_{H}$ $(-)^{-1}_{N}$ $N$ $H$ $H$ $N$ $(-)^{-1}$ $H$ $N$ $H$ $N$ $$

Similarly, by Remark 3.4 of Cornish (), we can write the action $\alpha_{X}$ on $X$ in the form

$$ $=$ $(-)^{-1}_{H}$ $(-)^{-1}_{N}$ $N$ $H$ $H$ $N$ $(-)^{-1}$ $H$ $N$ $H$ $N$ $\alpha_{X,H}$ $\alpha_{X,N}$ $N$ $H$ $X$ $X$ $=$ $\alpha_{X}$ $N$ $H$ $X$ $X$ (16) $$

for some action $\alpha_{X,H}$ of $H$ on $X$, and some action $\alpha_{X,N}$ of $N$ and $X$.
We can do similarly for $\alpha_{Y}$.
In addition, as per Example 3.23 of Cornish (), the inclusion homomorphism $i_{H}$, which plays the role of the $i_{N}$-coset map in Example 6.3, can be written as

$$ $=$ $(-)^{-1}_{H}$ $(-)^{-1}_{N}$ $N$ $H$ $H$ $N$ $(-)^{-1}$ $H$ $N$ $H$ $N$ $\alpha_{X,H}$ $\alpha_{X,N}$ $N$ $H$ $X$ $X$ $=$ $\alpha_{X}$ $N$ $H$ $X$ $X$ $=$ $i_{H}$ $N$ $H$ $N$ $N$ $H$ $e_{N}$ $H$ $$

Finally, the coset action $\ast/H$ of $N\times H$ on $H$ is just

$$ $=$ $(-)^{-1}_{H}$ $(-)^{-1}_{N}$ $N$ $H$ $H$ $N$ $(-)^{-1}$ $H$ $N$ $H$ $N$ $\alpha_{X,H}$ $\alpha_{X,N}$ $N$ $H$ $X$ $X$ $=$ $\alpha_{X}$ $N$ $H$ $X$ $X$ $=$ $i_{H}$ $N$ $H$ $N$ $N$ $H$ $e_{N}$ $H$ $\ast_{H}$ $N$ $H$ $H$ $H$ $=$ $\ast/N$ $N$ $H$ $H$ $H$ $$

where $\ast_{H}$ denotes the group multiplication operation for $H$.
Notice that, like in equation 16, the right-hand side here is decomposed into a pair of an $N$-action (namely the trivial action) and an $H$-action (namely $\ast_{H}$).
By Remark 3.4 of Cornish (), it therefore follows that $\gamma:X\to H$ is equivariant with respect to $\alpha_{X}$ and $\ast/H$ if and only if it is

- •
$N$-equivariant with respect to $\alpha_{X,N}$ and the trivial action on $H$, and
- •
$H$-equivariant with respect to $\alpha_{X,H}$ and the action $\ast_{H}$

In other words, $\gamma$ must be $N$-invariant and $H$-equivariant.

Now suppose we have such a $\gamma$.
Substituting all the components above into the end-to-end procedure from Section 5.4 of Cornish (), we obtain $\operatorname{\mathsf{sym}}_{\gamma}(k):X\to Y$ as follows:

$$ $=$ $(-)^{-1}_{H}$ $(-)^{-1}_{N}$ $N$ $H$ $H$ $N$ $(-)^{-1}$ $H$ $N$ $H$ $N$ $\alpha_{X,H}$ $\alpha_{X,N}$ $N$ $H$ $X$ $X$ $=$ $\alpha_{X}$ $N$ $H$ $X$ $X$ $=$ $i_{H}$ $N$ $H$ $N$ $N$ $H$ $e_{N}$ $H$ $\ast_{H}$ $N$ $H$ $H$ $H$ $=$ $\ast/N$ $N$ $H$ $H$ $H$ $\alpha_{X,H}$ $(-)^{-1}_{H}$ $k$ $\alpha_{Y,H}$ $X$ $\gamma$ $e_{N}$ $(-)^{-1}_{N}$ $\alpha_{X,N}$ $\alpha_{Y,H}$ $N$ $H$ $Y$ $Y$ $\alpha_{X,H}$ $(-)^{-1}_{H}$ $k$ $\alpha_{Y,H}$ $X$ $\gamma$ $Y$ $H$ $=$ (17) $$

where on the right-hand side we have simplified the operations involving the identity term $e_{N}$ via the basic group axioms (e.g. $(-)^{-1}_{N}\circ e_{N}$ becomes simply $e_{N}$).
As shown by Cornish (), when $k:X\to Y$ is $H$-equivariant with respect to $\alpha_{X,N}$ and $\alpha_{Y,N}$, it always holds that $\operatorname{\mathsf{sym}}_{\gamma}(k)$ is equivariant with respect to the ($N\times H$)-actions $\alpha_{X}$ and $\alpha_{Y}$.

To prove Proposition 1 here, we now simply instantiate the above discussion in the Markov category $\mathsf{C}\coloneqq\mathsf{Stoch}$, so that all the components involved become Markov kernels.
Concretely, we take $X\coloneqq\mathcal{X}$, $Y\coloneqq\mathcal{Y}$, $N\coloneqq\mathcal{H}$, and $H\coloneqq\mathcal{G}$.
A procedure for sampling from the symmetrised kernel $\operatorname{\mathsf{sym}}(\gamma)(d\mathbf{y}|\mathbf{x})$ may now be read off from equation 17 as follows:

$$ $g\sim\gamma(dg|\mathbf{x})\qquad\mathbf{y}_{0}\sim k(d\mathbf{y}|g^{-1}\cdot \mathbf{x})\qquad\mathbf{y}\coloneqq g\cdot\mathbf{y}_{0},$ $$

as in the statement of this result.
∎

###### Proof.

### A.2 Proof of Proposition 2

###### Proof.

For simplicity, we assume $k(d\mathbf{y}|\mathbf{x})$ has a density $k(\mathbf{y}|\mathbf{x})$ with respect to the Lebesgue measure $\mu$ on $\mathcal{Y}=\mathbb{R}^{n}$. To derive the density of $k(d\mathbf{y}|g,\mathbf{x})$ where $g\in\mathcal{G}$ is some fixed group element, we note that for $\mathbf{y}\sim k(d\mathbf{y}|g,\mathbf{x})$, we have $\mathbf{y}=g\cdot\mathbf{y}_{0}$ where $\mathbf{y}_{0}\sim k(\mathbf{y}|g^{-1}\cdot\mathbf{x})$. Hence, by the change-of-variables formula, we can conclude that the density of $k(\mathbf{y}|g,\mathbf{x})$ exists and has the form:

$$ $k(\mathbf{y}|g,\mathbf{x})=k(g^{-1}\cdot\mathbf{y}|g^{-1}\cdot\mathbf{x})\left |\frac{\partial(g^{-1}\cdot\mathbf{y})}{\partial{\mathbf{y}}}\right|.$ $$

Hence, we see that when the action of $\mathcal{G}$ has unit Jacobian, the density of $k(\mathbf{y}|g,\mathbf{x})=k(g^{-1}\cdot\mathbf{y}|g^{-1}\cdot\mathbf{x})$.

Further, suppose we have $g\sim\gamma(dg|\mathbf{x})$ and $\mathbf{y}\sim k(\mathbf{y}|g,\mathbf{x})$ - i.e. $\mathbf{y}\sim\operatorname{\mathsf{sym}}_{\gamma}(k)(d\mathbf{y}|\mathbf{x})$.
It is the case that for an arbitrary (Borel) measurable set $A$ in $\mathcal{Y}$, we have

$$ $\operatorname{\mathsf{sym}}_{\gamma}(k)(\mathbf{y}\in A|\mathbf{x})=\int_{ \mathcal{G}}k(\mathbf{y}\in A|g,\mathbf{x})\,\gamma(dg|\mathbf{x}).$ $$

Since we have shown above that $k(d\mathbf{y}|g,\mathbf{x})$ has a density, we can express this as

$$ $\displaystyle\operatorname{\mathsf{sym}}_{\gamma}(k)(\mathbf{y}\in A|\mathbf{x})$ $\displaystyle=\int_{\mathcal{G}}\left(\int_{A}k(\mathbf{y}|g,\mathbf{x})\ \mu( d\mathbf{y})\right)\ \gamma(dg|\mathbf{x})$ $\displaystyle=\int_{A}\mathbb{E}_{\gamma(dg|\mathbf{x})}\left[k(\mathbf{y}|g, \mathbf{x})\right]\mu(d\mathbf{y}).$ $$

where we use Fubini’s theorem for the second line as all quantities are non-negative. Hence, we can conclude that the density of $\operatorname{\mathsf{sym}}_{\gamma}(k)$ exists and has the form $\operatorname{\mathsf{sym}}_{\gamma}(k)=\mathbb{E}_{\gamma(dg|\mathbf{x})}[k(
\mathbf{y}|g,\mathbf{x})]$.
∎

###### Proof.

### A.3 Proof of Proposition 3

###### Proof.

The standard diffusion objective with data augmentation distributed according to the Haar measure $\lambda$ is given by

$$ $\displaystyle\mathcal{L}_{t}^{\mathrm{aug}}$ $\displaystyle=\mathbb{E}_{q(\mathbf{z}_{0}),\epsilon\sim\mathcal{N}_{\mathcal{ U}}(0,\mathbf{I}),\lambda(dR)}\left[\frac{1}{2}w(t)\norm{\epsilon-\epsilon_{ \theta}(\alpha_{t}R\cdot\mathbf{z}_{0}+\sigma_{t}\epsilon)}^{2}\right]$ (18) $\displaystyle=\mathbb{E}_{q(\mathbf{z}_{0}),\epsilon\sim\mathcal{N}_{\mathcal{ U}}(0,\mathbf{I}),\lambda(dR)}\left[\frac{1}{2}w(t)\norm{\epsilon-\epsilon_{ \theta}(R\cdot(\alpha_{t}\mathbf{z}_{0}+\sigma_{t}R^{T}\cdot\epsilon))}^{2}\right]$ (19) $\displaystyle=\mathbb{E}_{q(\mathbf{z}_{0}),\epsilon^{\prime}\sim\mathcal{N}_{ \mathcal{U}}(0,\mathbf{I}),\lambda(dR)}\left[\frac{1}{2}w(t)\norm{R\cdot \epsilon^{\prime}-\epsilon_{\theta}(R\cdot(\alpha_{t}\mathbf{z}_{0}+\sigma_{t} \epsilon^{\prime}))}^{2}\right]$ (20) $\displaystyle=\mathbb{E}_{q(\mathbf{z}_{0}),\epsilon\sim\mathcal{N}_{\mathcal{ U}}(0,\mathbf{I}),\lambda(dR)}\left[\frac{1}{2}w(t)\norm{\epsilon-R^{T}\cdot \epsilon_{\theta}(R\cdot(\alpha_{t}\mathbf{z}_{0}+\sigma_{t}\epsilon))}^{2} \right],$ (21) $$

where we use the fact that $\epsilon^{\prime}=R^{T}\cdot\epsilon$ is distributed according to $\mathcal{N}_{\mathcal{U}}(0,\mathbf{I})$ as $R,\epsilon$ are independent in the expectation, and that the action of $R\in\operatorname{O}(3)$ preserves the L2 norm. To conclude, we note that it is a standard result that if $R\sim\lambda$, the inverse $R^{T}$ is also distributed according to the Haar measure $\lambda$. Hence, we see that $\mathcal{L}_{t}^{\mathrm{aug}}$ coincides with

$$ $\displaystyle\mathcal{L}_{t}=\mathbb{E}_{q(\mathbf{z}_{0}),\epsilon\sim \mathcal{N}_{\mathcal{U}}(0,\mathbf{I}),\lambda(dR)}\left[\frac{1}{2}w(t)\norm {\epsilon-R\cdot\epsilon_{\theta}(R^{T}\cdot\mathbf{z}_{t})}^{2}\right],$ (22) $$

where $\mathbf{z}_{t}=\alpha_{t}\mathbf{z}_{0}+\sigma_{t}\epsilon$.
∎

###### Proof.

## Appendix B Model Architecture

Below, we outline the architectures used for $\epsilon_{\theta}$ and $\gamma_{\theta}$. Both components rely on Diffusion Transformers (DiTs) (Peebles & Xie, ) using the official PyTorch implementation at [https://github.com/facebookresearch/DiT](https://github.com/facebookresearch/DiT). We also state the hyperparameters that we kept fixed for both our QM9 and GEOM-Drugs experiments. Any hyperparameters that differed between the datasets are discussed in their respective sections later in the Appendix.

We emphasise that our architecture choices were not extensively tuned as the main purpose of our experiments was to show that we can use generic architectures for equivariant diffusion models. We arrived at the below architecture through small adjustments from experimenting with DiT models in the context of molecular generation, which we stuck with in our final experiments.

### B.1 Architecture of ϵ θ subscript italic-ϵ 𝜃 \epsilon_{\theta} italic_ϵ start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT

As we need $\epsilon_{\theta}$ to be an $S_{N}$-equivariant architecture, we take $\epsilon_{\theta}$ to be parametrised in terms of a DiT model $\epsilon^{\prime}_{\theta}$ with $n_{\text{layers}}$ layers, $n_{\text{head}}$ attention heads and hidden size $n_{\text{size}}$ where to ensure the output lies in $\mathcal{U}$ we project the output via $\pi$. For the MLP layers, we use SwiGLU activations Shazeer () instead of the standard GELU, where the ratio of the hidden size of the SwiGLU to the hidden size of $\epsilon^{\prime}_{\theta}$ is 2. We do not use the default Fourier embeddings for the inputs and instead use linear layers to project our inputs to the correct size; for time embeddings, we use the default. We use Gaussian positional embeddings from Luo et al. () as additional features that we concatenate to our inputs. To compute this, we let

$$ $\psi^{k}_{(i,j)}=-\frac{1}{\sqrt{2\pi}|\sigma^{k}|}\exp\left(-\frac{1}{2}\left (\frac{\norm{\mathbf{x}^{(j)}-\mathbf{x}^{(j)}}-\mu^{k}}{|\sigma^{k}|}\right) \right),$ $$

where $k=1,\ldots,K$ is the number of basis kernels we use and $\mu^{k},\sigma^{k}\in\mathbb{R}$ are learnable parameters. We define $\mathbf{\psi}_{(i,j)}=(\psi_{(i,j)}^{1},\ldots,\psi_{(i,j)}^{K})^{T}\in\mathbb
{R}^{K\times 1}$. We then compute our positional embeddings by $\Psi_{i}=\frac{1}{N}\sum_{j=1}^{N}\mathbf{\psi}_{(i,j)}W_{D}$ where $W_{D}\in\mathbb{R}^{K\times n_{\text{emb}}}$ is a learnable matrix, and we concatenate these to form our embeddings $\mathbf{\Psi}=[\Psi_{1},\ldots,\Psi_{N}]\in\mathbb{R}^{N\times n_{\text{emb}}}$.

To compute $\epsilon_{\theta}$, given the input $\mathbf{z}=[\mathbf{x},\mathbf{h}]\in\mathbb{R}^{N\times(3+d)}$, we use a learnable linear layer $W_{I}\in\mathbb{R}^{(3+d)\times n_{z}}$ to project $\mathbf{z}$ to $\mathbf{z}^{\prime}=\mathbb{R}^{N\times n_{z}}$, and we compute the positional embeddings $\mathbf{\Psi}\in\mathbb{R}^{N\times n_{\text{emb}}}$, as outlined above, where $n_{\text{emb}}=n_{\text{size}}-n_{z}$. We then pass $[\mathbf{z}^{\prime},\mathbf{\Psi}]\in\mathbb{R}^{N\times n_{\text{size}}}$ to our DiT model $\epsilon_{\theta}^{\prime}$.

### B.2 Architecture of γ θ subscript 𝛾 𝜃 \gamma_{\theta} italic_γ start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT

We construct $\gamma_{\theta}$ following the recursive setup in Section 3.3. We take $f_{\theta}$ to be comprised of a DiT model with $m_{\text{layers}}$ layers, $m_{\text{head}}$ attention heads and hidden size $m_{\text{size}}$, and the same SwiGLU setup as above. To ensure that the architecture is $S_{N}$-invariant, we use the approach in Zaheer et al. (), where we further take the token-wise average of the representations from the DiT model (before the final output layer) to give a vector in $\mathbb{R}^{m_{\text{size}}}$. We then pass this into a small MLP $\mathbf{x}\mapsto\text{GELU}(\mathbf{x}W_{1})W_{2}$ where $W_{1}\in\mathbb{R}^{m_{\text{size}}\times m_{\text{size}}^{\prime}},W_{2}\in
\mathbb{R}^{m_{\text{size}}^{\prime}\times(3\times 3)}$ are learnable linear layers, which produces a matrix $\mathbb{R}^{3\times 3}$. We apply the QR decomposition to convert this into an orthogonal matrix which we then return.

To compute $\gamma_{\theta}$ given the input $\mathbf{z}=[\mathbf{x},\mathbf{h}]\in\mathbb{R}^{N\times(3+d)}$, we discard $\mathbf{h}$ and compute $\mathbf{\Psi}\in\mathbb{R}^{N\times n_{\text{emb}}}$ with the same parameters used in $\epsilon_{\theta}$ from $\mathbf{x}$. We also sample $\eta\sim\mathcal{N}_{\mathcal{U}}(0,\mathbf{I})$ of size $\mathbb{R}^{N\times m_{\text{noise}}}$ and an orthogonal matrix $R\in\mathbb{R}^{3\times 3}$ from the Haar measure on $\operatorname{O}(3)$. We then project the concatenated input $[R^{T}\cdot\mathbf{x},\eta,\mathbf{\Psi}]\in\mathbb{R}^{N\times(3+m_{\text{
noise}}+n_{\text{emb}})}$ by a learnable linear layer $W_{G}\in\mathbb{R}^{(3+m_{\text{noise}}+n_{\text{emb}})\times m_{\text{size}}}$ to $\mathbf{x}^{\prime}=\mathbb{R}^{N\times m_{\text{size}}}$ (we note that $\mathbf{\Psi}$ is invariant to $\operatorname{O}(3)$). We then pass $\mathbf{x}^{\prime}$ through $f_{\theta}$ and apply $R$ to the output - i.e. we return $R\cdot f_{\theta}(\mathbf{x}^{\prime})$.

### B.3 Hyperparameters for γ θ subscript 𝛾 𝜃 \gamma_{\theta} italic_γ start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT and ϵ θ subscript italic-ϵ 𝜃 \epsilon_{\theta} italic_ϵ start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT

For both QM9 and GEOM-Drugs we fixed the following hyperparameters. For $\epsilon_{\theta}$, we set $n_{\text{size}}=K\approx\frac{1}{2}n_{\text{size}}$. For $\gamma_{\theta}$, we set $m_{\text{noise}}=3$ and $m^{\prime}_{\text{size}}=4m_{\text{size}}$.

## Appendix C Experimental details

### C.1 Zeroth likelihood term ℒ 0 subscript ℒ 0 \mathcal{L}_{0} caligraphic_L start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT

We have a presented a framework for parametrising and optimising $p_{\theta}(\mathbf{z}_{t-1}|\mathbf{z}_{t})$ for $t>1$ in Section 3.3 where $p_{\theta}$ is obtained via stochastic symmetrisation. This corresponds to the $\mathcal{L}_{t}$ terms where $t>0$ in our objective in equation 9. We note however that standard diffusion models usually choose a different parametrisation for $p_{\theta}(\mathbf{z}_{0}|\mathbf{z}_{1})$ as this corresponds to the final generation step. Depending on the modelling task, this requires a different approach compared to the other reverse kernels. For example, in Hoogeboom et al. (), $p_{\theta}(\mathbf{z}_{0}|\mathbf{z}_{1})$ is defined as the product of densities $p_{\theta}^{\text{cont}}(\mathbf{x}_{0}|\mathbf{z}_{1})p_{\theta}^{\text{disc}
}(\mathbf{h}_{0}|\mathbf{z}_{1})$ where $p_{\theta}^{\text{disc}}$ implements a quantisation step converting the continuous latent $\mathbf{z}_{1}$ to discrete values $\mathbf{h}_{0}$, while $p_{\theta}^{\text{cont}}$ is still a Gaussian distribution generating continuous geometric features $\mathbf{x}_{0}$ from $\mathbf{z}_{1}$. In particular, we have that $p_{\theta}^{\text{cont}}(\mathbf{x}_{0}|\mathbf{z}_{1})=\mathcal{N}_{\mathcal{
U}}(\mathbf{x}_{0};\mathbf{x}_{1}/\alpha_{1}-\sigma_{1}/\alpha_{1}\epsilon_{
\theta}^{(\mathbf{x})}(\mathbf{z}_{1}),\sigma_{1}^{2}/\alpha_{1}^{2}\mathbf{I})$ where $\epsilon_{\theta}:\mathcal{Z}\to\mathcal{Z}$ is some $S_{N}$-invariant neural network and $\epsilon_{\theta}^{(\mathbf{x})}$ denotes the $\mathbf{x}$ component of the output of $\epsilon_{\theta}$.

We note that our above methodology can still account for this case by defining the symmetrised kernel in terms of

$$ $p_{\theta}(\mathbf{z}_{0}|\mathbf{z}_{1})=\operatorname{\mathsf{sym}}_{\gamma_ {\theta}}(k_{\theta})(\mathbf{z}_{0}|\mathbf{z}_{1}),\qquad k_{\theta}(\mathbf {z}_{0}|\mathbf{z}_{1})=p_{\theta}^{\text{cont}}(\mathbf{x}_{0}|\mathbf{z}_{1} )p_{\theta}^{\text{disc}}(\mathbf{h}_{0}|\mathbf{z}_{1})$ $$

We can follow the same discussion in Section 3.3 to conclude that

$$ $p_{\theta}(\mathbf{z}_{0}|\mathbf{z}_{1})=\mathbb{E}_{\gamma_{\theta}(dR| \mathbf{z}_{1})}\left[k_{\theta}(\mathbf{z}_{0}|R,\mathbf{z}_{1})\right], \qquad k_{\theta}(\mathbf{z}_{0}|R,\mathbf{z}_{1})=k_{\theta}^{\text{cont}}( \mathbf{x}_{0}|R,\mathbf{z}_{1})k_{\theta}^{\text{disc}}(\mathbf{h}_{0}|R, \mathbf{z}_{1}),$ $$

where $k_{\theta}^{\text{cont}}(\mathbf{z}_{0}|R,\mathbf{z}_{1})=\mathcal{N}(\mathbf{
x}_{0};\mathbf{x}_{1}/\alpha_{1}-\sigma_{1}/\alpha_{1}R\cdot\epsilon_{\theta}(
R^{T}\cdot\mathbf{z}_{1}),\sigma_{1}^{2}/\alpha_{1}^{2}\mathbf{I})$ and $k_{\theta}^{\text{disc}}(\mathbf{h}_{0}|R,\mathbf{z}_{1})=p_{\theta}^{\text{
disc}}(\mathbf{h}_{0}|R^{T}\cdot\mathbf{z}_{1})$. This allows us to decompose $\mathcal{L}_{0}$ into the form in equation 9 and to tractably optimise this objective since we have access to the density $k_{\theta}(\mathbf{z}_{0}|R,\mathbf{z}_{1})$.

### C.2 QM9

#### C.2.1 Model hyperparameters

##### SymDiff

Table 5 shows the hyperparameters for the $\epsilon_{\theta}$ backbone of the $k_{\theta}$ component of the SymDiff models used for QM9. The remaining hyperparameters were kept the same as in Appendix B.3.

**Table 5: Choice of $n_{\text{size}}$, $n_{\text{layers}}$, $n_{\text{heads}}$ for the $\epsilon_{\theta}$ of the SymDiff models used for QM9.**
| Model | # Parameters | $n_{\text{size}}$ | $n_{\text{layers}}$ | $n_{\text{heads}}$ |
| --- | --- | --- | --- | --- |
| SymDiff^∗ | 115.6M | 768 | 12 | 12 |
| SymDiff^+++ | 29M | 384 | 12 | 6 |
| SymDiff^++ | 21.3M | 360 | 12 | 6 |
| SymDiff^+ | 11.3M | 294 | 12 | 6 |

For the $f_{\theta}$ backbone of the $\gamma_{\theta}$ component, we set $m_{\text{size}}=128$, $m_{\text{layers}}=12$, $m_{\text{heads}}=6$.

##### EDM

Table 6 shows the hyperparameters for the EDM models that we used for our QM9 experiments. The remaining model hyperparameters were kept the same as those in Hoogeboom et al. ().

**Table 6: Choices of the hyperparameters $n_{f}$ (# features per layer), $n_{l}$ (number of layers) for the EDM models used for QM9.**
| Model | # Parameters | $n_{f}$ | $n_{l}$ |
| --- | --- | --- | --- |
| EDM^+++ | 12.4M | 332 | 12 |
| EDM^++ | 9.5M | 256 | 16 |
| EDM^+ | 5.3M | 256 | 9 |

#### C.2.2 Optimisation

For our optimisation, we followed Peebles & Xie () and used AdamW () with a batch size of 256. We chose a learning rate of $10^{-4}$ and weight decay of $10^{-12}$ for our 31.2M parameter model by searching over a small grid of 3 values for each. To match the same number of steps as in Hoogeboom et al. (), we trained our model for 4350 epochs.

We applied the same optimization hyperparameters from our 31.2M model to all other SymDiff models. For the EDM models, we followed the default hyperparameters from Hoogeboom et al. (). In our augmentation experiments, we first tuned the learning rate and weight decay for the DiT model, keeping all other optimization hyperparameters unchanged. These tuned values were then applied to DiT-Aug.

#### C.2.3 Comparisons with SymDiff and EDM

To ensure a fair comparison between the 5.3M parameter EDM and our 31.2M parameter SymDiff model we compared the performance of our SymDiff model without the $\gamma_{\theta}$ symmetrisation component (DiT) against the EDM.

**Table 7: Test NLL, atom stability, molecular stability, validity and uniqueness on QM9 for 10,000 samples and 3 evaluation runs.**
| Method | NLL $\downarrow$ | Atm. stability (%) $\uparrow$ | Mol. stability (%) $\uparrow$ | Val. (%) $\uparrow$ | Uniq. (%) $\uparrow$ |
| --- | --- | --- | --- | --- | --- |
| EDM | -110.70${\scriptstyle\pm{1.50}}$ | 98.70${\scriptstyle\pm{0.10}}$ | 82.00 ${\scriptstyle\pm{0.40}}$ | 91.90${\scriptstyle\pm{0.50}}$ | 90.70${\scriptstyle\pm{0.60}}$ |
| DiT | -127.78${\scriptstyle\pm{2.49}}$ | 98.23${\scriptstyle\pm{0.04}}$ | 81.03${\scriptstyle\pm{0.25}}$ | 94.71${\scriptstyle\pm{0.31}}$ | 97.98${\scriptstyle\pm{0.12}}$ |

From Table 7 we see that the EDM method performs similarly to the DiT model, with no clear indication of a significant difference between the two. In particular, we see that our DiT model underperforms on the main molecular metrics of atom stability and molecular stability. This suggests that the backbone of our SymDiff model was not expressive enough by its own to surpass the EDM, making the comparison between the symmetrised SymDiff model and the EDM a fair one.

#### C.2.4 Additional results

Here we restate the results from Table 8 but now include the results for SymDiff^∗ when trained for an additional 500 epochs (we refer to this as SymDiff-Long^∗). We see that SymDiff-Long^∗ outperforms all our baselines bar uniqueness. For NLL it is within variance with that of MUDiff.

**Table 8: Test NLL, atom stability, molecular stability, validity and uniqueness on QM9 for 10,000 samples and 3 evaluation runs. We omit the results for NLL where not available.**
| Method | NLL $\downarrow$ | Atm. stability (%) $\uparrow$ | Mol. stability (%) $\uparrow$ | Val. (%) $\uparrow$ | Uniq. (%) $\uparrow$ |
| --- | --- | --- | --- | --- | --- |
| GeoLDM | – | 98.90 ${\scriptstyle\pm{0.10}}$ | 89.40 ${\scriptstyle\pm{0.50}}$ | 93.80 ${\scriptstyle\pm{0.40}}$ | 92.70 ${\scriptstyle\pm{0.50}}$ |
| MUDiff | -135.50${\scriptstyle\pm{2.10}}$ | 98.80${\scriptstyle\pm{0.20}}$ | $89.90{\scriptstyle\pm{1.10}}$ | 95.30${\scriptstyle\pm{1.50}}$ | 99.10${\scriptstyle\pm{0.50}}$ |
| END | – | 98.90 ${\scriptstyle\pm{0.00}}$ | 89.10 ${\scriptstyle\pm{0.10}}$ | 94.80 ${\scriptstyle\pm{0.10}}$ | 92.60 ${\scriptstyle\pm{0.20}}$ |
| EDM | -110.70${\scriptstyle\pm{1.50}}$ | 98.70${\scriptstyle\pm{0.10}}$ | 82.00 ${\scriptstyle\pm{0.40}}$ | 91.90${\scriptstyle\pm{0.50}}$ | 90.70${\scriptstyle\pm{0.60}}$ |
| SymDiff-Long^∗ | -135.04${\scriptstyle\pm{0.64}}$ | 98.98${\scriptstyle\pm{0.05}}$ | 90.26${\scriptstyle\pm{0.38}}$ | 96.38${\scriptstyle\pm{0.27}}$ | 97.60${\scriptstyle\pm{0.10}}$ |
| SymDiff^∗ | -133.79${\scriptstyle\pm{1.33}}$ | 98.92${\scriptstyle\pm{0.03}}$ | 89.65${\scriptstyle\pm{0.10}}$ | 96.36${\scriptstyle\pm{0.27}}$ | 97.66${\scriptstyle\pm{0.22}}$ |
| Data |  | 99.00 | 95.20 | 97.8 | 100 |

### C.3 GEOM-Drugs

Like with QM9, we compared our method with data augmentation again and investigated its scaling behavior. We also looked to verify that the comparison between our 6.4M SymDiff model and the EDM from Hoogeboom et al. () was a fair one.

#### C.3.1 Model hyperparameters

Table 9 shows the hyperparameters for the $\epsilon_{\theta}$ used for the GEOM-Drugs experiments. The SymDiff^+ model corresponds to the 6.4M parameter EDM used in Table 4. The remaining hyperparameters were kept the same as B.3.

**Table 9: Choice of $n_{\text{size}}$, $n_{\text{layers}}$, $n_{\text{heads}}$ for the $\epsilon_{\theta}$ of the SymDiff models used for the GEOM-Drugs.**
| Model | # Parameters | $n_{\text{size}}$ | $n_{\text{layers}}$ | $n_{\text{heads}}$ |
| --- | --- | --- | --- | --- |
| SymDiff^+++ | 19.9M | 384 | 8 | 6 |
| SymDiff^++ | 9.1M | 264 | 8 | 6 |
| SymDiff^+ | 5.3M | 228 | 6 | 6 |

For the $f_{\theta}$ backbone of the $\gamma_{\theta}$ component, we set $m_{\text{size}}=128$, $m_{\text{layers}}=4$, $m_{\text{heads}}=4$ and kept the remaining hyperparameters the same as in Appendix B.3. Thee backbone consisted of 1.1M parameters.

#### C.3.2 Comparisons with SymDiff and EDM

To ensure a fair comparison between the EDM and our SymDiff model we compared the performance of our SymDiff model without the $\gamma_{\theta}$ symmetrisation component (DiT) against the EDM.

**Table 10: Test NLL, atom stability, molecular stability, validity and uniqueness on Geom-Drugs for 10,000 samples and 3 evaluation runs.**
| Method | NLL $\downarrow$ | Atm. stability (%) $\uparrow$ | Val. (%) $\uparrow$ |
| --- | --- | --- | --- |
| EDM | -137.10 | 81.30 | – |
| DiT | -200.57${\scriptstyle\pm{2.28}}$ | 83.38${\scriptstyle\pm{0.10}}$ | 95.66${\scriptstyle\pm{0.05}}$ |

In this case, we observe that our DiT model outperformed the EDM. Although this would imply a potentially unfair comparison between our approach and EDM, we still used the results in this case to highlight the relative improvements that we can gain from symmetrisation.

#### C.3.3 Scalability

We investigated the scalability of our approach by training our SymDiff model for the models in Table 9. Due to computational constraints, we did not train the EDM models and only make a comparison for the relative scaling of our SymDiff approach. For all our models, we used the same hyperparameters as for the 6.4M parameter model.

**Table 11: Total wall time, iterations per second sampling time for different sizes of SymDiff.**
| Method | # Parameters | Total wall-time | Iterations per second | Sampling time (s) |
| --- | --- | --- | --- | --- |
| SymDiff^+++ | 21M | 1 day 8h | 7.72 | 0.35 |
| SymDiff^++ | 10.2M | 1 day 1h | 9.84 | 0.26 |
| SymDiff^+ | 6.4M | 20h | 12.28 | 0.21 |

From Table 11 we see that our method scales well. Additionally, from Table 12 we see that as we scale our model, our performance becomes more competitive with the more sophisticated baselines. Like with QM9, we conjecture that similar performance improvements can be achieved by using our SymDiff approach as a drop-in replacement for the reverse kernels in more sophisticated methods.

**Table 12: Test NLL, atom stability and validity on GEOM-Drugs for 10,000 samples and 3 evaluation runs.**
| Method | NLL $\downarrow$ | Atm. stability (%) $\uparrow$ | Val. (%) $\uparrow$ |
| --- | --- | --- | --- |
| GeoLDM | – | 84.4 | 99.3 |
| END | – | 87.8${\scriptstyle\pm{0.99}}$ | 92.9${\scriptstyle\pm{0.3}}$ |
| EDM | -137.1 | 81.3 | – |
| SymDiff^+++ | -266.88${\scriptstyle\pm{0.43}}$ | 86.16${\scriptstyle\pm{0.04}}$ | 98.18${\scriptstyle\pm{0.08}}$ |
| SymDiff^++ | -230.12${\scriptstyle\pm{1.45}}$ | 85.34${\scriptstyle\pm{0.08}}$ | 96.78${\scriptstyle\pm{0.09}}$ |
| SymDiff^+ | -196.77${\scriptstyle\pm{0.03}}$ | 84.15${\scriptstyle\pm{0.07}}$ | 95.99${\scriptstyle\pm{0.25}}$ |
| Data |  | 86.50 | 99.9 |

#### C.3.4 Comparisons with data augmentation

To compare our approach with data augmentation, we trained an unsymmetrised reverse process with a 5.3M parameter DiT backbone (DiT), as well as the same model using data augmentation (DiT-Aug).
We also compared our SymDiff model to the case where we simply let $\gamma_{\theta}$ be obtained using the Haar measure on $\operatorname{O}(3)$, rather than learning this component (SymDiff-H).

**Table 13: Test NLL, atom stability and validity on GEOM-Drugs for 10,000 samples and 3 evaluation runs. We omit the results for NLL and validity where not available.**
| Method | NLL $\downarrow$ | Atm. stability (%) $\uparrow$ | Val. (%) $\uparrow$ |
| --- | --- | --- | --- |
| SymDiff | -196.77${\scriptstyle\pm{0.03}}$ | 84.15${\scriptstyle\pm{0.07}}$ | 95.99${\scriptstyle\pm{0.25}}$ |
| SymDiff-H | -189.48 ${\scriptstyle\pm{1.26}}$ | 84.06${\scriptstyle\pm{0.03}}$ | 94.51${\scriptstyle\pm{0.17}}$ |
| DiT-Aug | -191.13${\scriptstyle\pm{0.70}}$ | 83.95${\scriptstyle\pm{0.11}}$ | 94.26${\scriptstyle\pm{0.31}}$ |
| DiT | -200.57${\scriptstyle\pm{2.28}}$ | 83.38${\scriptstyle\pm{0.10}}$ | 95.66${\scriptstyle\pm{0.05}}$ |
| Data |  | 86.50 | 99.9 |

From Table 13, we see that our SymDiff approach outperforms both our augmentation baselines, DiT-Aug and SymDiff-H. This further reinforces the benefit of our approach as a more effective way of incorporating $\operatorname{O}(3)$ equivariance than data augmentation.

### C.4 Pretrain-finetuning

To further explore the flexibility of our approach, we experimented with using it in the pretrain-finetune framework, similar to Mondal et al. (). Using QM9, we took the trained DiT model from Table 1 and substituted it as the $\epsilon_{\theta}$ for our SymDiff model, while keeping the same architecture and hyperparameters for $f_{\theta}$. We tested two setups: finetuning both $epsilon_{\theta}$ and $f_{\theta}$ (DiT-FT) and freezing $\epsilon_{\theta}$ while tuning only $f_{\theta}$ (DiT-FT-Freeze). The same training procedure and optimization hyperparameters were used, except we now trained our models for 800 epochs and used a larger grid for learning rate and weight decay tuning. Specifically, we searched first for the optimal learning rate in $[10^{-3},8\times 10^{-4},2\times 10^{-4},10^{-4}]$ and for the optimal weight decay in $[0,10^{-12},2\times 10^{-12}]$. We found the optimal learning rate and weight decay to be $10^{-3}$ and $2\times 10^{-12}$.

**Table 14: Test NLL, atom stability, molecular stability, validity and uniqueness on QM9 for 10,000 samples and 3 evaluation runs.**
| Method | NLL $\downarrow$ | Atm. stability (%) $\uparrow$ | Mol. stability (%) $\uparrow$ | Val. (%) $\uparrow$ | Uniq. (%) $\uparrow$ |
| --- | --- | --- | --- | --- | --- |
| SymDiff | -129.35${\scriptstyle\pm{1.07}}$ | 98.74${\scriptstyle\pm{0.03}}$ | 87.49${\scriptstyle\pm{0.23}}$ | 95.75${\scriptstyle\pm{0.10}}$ | 97.89${\scriptstyle\pm{0.26}}$ |
| DiT-FT | -111.66${\scriptstyle\pm{1.22}}$ | 98.43${\scriptstyle\pm{0.03}}$ | 83.27 ${\scriptstyle\pm{0.39}}$ | 94.19${\scriptstyle\pm{0.16}}$ | 98.17${\scriptstyle\pm{0.26}}$ |
| DiT-FT-Freeze | -43.29${\scriptstyle\pm{3.73}}$ | 95.68${\scriptstyle\pm{0.02}}$ | 55.02${\scriptstyle\pm{0.38}}$ | 90.48${\scriptstyle\pm{0.24}}$ | 99.06${\scriptstyle\pm{0.13}}$ |
| DiT | -127.78${\scriptstyle\pm{2.49}}$ | 98.23${\scriptstyle\pm{0.04}}$ | 81.03${\scriptstyle\pm{0.25}}$ | 94.71${\scriptstyle\pm{0.31}}$ | 97.98${\scriptstyle\pm{0.12}}$ |

From Table 14, we observe that finetuning both $\epsilon_{\theta}$ and $f_{\theta}$ improves performance over the DiT model, even with our minimal optimization tuning. However, finetuning only $f_{\theta}$ leads to worse results, indicating that end-to-end training or finetuning the whole model is necessary. This underscores the flexibility of our approach and its potential for easy and efficient symmetrisation of pretrained DiT models with an unconstrained $f_{\theta}$.

## Appendix D Extension to Score Matching and Flow Matching

In this section, we discuss how to extend stochastic symmetrisation to score and flow-based generative models to give an analogue of SymDiff to these paradigms. For clarity of presentation, we consider all models to be defined for $N$-body systems living in the full space $\mathcal{Z}=\mathbb{R}^{N\times 3}$ where we wish to obtain a $S_{N}\times\operatorname{O}(3)$-equivariant model - i.e. we do not consider non-geometric features or translation invariance. Although, we note that the below discussion can be extended to such settings in the natural way as presented for diffusion models above.

### D.1 Score Matching

Score-based generative models (SGMs) (Song et al., ) are the continuous-time analogue of diffusion models. SGMs consider the forward noising process $\mathbf{x}_{t}\sim p_{t}$ for $t\in[0,T]$ defined by the following stochastic differential equation (SDE) with the initial condition $\mathbf{x}_{0}\sim p_{\text{data}}$:

$$ $\displaystyle d\mathbf{x}_{t}=f(\mathbf{x}_{t},t)dt+g(t)d\mathbf{w},$ (23) $$

for some choice of functions $f:\mathcal{Z}\times[0,T]\to\mathcal{Z}$ and $g:[0,T]\to\mathbb{R}$, and where $\mathbf{w}$ is a standard Weiner process.

The corresponding backward process is shown in Anderson () to take the form:

$$ $\displaystyle d\mathbf{x}_{t}=\left[f(\mathbf{x}_{t},t)-g(t)^{2}\nabla_{ \mathbf{x}}\log p_{t}(\mathbf{x}_{t})\right]dt+g(t)d\bar{\mathbf{w}},$ (24) $$

where $\bar{\mathbf{w}}$ is a standard Weiner process and time runs backwards from $T$ to 0 0. Hence, given samples $\mathbf{x}_{T}\sim p_{T}$ and access to the score of the marginal distributions $\nabla_{\mathbf{x}}\log p_{t}(\mathbf{x}_{t})$, we can obtain samples from $p_{\text{data}}$ by simulating the backward process in equation 24.

By considering the Euler–Maruyama discretisation of equation 24, we can represent the sampling scheme of a SGM in terms of the Markov chain $p_{T}(\mathbf{x}_{T})\prod_{i=1}^{n}p(\mathbf{x}_{t_{i-1}}|\mathbf{x}_{t_{i}})$, where the time-points $t_{i}$ are uniformly spaced in $[0,T]$ - i.e. $t_{i}=i\Delta t$ where $\Delta t=T/n$ - and the reverse transition kernels are given by:

$$ $\displaystyle p(\mathbf{x}_{t_{i-1}}|\mathbf{x}_{t_{i}})=\mathcal{N}\left( \mathbf{x}_{t_{i-1}};\mathbf{x}_{t_{i}}+\Delta t\left\{f(\mathbf{x}_{t_{i}},t_ {i})-g(t_{i})^{2}\nabla_{x}\log p_{t_{i}}(\mathbf{x}_{t_{i}})\right\},g(t_{i}) ^{2}\Delta t\mathbf{I}\right).$ $$

In what follows, we additionally assume that $f(\cdot,t)$ is $S_{N}\times\operatorname{O}(3)$-equivariant for all $t\in[0,T]$. This is true for common choices of $f$ which take $f$ to be linear in $\mathbf{x}_{t}$.

##### Stochastic symmetrisation

In order to learn an approximation to the transition kernels via stochastic symmetrisation, we can parametrise the reverse transition kernels, in a similar fashion as for diffusion models, by

$$ $\displaystyle p_{\theta}(\mathbf{x}_{t_{i-1}}|\mathbf{x}_{t_{i}})= \operatorname{\mathsf{sym}}_{\gamma_{\theta}}(k_{\theta})(\mathbf{x}_{t_{i-1}} |\mathbf{x}_{t_{i}}),$ $$

where we take $k_{\theta}(\mathbf{x}_{t_{i-1}}|\mathbf{x}_{t_{i}})=\mathcal{N}(\mathbf{x}_{t_
{i-1}};\mu_{\theta}(\mathbf{x}_{t_{i}}),g(t_{i})^{2}\Delta t\mathbf{I})$. We define $\mu_{\theta}$ by the following parametrisation(^4^44Similar to our discussion on diffusion models, we leave the time dependency implicit in here.):

$$ $\displaystyle\mu_{\theta}(\mathbf{x}_{t_{i}})=\mathbf{x}_{t_{i}}+\Delta t\left \{f(\mathbf{x}_{t_{i}},t_{i})-g(t_{i})^{2}s_{\theta}(\mathbf{x}_{t_{i}})\right\},$ $$

where we take $s_{\theta}:\mathcal{Z}\to\mathcal{Z}$ to be a $S_{N}$-equivariant neural network which aims to learn an approximation to the true score $\nabla_{\mathbf{x}}\log p_{t}(\mathbf{x}_{t})$. This ensures that $k_{\theta}$ is $S_{N}$-invariant. Additionally, we assume that $\gamma_{\theta}:\mathcal{Z}\to\operatorname{O}(3)$ is some choice of a $S_{N}$-invariant and $\operatorname{O}(3)$-equivariant Markov kernel. Hence, we can conclude that $p_{\theta}:\mathcal{Z}\to\mathcal{Z}$ is a $S_{N}\times\operatorname{O}(3)$-equivariant Markov kernel by Proposition 1. We can also guarantee that $p_{\theta}$ admits a density by Proposition 2.

##### Training

To learn $\theta$ for $p_{\theta}(\mathbf{x}_{t_{i-1}}|\mathbf{x}_{t_{i}})$,
a natural objective is to minimise the KL divergence between the true reverse kernels and our parametrised reverse kernels

$$ $\displaystyle\mathcal{L}(\theta)=\sum_{i=1}^{n}\lambda_{0}(t_{i})\mathcal{L}_{ i}(\theta),\qquad\mathcal{L}_{i}(\theta)=\mathbb{E}_{p_{t_{i}}(\mathbf{x}_{t_{ i}})}\left[D_{\mathrm{KL}}(p(\mathbf{x}_{t_{i-1}}|\mathbf{x}_{t_{i}})||p_{ \theta}(\mathbf{x}_{t_{i-1}}|\mathbf{x}_{t_{i}}))\right],$ $$

where $\lambda_{0}$ is some time weighting function.
We note that we run into the same issue as with SymDiff in that we do not have access to $p_{\theta}(\mathbf{x}_{t_{i-1}}|\mathbf{x}_{t_{i}})$ in closed-form, since this is expressed in terms of an expectation.
However, as $\mathcal{L}_{i}(\theta)$ is a linear function of $-\log p_{\theta}(\mathbf{x}_{t_{i-1}}|\mathbf{x}_{t_{i}})$ and $-\log$ is a convex function, we can apply Jensen’s inequality again to provide the following upper bound to our original objective

$$ $\displaystyle\mathcal{L}^{\prime}(\theta)=\sum_{i=1}^{n}\lambda_{0}(t_{i}) \mathcal{L}_{i}^{\prime}(\theta),\qquad\mathcal{L}_{i}^{\prime}(\theta)= \mathbb{E}_{p_{t_{i}}(\mathbf{x}_{t_{i}}),\gamma_{\theta}(dR|\mathbf{x}_{t_{i} })}\left[D_{\mathrm{KL}}(p(\mathbf{x}_{t_{i-1}}|\mathbf{x}_{t_{i}})||k_{\theta }(\mathbf{x}_{t_{i-1}}|R,\mathbf{x}_{t_{i}}))\right],$ $$

which we can now use to train $\theta$. To further simplify $\mathcal{L}_{i}^{\prime}(\theta)$, we note that $k_{\theta}(\mathbf{x}_{t_{i-1}}|R,\mathbf{x}_{t_{i}})=\mathcal{N}(\mathbf{x}_{
t_{i-1}};R\cdot\mu_{\theta}(R^{T}\cdot\mathbf{x}_{t_{i}}),g(t_{i})^{2}\Delta t
\mathbf{I})$ with a similar derivation as before.
This allows us to evaluate the KL divergences in closed form since $p,k_{\theta}$ are defined in terms of Gaussians. We can show that this gives

$$ $\displaystyle\mathcal{L}_{i}^{\prime}(\theta)=\mathbb{E}_{p_{t_{i}}(\mathbf{x} _{t_{i}}),\gamma_{\theta}(dR|\mathbf{x}_{t_{i}})}\left[\frac{1}{2}g(t_{i})^{2} \Delta t\norm{R\cdot s_{\theta}(R^{T}\cdot\mathbf{x}_{t_{i}})-\nabla_{\mathbf{ x}}\log p_{t_{i}}(\mathbf{x}_{t_{i}})}^{2}\right],$ (25) $$

where we use the fact that $f(\cdot,t)$ is $S_{N}\times\operatorname{O}(3)$-equivariant.
To express equation 25 in a tractable form (as we do not have access to the true score), we can apply the standard technique of employing the score matching identity (Vincent, ) to give

$$ $\displaystyle\mathcal{L}_{i}^{\prime}(\theta)=\mathbb{E}_{p(\mathbf{x}_{0}),p( \mathbf{x}_{t_{i}}|\mathbf{x}_{0}),\gamma_{\theta}(dR|\mathbf{x}_{t_{i}})} \left[\frac{1}{2}g(t_{i})^{2}\Delta t\norm{R\cdot s_{\theta}(R^{T}\cdot\mathbf {x}_{t_{i}})-\nabla_{\mathbf{x}}\log p_{t_{i}}(\mathbf{x}_{t_{i}}|\mathbf{x}_{ 0})}^{2}\right]+C_{i},$ $$

where $p(\mathbf{x}_{t_{i}}|\mathbf{x}_{0})$ denotes the conditional distribution of $\mathbf{x}_{t_{i}}$ given $\mathbf{x}_{0}$ under the forward noising process $p$, and $C_{i}$ is some constant. In practice, the choice of forward noising SDE in equation 23 is made to ensure that we have access to $p(\mathbf{x}_{t_{i}}|\mathbf{x}_{0})$ in closed-form and that the distribution is easy to sample from.

By making the choice that $\lambda_{0}(t)=2\lambda(t)/(g(t_{i})^{2}T)$ for some suitable time weighting function $\lambda$, we can show as $\Delta t\to 0$ that our objective $\mathcal{L}^{\prime}(\theta)$ will converge to (modulo some constant)

$$ $\displaystyle\mathbb{E}_{p(\mathbf{x}_{0}),t\sim U(0,T),p(\mathbf{x}_{t}| \mathbf{x}_{0}),\gamma_{\theta}(dR|\mathbf{x}_{t})}\left[\lambda(t)\norm{R \cdot s_{\theta}(R^{T}\cdot\mathbf{x}_{t})-\nabla_{\mathbf{x}}\log p(\mathbf{x }_{t}|\mathbf{x}_{0})}^{2}\right],$ (26) $$

where $U(0,T)$ denotes the uniform distribution on $[0,T]$. We see that our final objective in equation 26 now resembles the standard score matching objective.

### D.2 Flow Matching

Continuous normalising flows (CNFs) (Chen et al., ) construct a generative model of data $\mathbf{x}_{1}\sim q=p_{\text{data}}$ by the pushforward of an ordinary differential equation (ODE) taking the form

$$ $\displaystyle\frac{d}{dt}\phi_{t}(\mathbf{x})=u_{t}(\phi_{t}(\mathbf{x})), \qquad\phi_{0}(\mathbf{x})=\mathbf{x},$ (27) $$

where $u_{t}:\mathcal{Z}\times[0,T]\to\mathcal{Z}$ is the vector field function defining the ODE, and $\phi_{t}:\mathcal{Z}\times[0,T]\to\mathcal{Z}$ denotes the flow implicitly defined by solutions to the above ODE. By letting $p_{0}$ be some simple prior distribution, the above ODE defines a generative model $\mathbf{x}_{t}\sim p_{t}$ by the pushforward of $p_{0}$ through the flow $\phi_{t}$

$$ $\displaystyle p_{t}=[\phi_{t}]_{\#}p_{0}$ (28) $$

If $u_{t}$ is chosen in such a way that $p_{1}\approx q=p_{\text{data}}$, we can then generate samples from $p_{\text{data}}$ by sampling some $\mathbf{x}_{0}\sim p_{0}$, then solving the ODE in equation 27 with this initial condition(^5^55The use of time here reverse the convention used in the diffusion literature.). Furthermore, as in previous work (Klein et al., ), we assume that $u_{t}$ is $S_{N}\times\operatorname{O}(3)$-equivariant.

By considering the Euler discretisation of equation 27, we can represent the generation process by the Markov chain $p_{0}(\mathbf{x}_{0})\prod_{i=1}^{T}p(d\mathbf{x}_{t_{i}}|\mathbf{x}_{t_{i-1}})$ where the time-points $t_{i}$ are uniformly spaced in $[0,T]$ - i.e. $t_{i}=i\Delta t$ where $\Delta t=T/n$ - and the transition kernels are given by the Markov kernels

$$ $\displaystyle p(d\mathbf{x}_{t_{i}}|\mathbf{x}_{t_{i-1}})=\delta(\mathbf{x}_{t _{i-1}}+u_{t_{i-1}}(\mathbf{x}_{t_{i-1}})\Delta t),$ (29) $$

where $\delta(\cdot)$ denotes the Dirac measure at some point.

##### Stochastic symmetrisation

To learn the transition kernels induced by the vector field $u_{t}$ via stochastic symmetrisation, we parametrise our transition kernels by

$$ $\displaystyle p_{\theta}(d\mathbf{x}_{t_{i}}|\mathbf{x}_{t_{i-1}})= \operatorname{\mathsf{sym}}_{\gamma_{\theta}}(k_{\theta})(\mathbf{x}_{t_{i}}| \mathbf{x}_{t_{i-1}}),$ (30) $$

where we take $k_{\theta}(d\mathbf{x}_{t_{i}}|\mathbf{x}_{t_{i-1}})=\delta(\mathbf{x}_{t_{i-1
}}+v^{\theta}_{t_{i-1}}(\mathbf{x}_{t_{i-1}})\Delta t)$ in which $v_{t}^{\theta}:\mathcal{Z}\to\mathcal{Z}$ is some $S_{N}$-equivariant neural network which aims to learn an approximation to the true vector field $u_{t}$. We further assume $\gamma_{\theta}:\mathcal{Z}\to\operatorname{O}(3)$ is some $S_{N}$-invariant and $\operatorname{O}(3)$-equivariant Markov kernel. We can again conclude that $p_{\theta}:\mathcal{Z}\to\mathcal{Z}$ is a $S_{N}\times\operatorname{O}(3)$-equivariant Markov kernel by Proposition 1.

##### Training

A natural objective to learn $p_{\theta}$ is to minimise the 2-Wasserstein distance $\mathcal{W}_{2}$ (Peyré et al., ) between $p(d\mathbf{x}_{t_{i}}|\mathbf{x}_{t_{i-1}})$ and $p_{\theta}(d\mathbf{x}_{t_{i}}|\mathbf{x}_{t_{i-1}})$ since these are defined in terms of Dirac measures. We can write our objective as

$$ $\displaystyle\mathcal{L}(\theta)=\sum_{i=1}^{T}\lambda_{0}(t_{i-1})\mathcal{L} _{i}(\theta),\qquad\mathcal{L}_{i}(\theta)=\mathbb{E}_{p_{t_{i-1}}(\mathbf{x}_ {t_{i-1}})}\left[\mathcal{W}_{2}^{2}(p(d\mathbf{x}_{t_{i}}|\mathbf{x}_{t_{i-1} }),p_{\theta}(d\mathbf{x}_{t_{i}}|\mathbf{x}_{t_{i-1}}))\right]$ $$

where $\lambda_{0}$ is some time weighting function, and the 2-Wasserstein distance $\mathcal{W}_{2}$ is defined as $\mathcal{W}_{2}^{2}(\pi_{1},\pi_{2})=\inf_{\pi}\int\norm{\mathbf{x}-\mathbf{y}
}^{2}\ d\pi(\mathbf{x},\mathbf{y})$ where $\pi$ is taken over the space of possible couplings between the measures $\pi_{1},\pi_{2}$. We note that as $p(d\mathbf{x}_{t_{i}}|\mathbf{x}_{t_{i-1}})$ is Dirac, there only exists a single coupling between the two kernels given by the product of the Markov kernels. This allows us to evaluate $\mathcal{L}_{i+1}$ as

$$ $\displaystyle\mathcal{L}_{i+1}(\theta)=\mathbb{E}_{p_{t_{i}}(\mathbf{x}_{t_{i} }),\gamma_{\theta}(dR|\mathbf{x}_{t_{i}})}\left[\Delta t^{2}\norm{R\cdot v_{t_ {i}}^{\theta}(R^{T}\cdot\mathbf{x}_{t_{i}})-u_{t_{i}}(\mathbf{x}_{t_{i}})}^{2} \right].$ (31) $$

To express equation 31 in a tractable form (as we do not have access to $u_{t}$), we can take $u_{t}$ to be constructed by the same setup used in Flow Matching (Lipman et al., ). This framework allows us to express equation 31 in the now tractable form

$$ $\displaystyle\mathcal{L}_{i+1}(\theta)=\mathbb{E}_{q(\mathbf{x}_{1}),p_{t_{i}} (\mathbf{x}_{t_{i}}|\mathbf{x}_{1}),\gamma_{\theta}(dR|\mathbf{x}_{t_{i}})} \left[\Delta t^{2}\norm{R\cdot v_{t_{i}}^{\theta}(R\cdot\mathbf{x}_{t_{i}})-u_ {t_{i}}(\mathbf{x}_{t_{i}}|\mathbf{x}_{1})}\right]+C_{i+1},$ (32) $$

by the use of the Conditional Flow Matching objective, where $C_{i+1}$ is some constant. Here $p_{t}(\mathbf{x}_{t}|\mathbf{x}_{1})$ is a family of conditional distributions where $p_{0}(\mathbf{x}_{0}|\mathbf{x}_{1})=p_{0}(\mathbf{x}_{0})$ equals our prior distribution and $p_{1}(\mathbf{x}_{1}|\mathbf{x}_{1})\approx\delta(\mathbf{x}_{1})$, and for which $u_{t}(\mathbf{x}_{t}|\mathbf{x}_{1})$ is a vector field generating $p_{t}(\mathbf{x}_{t}|\mathbf{x}_{1})$ by an ODE of the form in equation 28. These are constructed to be easy to sample from and evaluate. The true vector field $u_{t}$, which provides a generative model of $q=p_{\text{data}}$, is then defined by some expectation of the conditional vector fields $u_{t}(\mathbf{x}_{t}|\mathbf{x}_{0})$ over $p_{t}(\mathbf{x}_{t}|\mathbf{x}_{0})$ and $q(\mathbf{x}_{1})$.

Hence, by taking $\lambda_{0}(t)=\frac{\lambda(t)}{T\Delta t}$ for some suitable time weighting function $\lambda$, we can show as $\Delta t\to 0$, our objective $\mathcal{L}(\theta)$ will converge to (modulo some constant)

$$ $\displaystyle\mathbb{E}_{q(\mathbf{x}_{1}),t\sim U(0,T),p_{t}(\mathbf{x}_{t}| \mathbf{x}_{1}),\gamma_{\theta}(dR|\mathbf{x}_{t})}\left[\lambda(t)\norm{R \cdot v_{t}^{\theta}(R^{T}\cdot\mathbf{x}_{t})-u_{t}(\mathbf{x}_{t}|\mathbf{x} _{1})}^{2}\right].$ (33) $$

We see that our final objective in equation 33 now resembles the standard flow matching objective.

## Peer Reviews

confidence: 4
contribution: 3
initial_rating: 6
presentation: 3
questions: I have put almost all my questions in the ''Weaknesses'' section above already, but I have one very minor additional question. In 4.2, the authors write that, given the size of GEOM-Drugs, they train a smaller model than on QM9. But the GEOM-Drugs dataset is larger and more complex than QM9. Wouldn't one want to train a larger, more expressive model in that case, and not a smaller one?
soundness: 3
strengths: **Clarity**: The paper is overall well written and the necessary background is introduced in an appropriate manner, such that one can follow the overall paper. However, the stochastic symmetrization framework is introduced in a fairly abstract manner and I believe more intuitions could be discussed and examples given, which would likely help the broad ICLR audience (see questions below).

**Originality**: Building equivariant diffusion models with stochastic symmetrization is an overall novel and original idea. Methodologically, the approach directly builds on the stochastic symmetrization framework from Cornish (2024). Hence, the novelty here lies in the application of the framework to diffusion models, specifically those for small molecule generation, a popular and important benchmark task.

**Experimental Results**: Overall, the paper is able to show favourable results compared to various baselines, thereby successfully validating the method, and I liked the many ablation experiments that were run.

**Quality and Significance**: Overall, this paper is of high quality and tackles an important problem. Enforcing equivariances and invariances in generative models can be crucial in applications with symmetries, as commonly encountered in scientific applications. However, it is difficult to directly build high-performance and scalable equivariant network architectures. Although quite related to data augmentation methods, stochastic symmetrization and SymDiff proposes a new way to address the problem, and I believe this paper has the potential to result in follow-up work. Hence, I believe this is sufficiently significant for ICLR.
summary: This paper proposes a new way to build equivariant diffusion models. Whereas traditional methods rely either on simple data augmentation strategies or neural network architectures that are directly constructed to be equivariant, *SymDiff* uses stochastic symmetrization to learn equivariant diffusion models. Stochastic symmetrization is a recently introduced framework by Cornish (2024). On a high level, stochastic symmetrization can intuitively be seen as an advanced data augmentation method, where the augmentation is sampled from an additional, learnable Markov kernel and where the main denoiser network does not need to be equivariant itself. Moreover, the augmentation is also applied during sampling. The paper in detail introduces the approach and rigorously shows that it is principled. SymDiff is then validated on molecule generation, where the permutation equivariance is directly built into the employed networks, whereas rotation equivariance is enforced through stochastic symmetrization. The paper compares to various baselines, showing favourable results, and runs many appropriate and insightful ablation experiments.
weaknesses: I also see some weaknesses.

- *More discussion of intuitions and more experimental analyses*: I think it would be great if the authors discussed some intuitions in more detail and analyzed their models in more detail. Specifically:
  - In practice, training SymDiff boils down to the training procedure in Algorithm 1. $\gamma_\theta$ is trained as part of the ELBO objective, using reparametrization. What sort of $\gamma_\theta$ and $f_\theta$ functions are we intuitively expected to learn in different situations? How does the model intuitively exploit its flexibility to produce more complex $R$ distributions compared to when only using the plain Haar measure? Can examples and intuitions be given?
  - Moreover, is there anything preventing the model from simply learning $R = R_0 f_\theta(...) = I$, this is, no ''augmentation'' at all? In other words, could $f_\theta$ learn to ''undo'' the $R_0$ drawn from the Haar measure so as to remove any rotations? In that case, the training would reduce to standard training without the rotation equivariance, which may be beneficial from the objective's perspective considering that the network itself isn't equivariant? But this defeats the purpose of the method. Such potential edge cases should be discussed.
  - To this end, can the authors maybe show and analyze what sort of $\gamma_\theta$ distributions SymDiff ends up learning in practice?
  - The approach uses $\gamma_\theta$ (or $f_\theta$) also during sampling. What is the intuition for that? And what would happen if we didn't apply $\gamma_\theta$ during sampling? This would be an interesting discussion or ablation. Moreover, in the case without learnable $\gamma_\theta$, i.e. when using only the Haar measure, I would expect that the model would still perform fine even if no $R$ is applied during sampling.

- *Relation to data augmentation*: Related to the last point above, overall, it seems that SymDiff can be seen as a sophisticated data augmentation strategy (as also pointed out by the authors). However, the authors carefully distinguish the case between training with data augmentation (DiT-Aug) and SymDiff with Haar measure sampling only (SymDiff-H), pointing out that the only difference in that case is that SymDiff-H also applies augmentations during sampling, if I understood things correctly. However, in all experiments the performance of DiT-Aug and SymDiff-H is extremely close together. This makes me believe the application of the augmentation during sampling isn't necessary for SymDiff-H (also see point above), and that SymDiff-H really just corresponds to standard data augmentation. Related to that, it is not sufficiently well discussed, explained and analyzed how exactly this situation *qualitatively* changes when $\gamma_\theta$ is learnt (see point above).

- *Experimental results*: The performance gains over the EDM baseline are meaningfully large, but at the same time compared to more sophisticated baseline models SymDiff does not perform much better. Further, the authors make the point that the comparison to EDM is fair. However, the very different network sizes still make a fair comparison somewhat difficult. How would performance between EDM and SymDiff compare for similar-sized neural networks?

**Conclusion:** SymDiff is an interesting and novel method, tackling a significant problem, and I believe it is of interest to the ICLR audience. However, there are various fundamental questions about the method and the experiments (see above), that I hope the authors can answer and address. For now, I am leaning towards suggesting acceptance.

---

confidence: 3
contribution: 4
initial_rating: 8
presentation: 4
questions: - Do you have insights into how many training epochs are necessary for SYMDIFF to effectively learn equivariance in the molecular generation task?
- How does the convergence speed of SYMDIFF compare to that of intrinsically equivariant models in learning equivariant properties?
- In the context of molecular generation, do you consider the equivariant component of the task to be a minor or major challenge compared to learning the underlying chemical laws necessary for accurate atom placement?
soundness: 4
strengths: This paper introduces SYMDIFF, a simple yet powerful method for constructing equivariant diffusion models without relying on complex intrinsically equivariant neural network architectures. The strengths of the paper can be assessed across the following dimensions:

Originality:
The paper presents a novel application of stochastic symmetrisation to generative modeling, specifically in the context of diffusion models. This approach is original because it departs from traditional methods that require specialized equivariant architectures with complex parameterizations and higher-order geometric features. By leveraging stochastic symmetrisation as a learned data augmentation technique applied during sampling, the paper creatively combines existing ideas to overcome limitations of prior models. This opens up new possibilities for using standard neural network architectures in equivariant settings.

Quality:
The paper demonstrates high quality in both theoretical and empirical aspects. It provides a solid theoretical foundation based on Markov kernels, offering a clear understanding of how SYMDIFF achieves equivariance without specialized architectures. The empirical evaluations are thorough, showcasing the method’s effectiveness on E(3)-equivariant molecular generation tasks using the well-known QM9 and GEOM-DRUGS datasets. The results indicate significant performance improvements and scalability advantages over traditional equivariant models. The experiments are well-designed, and the conclusions are well-supported by the data.

Clarity:
The paper is written with clarity and precision, making complex concepts accessible to the reader. It systematically introduces the theoretical background before delving into the methodological details of SYMDIFF. The explanations are clear, and algorithmic descriptions helps in understanding the implementation. The empirical results are presented in an organized manner, with insightful discussions that highlight the implications of the findings. Overall, the paper effectively communicates its contributions and significance.

Significance:
The significance of the paper lies in its potential to impact the field of equivariant machine learning models substantially. By eliminating the need for specialized equivariant architectures, SYMDIFF lowers the barrier to entry for researchers and practitioners interested in building equivariant models. It enables the use of arbitrary off-the-shelf neural network architectures, which are often more scalable and easier to implement. This flexibility can accelerate the development and deployment of equivariant models in various applications, including molecular generation, computer vision, and physics-based simulations. Moreover, being the first to apply stochastic symmetrisation to generative modeling, the paper paves the way for future research in this promising direction.
summary: The paper “SymDiff: Equivariant Diffusion via Stochastic Symmetrisation” introduces SYMDIFF, a novel method for constructing equivariant diffusion models without the need for complex, intrinsically equivariant neural networks. Traditional approaches often rely on specialized architectures with complex parameterizations and higher-order geometric features to achieve equivariance, which can be computationally intensive and challenging to implement.

SYMDIFF leverages the framework of stochastic symmetrisation, effectively functioning as a learned data augmentation technique applied during the sampling phase. This approach is lightweight, computationally efficient, and can be easily implemented on top of arbitrary off-the-shelf models. By not requiring the neural network components to be intrinsically equivariant, SYMDIFF allows the use of highly scalable modern architectures, providing greater flexibility and performance.

The authors demonstrate the effectiveness and scalability of SYMDIFF through empirical experiments on E(3)-equivariant molecular generation tasks using the well-known QM9 and GEOM-DRUGS datasets. The results show significant improvements over traditional methods, highlighting the benefits of SYMDIFF in generating molecular structures that are equivariant under the E(3) group (which includes translations, rotations, and reflections in three-dimensional space).

This work represents the first application of stochastic symmetrisation to generative modeling, suggesting that SYMDIFF has the potential for broader applications in the domain of equivariant machine learning models.
weaknesses: While the paper introduces SYMDIFF as a flexible and efficient method for constructing equivariant diffusion models, it lacks specific guidelines to help practitioners decide when to use SYMDIFF over traditional intrinsically equivariant models. This absence of decision-making support makes it challenging to understand the circumstances under which SYMDIFF is more advantageous.

To address this, I suggest the authors include a comparative analysis between SYMDIFF and traditional methods. This could be presented as a table or flowchart highlighting key factors such as:

- Dataset Characteristics: Situations where SYMDIFF performs better or worse.
- Model Complexity: Differences in architecture and scalability.
- Computational Resources: Comparisons of memory usage, training time, and inference speed.
- Performance Metrics: Quantitative comparisons on accuracy and sample quality.
- Ease of Implementation: Practical considerations like code complexity.

Including such a comparison would provide clear guidance for practitioners, making the paper more actionable and useful.

Additionally, while the paper claims that SYMDIFF is computationally efficient, it doesn’t provide empirical evidence to support this. I recommend the authors include:

- Empirical Benchmarks: Runtime analysis and resource consumption statistics compared to intrinsically equivariant models.
- Convergence Plots: Performance versus number of training epochs for both SYMDIFF and baseline methods.
- Ablation Studies: Impact of stochastic symmetrisation on computational overhead and model performance.

By incorporating these elements, the paper would strengthen its claims and offer actionable insights, enabling readers to make informed decisions about adopting SYMDIFF.

---

confidence: 5
contribution: 3
initial_rating: 5
presentation: 3
questions: - Why SymDiff require larger model size to match the performance of EDM / GeoLDM etc.
soundness: 3
strengths: - Originality: The paper proposes SymDiff, a novel method using stochastic symmetrisation for equivariant diffusion models, extending symmetrisation to generative modelling and differing from prior approaches.

- Clarity: The paper is well-structured, clearly explaining concepts with examples and presenting algorithms, aiding understanding and implementation.

- Quality: The paper presents a comprehensive framework with solid derivations and conducts rigorous experiments on multiple datasets.
summary: This paper proposes SymDiff, a novel method for constructing equivariant diffusion models via stochastic symmetrisation. It addresses limitations of prior work on geometric data by allowing more flexibility than using intrinsically equivariant neural networks. Experiments on QM9 and GEOM-Drugs datasets for molecular generation show SymDiff outperforms EDM and is competitive with other baselines, with better scalability. SymDiff can be seen as a generalisation of data augmentation, with a learned and deployed augmentation process ensuring stochastic equivariance.
weaknesses: - Performance gap: from Tab2 and 3, seems the model is not parameter efficient compared with baseline models. When reducing the model size, the performance will have a larger drop than the baseline model.

- Unfair comparison: I didn't fully understand for SymDiff+++, ++, and +, why you cannot choose the same number of parameters of EDM. It seems your model SymDiff+ is already larger than EDM+++, but with worse performance.

- Missing comparison: I think the author also needs to incorporate GeoLDM with different sizes into Tab 2 and 3, which potentially can work better since it acts on low-dimensional space.

Meta-review: This paper introduces SymDiff, a novel method for constructing equivariant diffusion models through stochastic symmetrization. Unlike prior approaches that rely on intrinsically equivariant neural networks, SymDiff offers greater flexibility for geometric data. Experiments on the QM9 and GEOM-Drugs datasets demonstrate that SymDiff outperforms EDM and achieves competitive performance with other baselines, while offering improved scalability. Conceptually, SymDiff generalizes data augmentation by learning and applying a stochastic augmentation process that enforces equivariance.

All reviewers agree that the work is solid and interesting.

Decision: Accept (Poster)