1 Introduction↩︎

The application of machine learning (ML) to the development of surrogate models of complex physical systems is an active area of research. This is due to the intrinsic nature of physical systems which are often so complex that conventional surrogate models fail to describe them adequately. Prolific research areas for machine learning surrogates include quantum-accurate force fields for molecular dynamics [1], coarse-grained molecular dynamics [2], and more generally of learning dynamical systems from data [3].

The focus of the present work is on the role of symmetries in machine learning model architectures. These can be found everywhere in the physical, chemical and biological world, be it in atomic potentials, molecular structures or protein interactions. Two competing approaches to handling symmetry in machine learning models are (i) using a generic model that is unaware of symmetry but augmenting training data during learning [4], [5]; (ii) exactly enforcing known symmetry in the model [6], [7]. Each possibility can be achieved through a number of available techniques. It remains an ongoing debate (cf. foregoing references) which of these two approaches are preferable, and how they are to be optimally implemented. Generally speaking, striking a balance between enforcing structure in model-design and learning it from data remains a challenge in scientific machine learning. The purpose of our work is to rigorously analyze two natural data augmentation approaches to learn a known model symmetry and to contrast them against each other as well as against an analogous model with built-in symmetry.

An argument often advanced in favour of hard-coding symmetry into the model architecture is that it reduces the number of model parameters. This is generally true, but in the context of deep learning, the number of parameters is rarely an adequate measure of a model’s performance or robustness. The properties of training dynamics is generally more important and some evidence is accumulating that symmetry-unconstrained model may indeed be easier to train and are equally robust and accurate out-of-distribution [8], [9]. This has already lead to the development of highly efficient unconstrained models for static tasks such as geometric optimization [10], [11]

Our concern with preserving symmetry instead focuses on conserving important qualitative physics of a model. According to Noether’s theorem, symmetries and conserved quantities are intrinsically linked; the presence of a continuous symmetry is equivalent to the existence of a first integral of the motion. For example, translation invariance results in conservation of linear momentum while rotation invariance results in conservation of angular momentum. For many applications, it is desired that these conserved quantities are preserved by the numerical scheme. This led, for example, to the development of geometric integration, in particular symplectic integration [12]. A symplectic integrator is an integrator such that its time-stepping operation is a symplectomorphism, a consequence being that any first integrals of the simulated system will be preserved [13] for a separable Hamiltonian.

In section 2, we will show that this the impact of having only approximate symmetry in a surrogate model generically has a substantial impact on the conservation of the first integrals. We are particularly interested in the case of invariance under \(\mathrm{SO}(n)\) whose associated conserved quantity in two and three dimensions is angular momentum. We review established theory explaining that — except in the case of fully integrable systems — one will observe an exponential drift in the time evolution of these conserved quantities. This motivates the need for surrogate models models to preserve conserved quantities exactly or to within very high precision.

In section 3 we then present and study two different methods of data augmentation in order to improve the symmetry properties of a symmetry-unconstrained approximation scheme. These are augmenting the data through a quadrature rule for the Haar measure on the compact Lie group whose representation is the continuous symmetry of interest and a random sampling of that Haar measure. We rigorously prove convergence (rates) for both of these schemes. Finally, in section 4, we provide a more empirical analysis to illustrate our theory and provide further details of the performance of data augmentation schemes approximation schemes through a Least-Squares toy model - we will be considering the simplest non-trivial case of three particles on \(S^1\) and \(S^2\). In section 3 we then present and study two different methods of data augmentation in order to improve the symmetry properties of a symmetry-unconstrained approximation scheme. These are augmenting the data through a quadrature rule for the Haar measure on the compact Lie group whose representation is the continuous symmetry of interest and a random sampling of that Haar measure. We rigorously prove convergence (rates) for both of these schemes.

Finally, in section 4, we provide empirical numerical tests to illustrate our rigorous theory and to provide further details of the performance of the data augmentation schemes.

2 Approximate conservation of angular momentum↩︎

In this section, we review results concerning approximate conservation properties of symplectic integrators and their breakdown when a system is only approximately symmetric. \(H=T+V\) denotes a separable Hamiltonian and \((p,q)\) the canonical coordinates in phase space. This section gives a formal discussion of key ideas; we refer to [12] for technical details.

A numerical integrator for a Hamiltonian system is symplectic if it is compatible with the manifold structure; cf. [12], [14]. The canonical example of a symplectic integrator is the order 2 Velocity Verlet integrator we chose for our numerical applications [15]. While the dynamics obtained through the integration scheme are not exact, they correspond to a Hamiltonian that is close to the original one, known as the shadow Hamiltonian \(H_{\mathrm{eff}}\). For any symplectic integrator of order \(p\) with stepsize \(h\) \[H_{\mathrm{eff}}=H+O(h^p).\]

One of \(H_{\mathrm{eff}}\)’s constructions may be found in [13]. It implies that any quantity that commutes with both \(T\) and \(V\) also commutes with \(H_{\mathrm{eff}}\) as all its higher order terms are commutators of \(T\) and \(V\). In practice this is a strong result on the conservation of first integrals, as \(J\) is a first integral iff \(\{J,H\}=0\). Thus, \(H\) and \(H_{\mathrm{eff}}\) share the same conserved quantities.

We are interested in the case of a rotation-invariant potential \(V\), resulting in conservation of angular momentum under Hamiltonian dynamics. More generally, we can consider continuous symmetries.

Suppose now, that during an approximation procedure (see later section) exact invariance is violated, but the resulting Hamiltonian remains approximately invariant, that is, \(H=H_{\mathrm{I}}+\varepsilon H_{\mathrm{P}}\), where \(H_{\mathrm{I}}\) has a continuous symmetry, \(\varepsilon > 0\) and \(H_{\mathrm{P}}\) is a Hamiltonian. In this case, approximate conservation properties are not guaranteed. Two scenarios best illustrate the expected behavior: either trajectories stay controlled or the system contains some amount of chaos. The general setting belongs to the latter case; even under analytic assumptions, we observe an exponential growth in the error. This is expected, as it is well-known that if a system is not fully integrable, then it is likely for chaos to be present [16]–[20].

Although symplectic integrators do not alter this behaviour, they do guarantee that conserved quantities are preserved correctly. To observe this, we substitute \(X_{\mathrm{V}} = X_{\mathrm{i}}+\epsilon X_\mathrm{p}\) into \(X_{\mathrm{eff}}\)’s expression, where \(X_{\mathrm{i}}\) respectively \(X_{\mathrm{p}}\) are the Liouville operators [13], [14] for \(H_{\mathrm{i}}\) respectively \(H_{\mathrm{p}}\). This yields

\[\begin{align} X_{\mathrm{eff}} &= X_{\mathrm{T}}+ X_{\mathrm{i}}+\epsilon X_{\mathrm{p}} + \frac{h}{2}[X_{\mathrm{T}},X_{\mathrm{i}}]+\frac{h\epsilon}{2}[X_{\mathrm{T}},X_{\mathrm{p}}]+\cdots\\ J\circ \Phi_t&= e^{tX_{\mathrm{eff}}}J\\ &= e^{t\epsilon X_{\mathrm{eff,p}}}e^{t(X_{\mathrm{T}}+X_{\mathrm{i}})}J\\ &=e^{t\epsilon X_{\mathrm{eff,p}}}J, \end{align}\] where \(\Phi_t\) is the Hamiltonian flow, that is the map such that \(\Phi_t(p(0),q(0))=(p(t),q(t))\). We remark that \(e^{t\epsilon X_{\mathrm{eff,p}}}J=J\) iff \(X_{\mathrm{p}}\) commutes with \(J\).

While exponential deviation is the generic situation, more rigid behaviour appears in integrable systems, those with \(n\) (the dimension of configuration space) functionally independent constants of the motion. KAM theory [21] shows that in a completely integrable system, sufficiently irrational non-degenerate tori survive perturbation, implying that conserved quantities will be approximately conserved in the perturbed system. This is shown in Chapter XI of [12].

We illustrate this behaviour through the following simulations. We consider 3 particles on the circle \(S^1\) under a potential \(U=F+\epsilon G\), where \(F\) is rotation-invariant and analytic, while \(G\) is chosen to be analytic but not invariant under rotations. Due to Noether’s theorem, a continuous symmetry in the potential is equivalent to the presence of a first integral; rotation invariance corresponds to the conservation of angular momentum. The initial conditions are such that the initial total angular momentum vanishes. We perform the simulations for different values of \(\epsilon\). To compare approximately symmetric potentials, we define the symmetrization error.

When \(f\) is \(G\)-invariant \(\epsilon_{\mathrm{sym}}(f)=0\). The choice of norm can in general have strong consequences. However, in our applications we consider only very smooth function classes such that the \(L^2\) and \(L^{\infty}\) norms remain comparable up to moderate constants.

Figure \(\ref{fig:ApproxJ2d}\) illustrates the rapid breakdown of the conservation of angular momentum for approximately rotation-invariant potentials. Thus, if one wishes to preserve conserved quantities, either one should be working with a fully integrable system or have an approximation scheme that exactly preserves symmetries.

To conclude, we motivated that except in the very fringe case of full integrability, the Hamiltonian one gives to their numerical scheme should preserves symmetries as exactly as possible. Indeed, outside of KAM theory, only an exactly symmetric Hamiltonian will provably always have simulated dynamics whose outputs preserve the desired first integrals.

2.1 Geometry Optimization↩︎

Our foregoing discussion makes an argument that symmetry must be preserved to within very high precision in order to prevent loss of important conservation laws and the resulting (potentially catastrophic) failure of molecular dynamics simulations. Whether such failure occurs would depend on the specific simulation undertaken, e.g., the importance of conservation of angular momentum is different in a gas, liquid or solid.

In addition, this does not prevent a non-symmetry preserving model from being well suited for other simulation tasks. For example, in “geometry optimization” (computation of energy minima, saddle points, etc) we generally seek an equilibrium point of \(V\), i.e., solve the nonlinear system \[\nabla V(\bar{q}) = 0.\] Suppose now that an approximate model \(V_\epsilon\) exists such that \(|V - V_\epsilon| \leq \epsilon^2\), and \(\|\nabla V - \nabla V_\epsilon \|_2 = \epsilon\) in a neighbourhood of such an equilibrium point \(\bar{q}\). Under suitable and natural additional assumptions on the stability of the equilibrium point, the inverse function theorem guarantees existence of an approximate equilibrium \(\bar{q}_\epsilon\) such that \[\nabla V_\epsilon(\bar{q}_\epsilon) = 0, \qquad \| \bar{q} - \bar{q}_\epsilon \|_2 \leq C \epsilon \quad \text{and} \quad |V(\bar{q}) - V_\epsilon(\bar{q}_\epsilon)| \leq C \epsilon^2.\] Obtaining this conclusion required no information about the symmetry of \(V_\epsilon\). Thus, the requirement to exactly (or, to within high precision) preserve symmetries is necessarily highly task dependent. This is consistent with numerous empirical observations on non-trivial benchmark tasks [11].

3 Analysis of Data Augmentation↩︎

We first summarize the notation used throughout the rest of the paper. \(\mathbb{N}\) are the natural numbers including 0. \(S^{d},d\in\mathbb{N}\setminus\{0\}\) is the \(d\)-sphere and \(S^{N,d}\) its \(N\)-fold product. We denote by \(L^2(S^{N,d})\) the space of complex-valued, measurable, square integrable functions. \(R = (\mathbf{r}_j)_{j=1}^N\in S^{N,d}\) is a collection of points on the general \(d\)-sphere, for the special case of \(d=1\), we denote by \(\Theta=(\theta_1,\cdots,\theta_N)\in S^{N,1}\) a collection of angles representing points on the unit circle. The symbols \(\mathbf{l},\mathbf{m},\mathbf{k}\) denote multi-indices, that is, tuples of integers, for instance \(\mathbf{k}=(k_1,\cdots,k_N)\in\mathbb{Z}^N\). \(Q\in \mathrm{SO}(d+1)\) is a rotation. The action of \(Q\) on \(\Theta\) or \(R\) is written as \(Q\cdot \Theta\) or \(Q \cdot R\) and stands for element-wise rotation, i.e., \[Q\cdot R = (Q\mathbf{r}_1,\cdots, Q\mathbf{r}_N). \label{eq:defAction}\tag{1}\]

3.1 Polynomial Regression↩︎

We study different ways of approximating (or, “learning”) a rotation-invariant potential \(f\colon S^{N,d}\to\mathbb{C}\). We will focus on approximations using polynomial tensor product bases. Their elements are, \[\begin{align} \varphi_{\mathbf{k}}(\theta_1,\cdots, \theta_N)&=\prod_{j=1}^N e^{ik_j\theta_j}, \qquad \text{ for d=1}, \\ \varphi_{\mathbf{l,m}}(\mathbf{r}_1,\cdots, \mathbf{r}_N)&=\prod_{i=1}^N Y^{m_i}_{l_i}(\mathbf{r}_i), \qquad \text{ for d=2}, \label{eq:tensorbasis} \end{align}\tag{2}\] where the \(Y^{m}_{l}\) are complex spherical harmonics: for \(l\in\mathbb{N}\), \(-l\leq m\leq l\), \((\theta,\varphi)\in S^2\), they are defined by \[\begin{align} Y^m_l(\theta,\varphi) &= (-1)^m\sqrt{\frac{(2l+1)(l-m)!}{4\pi(l+m)!}}P^m_l(\cos(\theta))e^{im\varphi}, \\ P^m_l(x) &= \frac{(-1)^m}{2^l l!}(1-x^2)^{\frac{m}{2}}\frac{d^{l+m}}{dx^{l+m}}(x^2-1)^l. \end{align}\] In practice, this definition is never used as we rely on numerically stable recursive evaluation algorithms [22]. These functions are polynomials in the cartesian coordinates of \(\mathbb{R}^{d+1}\). We denote by \[V^{N,d}_{K}\subset L^2(S^{N,d}) \label{eq:defV}\tag{3}\] the subspace of polynomial functions (in cartesian coordinates) on \(S^{N,d}\) of total degree at most \(K\); specifically, \[\begin{align} V_K^{N,1} &:= \mathop{\mathrm{Span}}\left\{ \varphi_{\mathbf{k}}\;\colon \; \mathbf{k}\in \mathbb{Z}^N, \norm{\mathbf{k}}_1\leq K \right\}, \, \quad \text{and} \\ V^{N,2}_{K} &:= \text{Span}\left\{\varphi_{\mathbf{l},\mathbf{m}}\colon {\boldsymbol{l}} \in \mathbb{N}^N, {\boldsymbol{m}} \in \mathbb{Z}^N, |m_i| \leq l_i, \| {\boldsymbol{l}} \|_1 \leq K \right\}. \end{align}\] The total degree of \(\varphi_{\boldsymbol{l} \boldsymbol{m}}\) is taken to be \(\| {\boldsymbol{l}}\|_1\) to ensure that \(V^{N,2}_{K}\) is closed under rotations ([23]; see also § 3.2). It is well-known that \(\bigcup_{K} V^{N,d}_K\) - the space of all polynomials - is dense in \(L^2\) [24].

For notational convenience we specify an arbitrary ordering, \({\boldsymbol{k}}_j\) or \(({\boldsymbol{l}}_j, {\boldsymbol{m}}_j)\) of the basis elements and identifying \(\varphi_j = \varphi_{{\boldsymbol{k}}_j}\) resp. \(\varphi_j = \varphi_{{\boldsymbol{l}}_j {\boldsymbol{m}}_j}\). An arbitrary polynomial \(P \in V_K^{N, d}\) can then be written in the form \[P_{\beta} = \sum_{j}\beta_{j}\varphi_{j},\] with parameters \(\beta_{j}\).

We employ a Least-Squares Approximation scheme (LSQ) to construct approximants. Let \((R_i, f(R_i))_{i=1}^n\) be iid data sampled with respect to some distribution (they are specified in Section 4). From it, we construct the design matrix \(A\), \(A_{ij}=(\varphi_{j}(R_i))\), of degree \(K\) and the data vector \(Y\), with entries \(Y_i = f(R_i)\). The LSQ parameters are then determined by solving the linear least squares problem \[\beta = \mathrm{argmin}_{\beta}L(\beta), \qquad \text{where} \qquad L(\beta)=\frac{1}{2}\norm{A\beta - Y}^2. \label{eq:fullLSQ}\tag{4}\] It is well known that under the assumption that there is enough data, the LSQ approximates the \(L^2\) projection onto \(V^{N,d}_K\) [25]. However, no additional information is given about the symmetry of the estimated polynomial, and it is a priori unclear whether the deviation from exact invariance is such that it could lead to dramatic qualitative failures of a simulation (cf. Figure 1).

3.2 Regression in Invariant Subspaces↩︎

Arguably the canonical approach to learning an invariant function is restricting the approximation to the subspace of rotation-invariant polynomials on \(S^{N,d}\) of total degree at most \(K\): \[B^{N,d}_K := \big\{P\in V^{N,d}_K\colon P\circ Q=P,\forall Q\in \mathrm{SO}(d+1) \big\} \subset V^{N,d}_K,\] where \(V^{N,d}_K\), defined in 3 , is the space of polynomials of total degree at most \(K\) on \(S^{N,d}\) and the rotations have elementwise 1 action.

The proof of the Lemma is standard; we include it for the reader’s convenience.

Using Lemma 1 we can characterise the rotation-invariant subspaces. For \(d=1\), \[\begin{align} B_K^{N,1} & = \mathrm{Span}\left\{\varphi_{\mathbf{k}} \;\colon\; \varphi_{\mathbf{k}}\in V^{N,1}_K,\sum_{i=1}^N k_i=0\right\}. \end{align}\] This can be shown by applying \(\mathrm{sym}\) to the basis elements of \(V^{N,1}_K\). \[\begin{align} \notag \mathrm{sym}(\varphi_{\mathbf{k}}) &= \frac{1}{2\pi}\int_{0}^{2\pi}\prod_{j=1}^Ne^{ik_j(\theta_j+\alpha)}d\alpha \\ &=\frac{1}{2\pi}\int_{0}^{2\pi}e^{i\alpha\sum_j k_j}\prod_{j=1}^Ne^{ik_j\theta_j}d\alpha = \begin{cases} \varphi_{\boldsymbol{k}}, & \sum_i k_i = 0, \\ 0, & \text{otherwise.} \end{cases}. \end{align}\]

There are multiple natural choices of a basis for the rotation-invariant subspace of \(V^{N,2}_K\). Possible approaches include applying ?? directly, or to recursively decompose tensor products of pairs of irreducible representations using Clebsch–Gordan coefficients. Any approach results in invariant basis functions \[\varphi_{\mathbf{l},\mu}=\sum_{\mathbf{m\in\mathcal{M}_{\mathbf{l}}}}\mathcal{C}^{\mathbf{l}\mu}_{\mathbf{m}}\varphi_{\mathbf{l}\mathbf{m}},\] where \(\mu\) encodes all the possible invariant couplings and \(\mathcal{C}^{\mathbf{l}\mu}_{\mathbf{m}}\) are generalized (higher-order) Clebsch-Gordan coefficients. The details can be found in [23] and [26].

All our numerical tests are performed with \(N=3\) particles, in which case a basis of \(B^{3,2}_L\) is given by \[\begin{align} \tilde{\varphi}_{\mathbf{l}}(R) = \sum_{\substack{\boldsymbol{\mu}\in\mathcal{M}_{\mathbf{l}}\\ \mu_1+\mu_2+\mu_3=0}} C^{l_3,0}_{l_1,0,l_2,0}C^{l_3,\mu_1+\mu_2}_{l_1,\mu_1,l_2,\mu_2}C^{00}_{l_3,0,l_3,0} C^{00}_{l_3,-\mu_3,l_3,\mu_3}\prod_{i=1}^3 Y^{m_i}_{l_i}(\mathbf{r}_i), & \\ \text{where } \quad |l_1-l_2|\leq l_3 \leq l_1+l_2.& \label{eq:3dIbasis} \end{align}\tag{5}\]

We can now define the invariant design matrix, \((A_{\mathrm{I}})_{ij}=\tilde{\varphi}_{j}(x_i)\) where we have similarly specified an ordering of the basis elements of \(B^{N,d}_K\). The invariant regression problem is \[\bar\beta_{\mathrm{I}}= \mathrm{argmin}_{\beta_{\mathrm{I}}} L_{\mathrm{I}}(\beta_{\mathrm{I}}) \quad \text{where} \quad L_{\mathrm{I}}(\beta_{\mathrm{I}})=\frac{1}{2}\norm{A_{\mathrm{I}}\beta_{\mathrm{I}}-Y}^2 \; . \label{eq:invlsq}\tag{6}\]

This idea of enforcing structure in the architectures is widespread in the ML community. For instance, Physics Informed Neural Networks [27] are neural networks that force their outputs to be solutions to the equations they are trying to learn. In the case of symmetries, enforcing such a constraint in a ML model leads to the so called Equivariant Neural Networks (ENNs) [28]. These are widely used to study atomic potentials [6], [26] or other domains where symmetries are widespread such as lattice gauge theory [29]. Applications are not solely restricted to physics, for example, ENNs are also used in image processing [30] as two images that differ by some rigid body transformation are usually thought to be the same.

3.3 Data Augmentation↩︎

Enforcing symmetry will always result in exact symmetry of the learned model. Interestingly, however, it may not always be the optimal choice. In fact, recent studies have shown that enforcing symmetry may hurt performance [31] or make training harder [9]. To avoid these possible issues, one may want to chose another way of learning symmetry, for instance through data augmentation. For \(T\in\mathbb{N}\), let \((w_t, Q_t)_{t=1}^T\) be a collection of real weights and rotations. We will specify below how they are obtained. Then the symmetry-augmented linear least squares problem is \[\beta_{\rm aug} = \arg\min L_{\rm aug}(\beta) \qquad \text{where} \qquad L_{\mathrm{aug}}(\beta) =\frac{1}{2} \sum_{t=1}^Tw_t \big\|(A\circ Q_t)\beta-Y \big\|^2, \label{eq:lsqaug}\tag{7}\] where \(A\circ Q_t\) denotes the design matrix evaluated at the points \(Q_t\cdot R\), \(t = 1, \dots, n\). The notation \(A\circ Q\) comes from the fact that we may think of the design matrix as a function from data space to the space of matrices. The design matrix obtained from rotating data by \(Q_t\) then corresponds to a composition with \(Q_t\). Remark that \(Y\) does not need to be rotated as the target-function is chosen to be rotation-invariant.

This idea stems from the above observation that for a rotation-invariant function \(f\), its value at one point \(R\) is constant on the whole orbit \(G\cdot R\). We present two ways of augmenting data: a random and an optimal sampling of \(G\), leading to a near square root decay in \(\varepsilon_{\mathrm{sym}}\) respectively exact invariance. An important observation is that this type of augmentation does not correspond to a simple symmetrisation of the solution. The dependence of the output of our augmented LSQ system (cf 7 ) is nonlinear and the effects of augmentation must carefully analysed.

It is well known that \(V^{N,d}_K\) is closed under rotations. That is, for any \(v\in V^{N,d}_K\) and \(Q\in \mathrm{SO}(d+1)\), \(v\circ Q\in V^{N,d}_K\). Moreover, this action is linear on the space of polynomials. They are related with each other through the (unitary) Wigner \(D\)-matrices [23] denoted by \(D\). If we consider the spherical harmonic \(Y^m_l\), then its rotation by \(Q\) may be written as \[Y^m_l\circ Q = \sum_{m'=-l}^lY^{m'}_l D^l_{m'm}(Q).\]

These matrices \(\mathcal{D}\) can be thought of as generalised \(D\)-matrices. We refer to [23], [26] for details, but we include an explicit computation for the case of particles on the circle (\(d=1\)) for the sake of illustration and intuition. We recall that a basis element \(\varphi_{\mathbf{k}}\) is given by a product of trigonometric polynomials 2 . Acting on it by a rotation \(Q_{\alpha}\) of angle \(\alpha\), we obtain \[\begin{align} Q_{\alpha}\cdot \varphi_{\mathbf{k}}(\Theta)&=\prod_{j=1}^N e^{ik_j(\theta_j+\alpha)}\\ &=e^{i\alpha\sum_jk_j}\phi_{\mathbf{k}}(\Theta). \end{align}\] Thus, after reordering the basis to achieve the decomposition \(V^{N,1}_K=B^{N,1}_K\oplus (B^{N,1}_K)^{\perp}\), \(\mathcal{D}(Q_{\alpha})\) has block diagonal form: \[\mathcal{D}(Q_{\alpha})=\begin{pmatrix} \mathrm{Id}& 0\\ 0 & \mathrm{D}_{\alpha, \mathrm{N}} \end{pmatrix}, \label{eq:DQd611}\tag{8}\] where \(\mathrm{D}_{\alpha, \mathrm{N}}\) is a diagonal matrix with entries \(e^{i\alpha\sum_jk_j}\) on the row corresponding to \(\varphi_{\mathbf{k}}\).

For notational simplicity we introduce \(D_{t}:= \mathcal{D}(Q_t)\), then \(\nabla_{\beta}L_{\mathrm{aug}}(\beta)\) can be readily expressed as \[\begin{align} \nabla_{\beta}L_{\mathrm{aug}}(\beta) &=\sum_t w_tD_{t}^*(A^*AD_{t}\beta-A^*Y), \end{align}\] resulting in the normal equations \[\begin{align} \label{eq:normaleqns} \bigg(\sum_t w_tD_{t}^* A^*AD_{t}\bigg) \beta_{\rm aug} =\sum_t w_tD_{t}A^*Y. \end{align}\tag{9}\]

To understand the effect of data augmentation we need to analyze 9 , a task make non-trivial by the fact that the \(D_{t}\) matrices occur nonlinearly. We begin by decomposing \(V^{N,d}_K=B^{N,d}_K\oplus (B^{N,d}_K)^{\perp}\) as indicated above for \(d = 1\). We can select a new orthonormal basis that is compatible with this decomposition. With abuse of notation, this results in new design matrices \(A\), parameters \(\beta\) and \(D\)-matrices \(D_{t}\) that have block structure with \({\rm I, N}\), respectively, denoting the invariant and non-invariant parts: \[A=\begin{pmatrix} A_{\mathrm{I}}&A_{\mathrm{N}} \end{pmatrix},\; A^*=\begin{pmatrix} A_{\mathrm{I}}^*\\ A_{\mathrm{N}}^* \end{pmatrix},\; D_{t}=\begin{pmatrix} \mathrm{Id}&0\\ 0& D_{t,\mathrm{N}} \end{pmatrix},\; \beta = \begin{pmatrix} \beta_{\mathrm{I}}\\ \beta_{\mathrm{N}} \end{pmatrix}.\] In this decomposition, the shape of \(D_{t}\) follows as rotations act trivially on \(B^{N,d}_K\) which implies that \((B^{N,d}_K)^{\perp}\) must also be closed under this action. Indeed, acting by a rotation and its inverse \(QQ^{-1}\) should be the trivial action.

For a given \(t\), we obtain \[\begin{align} D_{t}^* A^*AD_{t}&= \begin{pmatrix} A_{\mathrm{I}}^*A_{\mathrm{I}}&A_{\mathrm{I}}^*A_{\mathrm{N}}D_{t,\mathrm{N}}\label{Eq:Augdm}\\ D_{t,\mathrm{N}}^*A_{\mathrm{N}}^*A_{\mathrm{I}}&D_{t,\mathrm{N}}^*A_{\mathrm{N}}^*A_{\mathrm{N}}D_{t,\mathrm{N}} \end{pmatrix}. \end{align}\tag{10}\] Lastly, using the fact that the weights of a quadrature rule sum to the measure of the integration space we obtain \[\sum_t w_tD_{t}^* A^*AD_{t}= \begin{pmatrix} A_{\mathrm{I}}^*A_{\mathrm{I}}& \sum_i w_t A_{\mathrm{I}}^*A_{\mathrm{N}}D_{t,\mathrm{N}}\\ \sum w_tD_{t,\mathrm{N}}^*A_{\mathrm{N}}^*A_{\mathrm{I}}&\sum_i w_tD_{t,\mathrm{N}}^*A_{\mathrm{N}}^*A_{\mathrm{N}}D_{t,\mathrm{N}} \end{pmatrix}. \label{eq:EzAug}\tag{11}\]

In the following computations, we will write 11 in block form \[\begin{pmatrix} B&C\\ C^* & D \end{pmatrix},\] and assume that the original, un-augmented, design matrix \(A\) has full column rank. In our numerical applications, \(A\) being of full rank is equivalent to having more data points than parameters.

This yields the existence of constants \(c_0,c_1>0\) such that for any vector \(u \in \mathbb{R}^{n_{\rm N}}\) (\(n_{\mathrm{N}}\) being the number of columns of \(A_{\mathrm{N}}\)) and any rotation \(Q\), \[\begin{align} c_0\norm{u}&\leq \norm{A_{\mathrm{N}}u}\leq c_1\norm{u} \qquad \text{and} \\ c_0\norm{u}&\leq \norm{A_{\mathrm{N}}D(Q)u}\leq c_1\norm{u}, \label{eq:LAbounds} \end{align}\tag{12}\] where the second line is a consequence of the unitarity of \(D(Q)\).

Combining Lemma 4 and Remark 1 results in the following upper bound: \[\varepsilon_{\mathrm{sym}}(P_{\mathrm{A}})\leq \frac{1}{c_2} \norm{A_{\mathrm{N}}\overline{D}_T}_{\mathrm{op}}\norm{A_{\mathrm{I}}\bar\beta_{\mathrm{I}}-Y}, \label{eq:SchurUB}\tag{13}\] \(\norm{\cdot}_{\mathrm{op}}\) being the operator norm induced by the 2-norm.

3.3.1 Quadrature Data Augmentation↩︎

From 11 we observe that if \((w_t,Q_t)_{t=1}^T\) is a quadrature rule for \(\mathrm{SO}(d+1)\), then the augmented system becomes block diagonal as quadrature rules integrate polynomials of lower degrees exactly.

This is a standard fact; the \(\mathcal{D}\) are finite dimensional representations of \(\mathrm{SO}(d+1)\) [33].

In order to benefit from the quadrature augmentation method one has to obtain the appropriate quadrature rules. For \(\mathrm{SO}(2)\), after identifying it with \([0,2\pi)\), we obtain that an equi-weight quadrature rule with degree of accuracy \(n-1\) is given by \(\big\{\frac{1}{n}, \frac{2\pi t}{n}\colon t\in[n-1]\big\}\) [34]. For \(\mathrm{SO}(3)\), quadrature rules have been developed in [35], which we use for all our numerical tests. These rules are significantly more computationally expensive both to obtain and to use.

We illustrate this with Table 1 where we provide the size of some quadrature rules [35], [36]. The rapid increase in the size of quadrature rules cannot be changed.

3.3.2 Random Data Augmentation↩︎

The common way of approximating this symmetrisation process is by augmenting via randomly sampled rotations [37], [38]. We show that this leads to a near \(T^{-1/2}\) decay in \(\varepsilon_{\mathrm{sym}}\).

While random augmentations are easy to implement, it has an asymptotic error rate lower bounded by \(T^{-1/2}\), which, given the discussion in Section 2, is insufficient.

4 Numerical Tests↩︎

We illustrate the theory developed in the previous section with a range of numerical tests. For these tests, we choose \(d \in \{2, 3\}\) and \(N=3\), i.e. systems consisting of three two- or three-dimensional particles. The choice \(N = 3\) corresponds to the simplest non-trivial invariant basis. Indeed, for \(N=1\), the only rotation-invariant functions are the constant functions, while for \(N=2\), the position of any one particle determines the entire system. Only from \(N=3\) onwards do we observe non-trivial couplings.

Restricting our tests to relatively few particles is also motivated by the fact that in much larger point cloud models, dimension reduction techniques are always applied to reduce the many particle interaction to moderate correlation order. For example, the highly successful MACE model [6] employs just two layers of \(N=3\) for modelling complex interatomic interactions.

An additional common choice in molecular dynamics especially for materials systems is that particles are rotated into a reference frame (natural unit cells or supercells). We therefore also consider non uniform distributions. However, these distributions have the property that when they are averaged over \(\mathrm{SO}(d+1)\) they become the uniform distribution. Except when otherwise stated, we only vary the distributions when fitting the model, the error will be computed using the uniform distribution. Our models are not permutation invariant, thus when sampling the non-uniform distributions, the choice of the constrained particles is fixed throughout the tests. It would have been possible to make our bases permutation invariant, but we chose not to as we wish to focus on the specific case of rotation invariance.

4.1 Results in Two Dimensions↩︎

The distributions used for two-dimensional tests are summarized in Table 2 and illustrated in Figure 2.

As target function we choose a high-degree polynomial \(f \in B_{30}^{3,1}\), i.e., such that its coefficients have an exponential decay of rate \(\alpha=2\), \[f(\Theta) = \sum_{\mathbf{k}\in\mathbb{Z}^N,\norm{\mathbf{k}}_{1}\leq \tilde{K}}c_\mathbf{k} e^{-\alpha||\mathbf{k}||_1} e^{\mathbf{k}\cdot \Theta},\] with \(c_\mathbf{k}\) drawn uniformly from \([-1,1]\). The exponential decay is a natural generalization of analyticity in the univariate setting. Due to the assumed decay, it is straightforward to see that there exists \(C>0\) such that \[||f-f_K||_{L_2} \leq Ce^{-\alpha K},\] where \(f_K \in B^{3,1}_K\) is the degree \(K\) truncation (\(L^2\)-orthogonal projection) of \(f\). For sufficient training, all LSQ methods should attain this approximation error rate.

Figure 3 shows our LSQ scheme’s performances when sampling the different distributions (cf Table 2); the rotation invariant LSQ 6 always achieves the optimal approximation error. This is also the case for the full basis LSQ 4 with \(UUU\) samples (see subfigure 3 (a)). This approximation error is not exactly optimal for the higher degrees, as we are only using a training set of size 8000. However, we can see that \(\mathrm{sym}\) improves this error to the optimal one. Subfigure 3 (b) shows that fitting the full basis is impossible given data sampled from \(\delta UU\) as it is lacking information on Particle 1. Figure 3 (c) shows that this is also the case when the constraint on Particle \(1\) has slightly been relaxed.

On Figure 4 we can see that a quadrature augmented LSQ 7 achieves exact symmetry learning for all sampling distributions. In fact, exact symmetry is achieved immediately after the quadrature rule’s degree overtakes that of the LSQ system. We observe that this numerical experiment also illustrates Remark 2 as there is no meaningful improvements in \(\varepsilon_{\mathrm{sym}}\) before the quadrature degree overcomes the approximation degree. If a quadrature augmentation symmetrised all the lower degree terms of the LSQ approximation, then we would expect an exponential decay in the symmetrisation error due to the decay in the polynomial coefficients - this follows from Remark 1.

In the proofs of Theorem 1 and \(\ref{thm:randaugment}\) we supposed that the design matrix had full rank for the sake of simplicity. However, this is not the case for the full basis LSQ 4 when sampling \(\delta UU\). And when sampling \(\delta^*UU\), while the resulting design matrix has full rank, it is ill-conditioned. For instance, for for an approximation of degree \(9\) and \(800\) data points in the learning set, the smallest singular value is on the order of \(10^{-6}\). This explains the need to regularise. However, data augmentation leads to better conditioned systems and quadrature augmentation will still lead to exact symmetry learning. This is illustrated in Figure 5.

In Figure 6 we have now decreased the training dataset to 100 points. This experiment demonstrates that even for a insufficiently sized dataset, quadrature augmentation will eventually lead to perfect symmetrisation.

Figure 7 shows us the result of random augmentation 7 . In the asymptotic regime we clearly observe the rate \(T^{-1/2}\) as predicted by 2. However, in the pre-asymptotic regime we observe a sharp drop in the symmetrization error which we hypothesize corresponds to the moment all parameters are fully determined. This is especially well illustrated by the concentrated distributions. Indeed, some basis elements are undetermined in the unaugmented dataset since particle 1 is always constrained to be at \(\theta=0\). This can also be observed in the relatively early drop for \(\delta^*UU\) that can be observed in Figure 7 (c) as in that case Particle 1 is not fully constrained and thus these basis elements have relatively more information. Additionally, around that drop the standard deviation is quite high, which can be explained by the fact that the LSQ system is on the cusp of having enough data to fully determine these few basis elements. However, there is one a priori counterintuitive observation. For a fixed augmentation \(\varepsilon_{\mathrm{sym}}\) decreases when the degree increases. One would expect that a lower degree approximation is easier to symmetrise since there are less parameters. The upper bound in 13 may explain this phenomenon. Due to the high regularity of our target function, for a well-fitted invariant LSQ 6 , the term \(\norm{Y-A_{\mathrm{I}}\overline{\beta_{\mathrm{I}}}}\) decays exponentially in \(K\). And so even if \(\norm{\overline{D}_T^*A_{\mathrm{N}}^*}\) were to grow for increasing \(K\), overall, \(\norm{\overline{D}_T^*A_{\mathrm{N}}^*}\norm{Y-A_{\mathrm{I}}\overline{\beta_{\mathrm{I}}}}\) will be decreasing in \(K\) for a sufficiently smooth target. We illustrate this behavior in the appendix (section 6).

Figure 8 compares the symmetrisation error \(\varepsilon_{\mathrm{sym}}\) of both random and quadrature augmentation 7 . While our proofs do not guarantee that quadrature augmentation will be significantly better than random augmentation before exact symmetry learning (cf. Remark 2), we observe that empirically quadrature augmentation is an order of magnitude better than random augmentations. Remark that the square root behaviour is not visible as the system has not yet reached the asymptotic regime.

4.2 Results in Three Dimensions↩︎

The data distributions used in the three-dimensional experiments are summarized in Table 3 and illustrated in Figure 9.

As we work with 3 particles, we may use 5 to build our approximation schemes as it enumerates an orthonormal basis for the space \(B^{3,2}_{K}\). Again, we choose our target polynomial to have an exponential decay of rate \(\alpha=2\) in its coefficients in order to mimic analyticity. That is, we may write our target function \(f\) in the following form: \[f = \sum_{\mathbf{l}}c_{\mathbf{\mathbf{l}}}e^{-\alpha \norm{\mathbf{l}}_1}\tilde{\varphi}_{\mathbf{l}},\] where the \(c_\mathbf{l}\) are chosen uniformly on \([-1,1]\). Throughout our tests we observe qualitatively very similar results as in our two-dimensional tests.

In Figure 10 we recover the spectral convergence rate of the optimal polynomial approximation of an analytic function. Again, we had to use some regularisation for \(\delta^*UU\) and \(\delta^*H_1^*U\). The optimal cutoffs of singular values were found to be, respectively, \(10^{-3.4}\) and \(10^{-5}\).

Next, to test quadrature augmentation, we need to obtain a quadrature rule for the group average over \(\mathrm{SO}(3)\). In general, obtaining quadrature rules on manifolds, for instance Lie Groups, is a non-trivial open problem. The specific case of \(\mathrm{SO}(3)\) was treated in [35] from which we obtained quadrature rules of degree 1 through 11.

In Figure 11 we observe (consistent with the 2D tests) that a quadrature augmentation 7 leads to perfect symmetry learning once the order of the quadrature exceeds the order of the model. Particularly noteworthy is that while \(\delta H_1U\) did not require prior regularisation for an unaugmented LSQ, when the data is not sufficiently augmented the LSQ becomes highly unstable similarly to 2D. We conjecture that the constrained data leaves some basis elements undetermined until the rotations bring enough variety into the data set.

For comparison, Figure 12 shows the evolution of the symmetrisation error under a random augmentation 7 . We observe the same behaviour as for \(d=1\). The expected square root decay is observed as well as the drops in \(\varepsilon_{\mathrm{sym}}\) once enough rotations have been applied.

Figure 13 compares the symmetrisation error \(\varepsilon_{\mathrm{sym}}\) of both random and quadrature augmentation 7 . Just as in Figure 8 we observe that empirically quadrature augmentation is an order of magnitude better than random augmentations. Additionally the square root behaviour can be seen in Figure 13 (a).

5 Conclusion↩︎

We considered the approximation on functions invariant under a group action, in a data-driven setting. We analyzed two methods for augmenting a dataset and least squares system in order to improve the symmetry properties of the resulting approximant: a deterministic quadrature based approach and a probabilistic approach. For both cases we provided rigorous convergence rates. The key take-away is that, when the original data distribution is highly concentrated relative to the symmetry group, then the quadrature approach can provide exact symmetrization but only at the cost of a significant enlargement of the least squares system. In the pre-asymptotic regime it provides no discernable improvement over the probabilistic approach, which unsurprisingly converges at the usual \(T^{-1/2}\) rate with \(T\) begin the number of augmentations (though the explanation for that rate was not immediately obvious). Neither of our approaches make data augmentation both accurate and computationally cheap at the same time.

Following our arguments presented in Section 2, when conservation properties are important in a downstream simulation task, then preserving symmetry as tightly as possible can be crucial and in those situations exact symmetrization, whether by enforcing it in the architecture, or through the proposed quadrature approach, is advisable.

We conclude with the caveat, that the present work is restricted to polynomial models and that the extension of our observations to nonlinear and deep learning methods are not obvious: further theory and empirical tests are required but empirical evidence suggests that quadrature augmentation has potential in that setting as well as it performs an order of magnitude better than random augmentations in our numerical tests.

Acknowledgements↩︎

We are grateful to stimulating discussions with Timon Gutleb and Isaac Holt at early stages of this research project. The authors were supported by the Natural Sciences and Engineering Research Council of Canada through NSERC Discovery grants and an NSERC-NSF Collaboration on quantum science and artificial intelligence grant. The research was also supported through computational resources and services provided by the Digital Research Alliance of Canada (alliancecan.ca) and by Advanced Research Computing at the University of British Columbia. Christoph Ortner is a partner in Symmetric Group LLP which licenses force fields commercially.

6 sec:Appendix↩︎

In this appendix we discuss the observed oddity in Figure 7. We ran simulations, sampling form the uniform distribution, for target functions of different regularity. We simulate targets with coefficients that asymptotically behave like \(K^{-1/2}, K^{-1}, K^{-2}, K^{-3}\), which corresponds to \(\frac{1}{2}\)-Hölder, \(C^1\) and \(C^2\) and \(C^3\) functions.

From Figure 14 we observe the effects of the upper bound 13 . Indeed, we can see that for less regular target functions, our aforementioned intuition was correct, it is easier to symmetrise a low degree approximant. But then, as the regularity increases, the higher degree approximants catch up and eventually overtake the lower degree ones. Showing that regularity makes symmetrisation easier.

References↩︎