Generalizing Convolutional Neural Networks for Equivariance to Lie Groups on Arbitrary Continuous Data

Introduction

Symmetry pervades the natural world. The same law of gravitation governs a game of catch, the orbits of our planets, and the formation of galaxies. It is precisely because of the order of the universe that we can hope to understand it. Once we started to understand the symmetries inherent in physical laws, we could predict behavior in galaxies billions of light-years away by studying our own local region of time and space. For statistical models to achieve their full potential, it is essential to incorporate our knowledge of naturally occurring symmetries into the design of algorithms and architectures. An example of this principle is the translation equivariance of convolutional layers in neural networks (LeCun et al., 1995): when an input (e.g. an image) is translated, the output of a convolutional layer is translated in the same way.

Group theory provides a mechanism to reason about symmetry and equivariance. Convolutional layers are equivariant to translations, and are a special case of group convolution. A group convolution is a general linear transformation equivariant to a given group, used in group equivariant convolutional networks (Cohen and Welling, 2016a).

In this paper, we develop a general framework for equivariant models on arbitrary continuous (spatial) data represented as coordinates and values $\{(x_{i},f_{i})\}_{i=1}^{N}$. Spatial data is a broad category, including ball-and-stick representations of molecules, the coordinates of a dynamical system, and images (shown in Figure 1). When the inputs or group elements lie on a grid (e.g., image data) one can simply enumerate the values of the convolutional kernel at each group element. But in order to extend to continuous data, we define the convolutional kernel as a continuous function on the group parameterized by a neural network.

We consider the large class of continuous groups known as Lie groups. In most cases, Lie groups can be parameterized in terms of a vector space of infinitesimal generators (the Lie algebra) via the logarithm and exponential maps. Many useful transformations are Lie groups, including translations, rotations, and scalings. We propose LieConv, a convolutional layer that can be made equivariant to a given Lie group by defining $\exp$ and $\log$ maps. We demonstrate the expressivity and generality of LieConv with experiments on images, molecular data, and dynamical systems. We emphasize that we use the same network architecture for all transformation groups and data types. LieConv achieves state-of-the-art performance in these domains, even compared to domain-specific architectures. In short, the main contributions of this work are as follows:

We propose LieConv, a new convolutional layer equivariant to transformations from Lie groups. Models composed with LieConv layers can be applied to non-homogeneous spaces and arbitrary spatial data.

We evaluate LieConv on the image classification benchmark dataset rotMNIST (Larochelle et al., 2007), and the regression benchmark dataset QM9 (Blum and Reymond, 2009; Rupp et al., 2012). LieConv outperforms state-of-the-art methods on some tasks in QM9, and in all cases achieves competitive results.

We apply LieConv to modeling the Hamiltonian of physical systems, where equivariance corresponds to the preservation of physical quantities (energy, angular momentum, etc.). LieConv outperforms state-of-the-art methods for the modeling of dynamical systems.

We make code available at https://github.com/mfinzi/LieConv

Related Work

One approach to constructing equivariant CNNs, first introduced in Cohen and Welling (2016a), is to use standard convolutional kernels and transform them or the feature maps for each of the elements in the group. For discrete groups this approach leads to exact equivariance and uses the so-called regular representation of the group (Cohen et al., 2019). This approach is easy to implement, and has also been used when the feature maps are vector fields (Zhou et al., 2017; Marcos et al., 2017), and with other representations (Cohen and Welling, 2016b), but only on image data where locations are discrete and the group cardinality is small. This approach has the disadvantage that the computation grows quickly with the size of the group, and some groups like 3D rotations cannot be easily discretized onto a lattice that is also a subgroup.

Another approach, drawing on harmonic analysis, finds a basis of equivariant functions and parametrizes convolutional kernels in that basis (Worrall et al., 2017; Weiler and Cesa, 2019; Jacobsen et al., 2017). These kernels can be used to construct networks that are exactly equivariant to continuous groups. While the approach has been applied on general data types like spherical images (Esteves et al., 2018; Cohen et al., 2018), voxel data (Weiler et al., 2018), and point clouds (Thomas et al., 2018; Anderson et al., 2019), the requirement of working out the representation theory for the group can be cumbersome and is limited to compact groups. Our approach reduces the amount of work to implement equivariance to a new group, enabling rapid prototyping.

There is also work applying Lie group theory to deep neural networks. Huang et al. (2017) define a network where the intermediate activations of the network are 3D rotations representing skeletal poses and embed elements into the Lie algebra using the $\log$ map. Bekkers (2019) use the $\log$ map to express an equivariant convolution kernel through the use of B-splines, which they evaluate on a grid and apply to image problems. While similar in motivation, their method is not readily applicable to point data and can only be used when the equivariance group acts transitively on the input space. Both of these issues are addressed by our work.

Background

A mapping $h(\cdot)$ is equivariant to a set of transformations $G$ if when we apply any transformation $g$ to the input of $h$, the output is also transformed by $g$. The most common example of equivariance in deep learning is the translation equivariance of convolutional layers: if we translate the input image by an integer number of pixels in $x$ and $y$, the output is also translated by the same amount (ignoring the regions close to the boundary of the image). Formally, if $h:A\rightarrow A$, and $G$ is a set of transformations acting on $A$, we say $h$ is equivariant to $G$ if $\forall a\in A$, $\forall g\in G$,

It is easy to construct invariant functions, where transformations on the input do not affect the output, by simply discarding information. Strict invariance unnecessarily limits the expressive power by discarding relevant information, and instead it is necessary to use equivariant transformations that preserve the information.

2 Groups of Transformations and Lie Groups

3 Group Convolutions

Adopting the convention of left equivariance, one can define a group convolution between two functions on the group, which generalizes the translation equivariance of convolution to other groups:

(Kondor and Trivedi, 2018; Cohen et al., 2019)

4 PointConv Trick

We approximate the integral using a discretization:

In PointConv, Wu et al. (2019) develop a trick where clever reordering of the computation cuts memory and computational requirements by $\sim 2$ orders of magnitude, allowing them to scale to the point cloud classification, segmentation datasets ModelNet40 and ShapeNet, and the image dataset CIFAR-10. We review and generalize the Efficient-PointConv trick in Appendix A.1, which we will use to accelerate our method.

Convolutional Layers on Lie Groups

We begin with a high-level overview of the method. In Section 4.1 we discuss transforming raw inputs $x_{i}$ into group elements $u_{i}$ on which we can perform group convolution. We refer to this process as lifting. Section 4.2 addresses the irregular and varied arrangements of group elements that result from lifting arbitrary continuous input data by parametrizing the convolutional kernel $k$ as a neural network. In Section 4.3, we show how to enforce the locality of the kernel by defining an invariant distance on the group. In Section 4.4, we define a Monte Carlo estimator for the group convolution integral in Eq. (2) and show that this estimator is equivariant in distribution. In Section 4.5, we extend the procedure to cases where the group does not act transitively on the input space (when we cannot map any point to any other point with a transformation from the group). Additionally, in Appendix A.2, we show that our method generalizes coordinate transform equivariance when $G$ is Abelian. At the end of Section 4.5 we provide a concise algorithmic description of the lifting procedure and our new convolution layer.

If $\mathcal{X}$ is a homogeneous space of $G$, then every two elements in $\mathcal{X}$ are connected by an element in $G$, and one can lift elements by simply picking an origin $o$ and defining $\textrm{Lift}(x)=\{u\in G:uo=x\}$: all elements in the group that map the origin to $x$. This procedure enables lifting tuples of coordinates and features $\{(x_{i},f_{i})\}_{i=1}^{N}\rightarrow\{(u_{ik},f_{i})\}_{i=1,k=1}^{N,K}$, with up to $K$ group elements for each input.When $f_{i}=f(x_{i})$, lifting in this way is equivalent to defining $f^{\uparrow}(u)=f(uo)$ as in Kondor and Trivedi (2018). To find all the elements $\{u\in G:uo=x\}$, one simply needs to find one element $u_{x}$ and use the elements in the stabilizer of the origin $H=\{h\in G:ho=o\}$, to generate the rest with $\textrm{Lift}(x)=\{u_{x}h\textrm{ for }h\in H\}$. For continuous groups the stabilizer may be infinite, and in these cases we sample uniformly using the Haar measure $\mu$ which is described in Appendix C.2. We visualize the lifting procedure for different groups in Figure 2.

2 Parameterization of the Kernel

The conventional method for implementing an equivariant convolutional network (Cohen and Welling, 2016a) requires enumerating the values of $k(\cdot)$ over the elements of the group, with separate parameters for each element. This procedure is infeasible for irregularly sampled data and problematic even for a discretization because there is no generalization between different group elements. Instead of having a discrete mapping from each group element to the kernel values, we parametrize the convolutional kernel as a continuous function $k_{\theta}$ using a fully connected neural network with Swish activations, varying smoothly over the elements in the Lie group.

3 Enforcing Locality

Important both to the inductive biases of convolutional neural networks and their computational efficiency is the fact that convolutional filters are local, $k_{\theta}(u_{i}-u_{j})=0$ for $\|u_{i}-u_{j}\|>r$. In order to quantify locality on matrix groups, we introduce the function:

where $\log$ is the matrix logarithm, and $F$ is the Frobenius norm. The function is left invariant, since $d(wu,wv)=\|\log(u^{-1}w^{-1}wv)\|_{F}=d(u,v)$, and is a semi-metric (it does not necessarily satisfy the triangle inequality). In Appendix A.3 we show the conditions under which $d(u,v)$ is additionally the distance along the geodesic connecting $u,v$ , a generalization of the well known formula for the geodesic distance between rotations $\|\log(R_{1}^{T}R_{2})\|_{F}$ (Kuffner, 2004).

4 Discretization of the Integral

5 More Than One Orbit?

In this paper, we consider groups both large and small, and we require the ability to enable or disable equivariances like translations. To achieve this functionality, we need to go beyond the usual setting of homogeneous spaces considered in the literature, where every pair of elements in $\mathcal{X}$ are related by an element in $G$. Instead, we consider the quotient space $Q=\mathcal{X}/G$, consisting of the distinct orbits of $G$ in $\mathcal{X}$ (visualized in Figure 4).When $\mathcal{X}$ is a homogeneous space and the quantity of interest is the quotient with the stabilizer of the origin $H$: $G/H\simeq\mathcal{X}$, which has been examined extensively in the literature. Here we concerned with the separate quotient space $Q=\mathcal{X}/G$, relevant when $\mathcal{X}$ is not a homogeneous space. Each of these orbits $q\in Q$ is a homogeneous space of the group, and when $\mathcal{X}$ is a homogeneous space of $G$ then there is only a single orbit. But in general, there will be many distinct orbits, and lifting should preserve the information on which orbit each point is on.

Since the most general equivariant mappings will use this orbit information, throughout the network the space of elements should not be $G$ but rather $G\times\mathcal{X}/G$, and $x\in\mathcal{X}$ is lifted to the tuples $(u,q)$ for $u\in G$ and $q\in Q$. This mapping may be one-to-one or one-to-many depending on the size of $H$, but will preserve the information in $x$ as $uo_{q}=x$ where $o_{q}$ is the chosen origin for each orbit. General equivariant linear transforms can depend on both the input and output orbit, and equivariance only constrains the dependence on group elements and not the orbits.

When the space of orbits $Q$ is continuous we can write the equivariant integral transform as

When $G$ is the trivial group $\{\textrm{id}\}$, this equation simplifies to the integral transform $h(x)=\int k(x,x^{\prime})f(x^{\prime})dx^{\prime}$ where each element in $\mathcal{X}$ is in its own orbit.

In general, even if $\mathcal{X}$ is a smooth manifold and $G$ is a Lie group it is not guaranteed that $\mathcal{X}/G$ is a manifold (Kono and Ishitoya, 1987). However in practice this is not an issue as we will only have a finite number of orbits present in the data. All we need is an invertible way of embedding the orbit information into a vector space to be fed into $k_{\theta}$. One option is to use an embedding of the orbit origin $o_{q}$, or simply find enough invariants of the group to identify the orbit. To give a few examples:

Discretizing (8) as we did in (7), we get

To recap, Algorithms 1 and 2 give a concise overview of our lifting procedure and our new convolution layer respectively. Please consult Appendix C.1 for additional implementation details.

Applications to Image and Molecular Data

First, we evaluate LieConv on two types of problems: classification on image data and regression on molecular data. With LieConv as the convolution layers, we implement a bottleneck ResNet architecture with a final global pooling layer (Figure 5). For a detailed architecture description, see Appendix C.3. We use the same model architecture for all tasks and achieve performance competitive with task-specific specialized methods.

2 Molecular Data

We first perform an ablation study on the Homo problem of predicting the energy of the highest occupied molecular orbital for the molecules. We apply LieConv with different equivariance groups, combined with SO($3$) data augmentation. The results are reported in Table 3. Of the three groups, our SE($3$) network performs the best. We then apply T($3$)-equivariant LieConv layers to the full range of tasks in the QM9 dataset and report the results in Table 2. We perform competitively with state-of-the-art methods (Gilmer et al., 2017; Schütt et al., 2018; Anderson et al., 2019), with lowest MAE on several of the tasks. See 4 for a demonstration of the equivariance property and efficiency with limited data.

Modeling Dynamical Systems

Accurate transition models for macroscopic physical systems are critical components in control systems (Lenz et al., 2015; Kamthe and Deisenroth, 2017; Chua et al., 2018) and data-efficient reinforcement learning algorithms (Nagabandi et al., 2018; Janner et al., 2019). In this section we show how to enforce conservation of quantities such as linear and angular momentum in the modeling of Hamiltonian systems through LieConv symmetries.

For dynamical systems, the equations of motion can be written in terms of the state $\operatorname{\mathbf{z}}$ and time $t$: $\dot{\operatorname{\mathbf{z}}}=F(\operatorname{\mathbf{z}},t)$. Many physically occurring systems have Hamiltonian structure, meaning that the state can be split into generalized coordinates and momenta $\operatorname{\mathbf{z}}=(\operatorname{\mathbf{q}},\operatorname{\mathbf{p}})$, and the dynamics can be written as

for some choice of scalar Hamiltonian $\mathcal{H}(\operatorname{\mathbf{q}},\operatorname{\mathbf{p}},t)$. $\mathcal{H}$ is often the total energy of the system, and can sometimes be split into kinetic and potential energy terms $\mathcal{H}(\operatorname{\mathbf{q}},\operatorname{\mathbf{p}})=K(\operatorname{\mathbf{p}})+V(\operatorname{\mathbf{q}})$. The dynamics can also be written compactly as $\dot{\operatorname{\mathbf{z}}}=J\nabla_{\operatorname{\mathbf{z}}}\mathcal{H}$ for $J=\begin{bmatrix}0&I\\ -I&0\\ \end{bmatrix}$.

2 Exact Conservation of Momentum

While equivariance is broadly useful as an inductive bias, it has a very special implication for the modeling of Hamiltonian systems. Noether’s Hamiltonian theorem states that each continuous symmetry in the Hamiltonian of a dynamical system has a corresponding conserved quantity (Noether, 1971; Butterfield, 2006). Symmetry with respect to the continuous transformations of translations and rotations lead directly to conservation of the total linear and angular momentum of the system, an extremely valuable property for modeling dynamical systems. In fact, all models that exactly conserve linear and angular momentum must have a corresponding translational and rotational symmetry. See Appendix A.5 for a primer on Hamiltonian symmetries, Noether’s theorem, and the implications in the current setting.

As showed in Section 4, we can construct models that are equivariant to a large variety of continuous Lie Group symmetries, and therefore we can exactly conserve associated quantities like linear and angular momentum. Figure 7 shows that using LieConv layers with a given T($2$) and/or SO($2$) symmetry, the model trajectories conserve linear and/or angular momentum with relative error close to machine epsilon, determined by the integrator tolerance. As there is no corresponding Noether conservation for discrete symmetry groups, discrete approaches to enforcing symmetry (Cohen and Welling, 2016a; Marcos et al., 2017) would not be nearly as effective.

3 Results

For evaluation, we compare a fully-connected (FC) Neural-ODE model (Chen et al., 2018), ODE graph networks (OGN) (Battaglia et al., 2016), Hamiltonian ODE graph networks (HOGN) (Sanchez-Gonzalez et al., 2019), and our own LieConv architecture on predicting the motion of point particles connected by springs as described in (Sanchez-Gonzalez et al., 2019). Figure 6 shows example rollout trajectories, and our quantitative results are presented in Figure 7. In the spring problem $N$ bodies with mass $m_{1},\dots,m_{N}$ interact through pairwise spring forces with constants $k_{1},\dots,k_{N\times N}$. The system preserves energy, linear momentum, and angular momentum. The behavior of the system depends both the values of the system parameters ($\mathbf{s}=(k,m)$) and the initial conditions $\operatorname{\mathbf{z}}_{0}$. The dynamics model must learn not only to predict trajectories across a broad range of initial conditions, but also infer the dependence on varied system parameters, which are additional inputs to the model. We compare models that attempt to learn the dynamics $F_{\theta}(\operatorname{\mathbf{z}},t)=d\mathbf{z}/dt$ directly against models that learn the Hamiltonian as described in section 6.1.

In Figure 7(a) and 7(b) we show that by changing the invariance of our Hamiltonian models, we have direct control over the conservation of linear and angular momentum in the predicted trajectories. Figure 7(c) demonstrates that our method outperforms HOGN, a SOTA architecture for dynamics problems, and achieves significant improvement over the naïve fully-connected (FC) model. We summarize the various models and their symmetries in Table 6. Finally, in Figure 8 we evaluate test MSE of the different models over a range of training dataset sizes, highlighting the additive improvements in generalization from the Hamiltonian, Graph-Network, and equivariance inductive biases successively.

Discussion

We presented a convolutional layer to build networks that can handle a wide variety of data types, and flexibly swap out the equivariance of the model. While the image, molecular, and dynamics experiments demonstrate the generality of our method, there are many exciting application domains (e.g. time-series, geostats, audio, mesh) and directions for future work. We also believe that it will be possible to benefit from the inductive biases of HLieConv models even for systems that do not exactly preserve energy or momentum, such as those found in control systems and reinforcement learning.

The success of convolutional neural networks on images has highlighted the power of encoding symmetries in models for learning from raw sensory data. But the variety and complexity of other modalities of data is a significant challenge in further developing this approach. More general data may not be on a grid, it may possess other kinds of symmetries, or it may contain quantities that cannot be easily combined. We believe that central to solving this problem is a decoupling of convenient computational representations of data as dense arrays from the set of geometrically sensible operations they may have. We hope to move towards models that can ‘see’ molecules, dynamical systems, multi-scale objects, heterogeneous measurements, and higher mathematical objects, in the way that convolutional neural networks perceive images.

MF, SS, PI and AGW are supported by an Amazon Research Award, Amazon Machine Learning Research Award, Facebook Research, NSF I-DISRE 193471, NIH R01 DA048764-01A1, NSF IIS-1910266, NSF 1922658 NRT-HDR: FUTURE Foundations, Translation, and Responsibility for Data Science, and by the United States Department of Defense through the National Defense Science & Engineering Graduate (NDSEG) Fellowship Program. We thank Alex Wang for helpful comments.

References

Appendix A Derivations and Additional Methodology

The matrix notation becomes very cumbersome for manipulating these higher order $n$-dimensional arrays, so we will instead use index notation with Latin indices $i,j,k$ indexing points, Greek indices $\alpha,\beta,\gamma$ indexing feature channels, and $c$ indexing the coordinate dimensions of which there are $d=3$ for PointConv and $d=\textrm{dim}(G)+2\textrm{ dim}(Q)$ for LieConv.$\textrm{dim}(Q)$ is the dimension of the space into which $Q$, the orbit identifiers, are embedded. As the objects are not geometric tensors but simply $n$-dimensional arrays, we will make no distinction between upper and lower indices. After expanding into indices, it should be assumed that all values are scalars, and that any free indices can range over all of the values.

In Wu et al. , it was observed that since $k_{ij}^{\alpha,\beta}$ is the output of an MLP, $k_{ij}^{\alpha,\beta}=\sum_{\gamma}W^{\alpha,\beta}_{\gamma}s_{i,j}^{\gamma}$ for some final weight matrix $W$ and penultimate activations $s_{i,j}^{\gamma}$ ($s_{i,j}^{\gamma}$ is simply the result of the MLP after the last nonlinearity). With this in mind, we can rewrite (12)

In practice, the intermediate number of channels is much less than the product of $c_{in}$ and $c_{out}$: $|\gamma|<|\alpha||\beta|$ and so this reordering of the computation leads to a massive reduction in both memory and compute. Furthermore, $b_{i}^{\gamma,\beta}=\sum_{j}s_{i,j}^{\gamma}f_{j}^{\beta}$ can be implemented with regular matrix multiplication and $h_{i}^{\alpha}=\sum_{\beta,\gamma}W^{\alpha,\beta}_{\gamma}b_{i}^{\gamma,\beta}$ can be also by flattening $(\beta,\gamma)$ into a single axis $\varepsilon$: $h_{i}^{\alpha}=\sum_{\varepsilon}W^{\alpha,\varepsilon}b_{i}^{\varepsilon}$.

The sum over index $j$ can be restricted to a subset $j(i)$ (such as a chosen neighborhood) by computing $f^{\beta}_{(\cdot)}$ at each of the required indices and padding to the size of the maximum subset with zeros, and computing $b_{i}^{\gamma,\beta}=\sum_{j}s_{i,j(i)}^{\gamma}f^{\beta}_{j(i)}$ using dense matrix multiplication. Masking out of the values at indices $i$ and $j$ is also necessary when there are different numbers of points per minibatch but batched together using zero padding. The generalized PointConv trick can thus be applied in batch mode when there may be varied number of points per example and varied number of points per neighborhood.

A.2 Abelian G𝐺G and Coordinate Transforms

This directly generalizes some of the existing coordinate transform methods for achieving equivariance from the literature such as log polar coordinates for rotation and scaling equivariance [Esteves et al., 2017], and using hyperbolic coordinates for squeeze and scaling equivariance.

Now writing out our Monte Carlo estimation of the integral:

Again $\mathcal{X}$ is a homogeneous space of $G$, and we choose a single origin $o=$. With a little algebra, it is clear that $M(r_{i},s_{i})o=p_{i}$ where $r=\sqrt{xy}$ and $s=\sqrt{x/y}$ are the hyperbolic coordinates of $p_{i}$.

Expressed in the basis $\mathcal{B}=[I,A]$ for the Lie algebra above, we see that

which is equivariant to squeezes and scalings.

As demonstrated, equivariance to groups that contain the input space in a single orbit and are abelian can be achieved with a simple coordinate transform; however our approach generalizes to groups that are both ’larger’ and ’smaller’ than the input space, including coordinate transform equivariance as a special case.

A.3 Sufficient Conditions for Geodesic Distance

$d(u,v)=0\Leftrightarrow\log(u^{-1}v)=0\Leftrightarrow u=v$

$d(u,v)=\|\log(v^{-1}u)\|=\|-\log(u^{-1}v)\|=d(v,u)$.

where $-T$ denotes inverse and transpose.

Specifically, if the subgroup $G$ is in the image of the $\exp:\mathfrak{g}\to G$ map and each infinitesmal generator commutes with its transpose: $[A,A^{T}]=0$ for $\forall A\in\mathfrak{g}$, then $d(u,v)=\|\log(v^{-1}u)\|_{F}$ is the geodesic distance between $u,v$.

Geodesic Equation: Geodesics of (16) satisfying $\nabla_{\dot{\gamma}}\dot{\gamma}=0$ can equivalently be derived by minimizing the energy functional

using the calculus of variations. Minimizing curves $\gamma(t)$, connecting elements $u$ and $v$ in $G$ ($\gamma(0)=v,\gamma(1)=u$) satisfy

Noting that $\delta(\gamma^{-1})=-\gamma^{-1}\delta\gamma\gamma^{-1}$ and the linearity of the trace,

Using the cyclic property of the trace and integrating by parts, we have that

where the boundary term \operatorname{Tr}(\dot{\gamma}\gamma^{-T}\gamma^{-1}\delta\gamma)\big{|}_{0}^{1} vanishes since $(\delta\gamma)(0)=(\delta\gamma)(1)=0$.

As $\delta\gamma$ may be chosen to vary arbitrarily along the path, $\gamma$ must satisfy the geodesic equation:

Solutions: When $A=\log(v^{-1}u)$ satisfies $[A,A^{T}]=0$, the curve $\gamma(t)=v\exp(t\log(v^{-1}u))$ is a solution to the geodesic equation (17). Clearly $\gamma$ connects $u$ and $v$, $\gamma(0)=v$ and $\gamma(1)=u$. Plugging in $\dot{\gamma}=\gamma A$ into the left hand side of equation (17), we have

Length of $\gamma$: The length of the curve $\gamma$ connecting $u$ and $v$ is $\|\log(v^{-1}u)\|_{F}$,

A.4 Equivariant Subsampling

Even if all distances and neighborhoods are precomputed, the cost of computing equation (6) for $i=1,...,N$ is still quadratic, $O(nN)=O(N^{2})$, because the number of points in each neighborhood $n$ grows linearly with $N$ as $f$ is more densely evaluated. So that our method can scale to handle a large number of points, we show two ways two equivariantly subsample the group elements, which we can use both for the locations at which we evaluate the convolution and the locations that we use for the Monte Carlo estimator. Since the elements are spaced irregularly, we cannot readily use the coset pooling method described in [Cohen and Welling, 2016a], instead we can perform:

Random Selection: Randomly selecting a subset of $p$ points from the original $n$ preserves the original sampling distribution, so it can be used.

Farthest Point Sampling: Given a set of group elements $S=\{u_{i}\}_{i=1}^{k}\in G$, we can select a subset $S_{p}^{*}$ of size $p$ by maximizes the minimum distance between any two elements in that subset,

farthest point sampling on the group. Acting on a set of elements, $\textrm{Sub}_{p}:S\mapsto S_{p}^{*}$, the farthest point subsampling is equivariant $\textrm{Sub}_{p}(wS)=w\textrm{Sub}_{p}(S)$ for any $w\in G$. Meaning that applying a group element to each of the elements does not change the chosen indices in the subsampled set because the distances are left invariant $d(u_{i},u_{j})=d(wu_{i},wu_{j})$.

Now we can use either of these methods for $\textrm{Sub}_{p}(\cdot)$ to equivariantly subsample the quadrature points in each neighborhood used to estimate the integral to a fixed number $p$,

Doing so has reduced the cost of estimating the convolution from $O(N^{2})$ to $O(pN)$, ignoring the cost of computing $\textrm{Sub}_{p}$ and $\{\textrm{nbhd}(u_{i})\}_{i=1}^{N}$.

A.5 Review and Implications of Noether’s Theorem

In the Hamiltonian setting, Noether’s theorem relates the continuous symmetries of the Hamiltonian of a system with conserved quantities, and has been deeply impactful in the understanding of classical physics. We give a review of Noether’s theorem, loosely following Butterfield .

As introduced earlier, the Hamiltonian is a function acting on the state $H(z)=H(q,p)$, (we will ignore time dependence for now) can be viewed more formally as a function on the cotangent bundle $(q,p)=z\in M=T^{*}C$ where $C$ is the coordinate configuration space, and this is the setting for Hamiltonian dynamics.

In general, on a manifold $\mathcal{M}$, a vector field $X$ can be viewed as an assignment of a directional derivative along $\mathcal{M}$ for each point $z\in\mathcal{M}$. It can be expanded in a basis using coordinate charts $X=\sum_{\alpha}X^{\alpha}\partial_{\alpha}$, where $\partial_{\alpha}=\frac{\partial}{\partial z^{\alpha}}$ and acts on functions $f$ by $X(f)=\sum_{\alpha}X^{\alpha}\partial_{\alpha}f$. In the chart, each of the components $X^{\alpha}$ are functions of $z$.

In Hamiltonian mechanics, for two functions on $M$, there is the Poisson bracket which can be written in terms of the canonical coordinates $q_{i},p_{i}$, Here we take the definition of the Poisson bracket to be negative of the usual definition in order to streamline notation.

The Poisson bracket can be used to associate each function $f$ to a vector field

which specifies, by its action on another function $g$, the directional derivative of $g$ along $X_{f}$: $X_{f}(g)=\{f,g\}$. Vector fields that can be written in this way are known as Hamiltonian vector fields, and the Hamiltonian dynamics of the system is a special example $X_{H}=\{H,\cdot\}$. This vector field in canonical coordinates $z=(p,q)$ is the vector field $X_{H}=F(z)=J\nabla_{z}H$ (i.e. the symplectic gradient, as discussed in Section 6.1). Making this connection clear, a given scalar quantity evolves through time as $\dot{f}=\{H,f\}$. But this bracket can be used to evaluate the rate of change of a scalar quantity along the flows of vector fields other than the dynamics, such as the flows of continuous symmetries.

A scalar function is invariant to the flow of a vector field if and only if the Lie Derivative is zero

by the antisymmetry of the Poisson bracket. So if $\phi^{X}_{\lambda}$ is a symmetry of $H$, then $X=X_{f}$ for some function $f$, and $H(\phi^{X_{f}}_{\lambda}(z))=H(z)$ implies

or in other words $f(z(t+\tau))=f(z(t))$ and $f$ is a conserved quantity of the dynamics.

This implication goes both ways, if $f$ is conserved then $\phi^{X_{f}}_{\lambda}$ is necessarily a symmetry of the Hamiltonian, and if $\phi^{X_{f}}_{\lambda}$ is a symmetry of the Hamiltonian then $f$ is conserved.

So far we have been discussing Hamiltonian symmetries, invariances of the Hamiltonian. But in the study of dynamical systems there is a related concept of dynamical symmetries, symmetries of the equations of motion. This notion is also captured by the Lie Derivative, but between vector fields. A dynamical system $\dot{z}=F(z)$, has a continuous dynamical symmetry $\phi^{X}_{\lambda}$ if the flow along the dynamical system commutes with the symmetry:

Meaning that applying the symmetry transformation to the state and then flowing along the dynamical system is equivalent to flowing first and then applying the symmetry transformation. Equation (20) is satisfied if and only if the Lie Derivative is zero:

where $[\cdot,\cdot]$ is the Lie bracket on vector fields.The Lie bracket on vector fields produces another vector field and is defined by how it acts on functions, for any smooth function $g$: $[X,F](g)=X(F(g))-F(X(g))$

For Hamiltonian systems, every Hamiltonian symmetry is also a dynamical symmetry. In fact, it is not hard to show that the Lie and Poisson brackets are related,

and this directly shows the implication. If $X_{f}$ is a Hamiltonian symmetry, $\{f,H\}=0$, and then

However, the converse is not true, dynamical symmetries of a Hamiltonian system are not necessarily Hamiltonian symmetries and thus might not correspond to conserved quantities. Furthermore even if the system has a dynamical symmetry which is the flow along a Hamiltonian vector field $\phi^{X}_{\lambda}$, $X=X_{f}=\{f,\cdot\}$, but the dynamics $F$ are not Hamiltonian, then the dynamics will not conserve $f$ in general. Both the symmetry and the dynamics must be Hamiltonian for the conservation laws.

This fact is demonstrated by Figure 9, where the dynamics of the (non-Hamiltonian) equivariant LieConv-T($2$) model has a T($2$) dynamical symmetry with the generators $\partial_{x},\partial_{y}$ which are Hamiltonian vector fields for $f=p_{x},f=p_{y}$, and yet linear momentum is not conserved by the model.

Consider a system of $N$ interacting particles described in Euclidean coordinates with position and momentum $q_{im},p_{im}$, such as the multi-body spring problem. Here the first index $i=1,2,3$ indexes the spatial coordinates and the second $m=1,2,...,N$ indexes the particles. We will use the bolded notation $\mathbf{q}_{m},\mathbf{p}_{m}$ to suppress the spatial indices, but still indexing the particles $m$ as in Section 6.1.

The total linear momentum along a given direction $\mathbf{n}$ is $\mathbf{n}\cdot\mathbf{P}=\sum_{i,m}n_{i}p_{im}=\mathbf{n}\cdot(\sum_{m}\mathbf{p}_{m})$. Expanding the Poisson bracket, the Hamiltonian vector field

which has the flow $\phi^{X_{\mathbf{n}\mathbf{P}}}_{\lambda}(\mathbf{q}_{m},\mathbf{p}_{m})=(\mathbf{q}_{m}+\lambda\mathbf{n},\mathbf{p}_{m})$, a translation of all particles by $\lambda\mathbf{n}$. So our model of the Hamiltonian conserves linear momentum if and only if it is invariant to a global translation of all particles, (e.g. T($2$) invariance for a 2D spring system).

The total angular momentum along a given axis $\mathbf{n}$ is

, where $\epsilon_{ijk}$ is the Levi-Civita symbol and we have defined the antisymmetric matrix $A$ by $A_{kj}=\sum_{i}\epsilon_{ijk}n_{i}$.

where the second line follows from the antisymmetry of $A$. We can find the flow of $X_{\mathbf{n}\mathbf{L}}$ from the differential equations $\dot{\mathbf{q}}_{m}=A\mathbf{q},\dot{\mathbf{p}}_{m}=A\mathbf{q}$ which have the solution

where $R_{\theta}$ is a rotation about the axis $\mathbf{n}$ by the angle $\theta$, which follows from the Rodriguez rotation formula. Therefore, the flow of the Hamiltonian vector field of angular momentum along a given axis is a global rotation of the position and momentum of each particle about that axis. Again, the dynamics of a neural network modeling a Hamiltonian conserve total angular momentum if and only if the network is invariant to simultaneous rotation of all particle positions and momenta.

Appendix B Additional Experiments

While (7) shows that the convolution estimator is equivariant, we have conducted the ablation study below examining the equivariance of the network empirically. We trained LieConv (Trivial, T($3$), SO($3$), SE($3$)) models on a limited subset of 20k training examples (out of 100k) of the HOMO task on QM9 without any data augmentation. We then evaluate these models on a series of modified test sets where each example has been randomly transformed by an element of the given group (the test translations in T($3$) and SE($3$) are sampled from a normal with stddev 0.5). In table 4 the rows are the models configured with a given group equivariance and the columns N/G denote no augmentation at training time and transformations from G applied to the test set (test translations in T($3$) and SE($3$) are sampled from a normal with stddev 0.5).

Notably, the performance of the LieConv-G models do not degrade when random G transformations are applied to the test set. Also, in this low data regime, the added equivariances are especially important.

B.2 RotMNIST Comparison

While the RotMNIST dataset consists of 12k rotated MNIST digits, it is standard to separate out 10k to be used for training and 2k for validation. However, in Ti-Pooling and E(2)-Steerable CNNs, it appears that after hyperparameters were tuned the validation set is folded back into the training set to be used as additional training data, a common approach used on other datasets. Although in table 1 we only use 10k training points, in the table below we report the performance with and without augmentation trained on the full 12k examples.

Appendix C Implementation Details

While the high-level summary of the lifting procedure (Algorithm 1) and the LieConv layer (Algorithm 2) provides a useful conceptual understanding of our method, there are some additional details that are important for a practical implementation.

According to Algorithm 2, $a_{ij}$ is computed in every LieConv layer, which is both highly redundant and costly. In practice, we precompute $a_{ij}$ once after lifting and feed it through the network with layers operating on the state \big{(}\{a_{ij}\}_{i,j}^{N,N},\{f_{i}\}_{i=1}^{N}\big{)} instead of $\{(u_{i},q_{i},f_{i})\}_{i=1}^{N}$. Doing so requires fixing the group elements that will be used at each layer for a given forwards pass.

We use the analytic forms for the exponential and logarithm maps of the various groups as described in Eade .

C.2 Sampling from the Haar Measure for Various groups

In this paper, the groups we use in which the lifting map is multi-valued are SE($2$), SO($3$), and SE($3$). The process is especially straightforward for SE($2$) and SE($3$) as these groups can be expressed as a semi-direct product of two groups $G=H\ltimes N$,

C.3 Model Architecture

We employ a ResNet-style architecture [He et al., 2016], using bottleneck blocks [Zagoruyko and Komodakis, 2016], and replacing ReLUs with Swish activations [Ramachandran et al., 2017]. The convolutional kernel $g_{\theta}$ internal to each LieConv layer is parametrized by a 3-layer MLP with 32 hidden units, batch norm, and Swish nonlinearities. Not only do the Swish activations improve performance slightly, but unlike ReLUs they are twice differentiable which is a requirement for backpropagating through the Hamiltonian dynamics. The stack of elementwise linear and bottleneck blocks is followed by a global pooling layer that computes the average over all elements, but not over channels. Like for regular image bottleneck blocks, the channels for the convolutional layer in the middle are smaller by a factor of 4 for increased parameter and computational efficiency.

Downsampling: As is traditional for image data, we increase the number of channels and the receptive field at every downsampling step. The downsampling is performed with the farthest point downsampling method described in Appendix A.4. For a downsampling by a factor of $s<1$, the radius of the neighborhood is scaled up by $s^{-1/2}$ and the channels are scaled up by $s^{-1/2}$. When an image is downsampled with $s=(1/2)^{2}$ that is typical in a CNN, this results in 2x more channels and a radius or dilation of 2x. In the bottleneck block, the downsampling operation is fused with the LieConv layer, so that the convolution is only evaluated at the downsampled query locations. We perform downsampling only on the image datasets, which have more points.

BatchNorm: In order to handle the varied number of group elements per example and within each neighborhood, we use a modified batchnorm that computes statistics only over elements from a given mask. The batch norm is computed per channel, with statistics averaged over the batch size and each of the valid locations.

C.4 Details for Hamiltonian Models

As the position vectors are mean centered in the model forward pass $\operatorname{\mathbf{q}}_{i}^{\prime}=\operatorname{\mathbf{q}}_{i}-\;\bar{\operatorname{\mathbf{q}}}$, HOGN and HLieConv-SO2* have additional T($2$) invariance, yielding SE($2$) invariance for HLieConv-SO2*. We also experimented with a HLieConv-SE2 equivariant model, but found that the exponential map for SE2 (involving taylor expands and masking) was not numerically stable enough for for second derivatives, required for optimizing through the Hamiltonian dynamics. So instead we benchmark the HLieConv-SO2 (without centering) and the HLieConv-SO2* (with centering) models separately. Layer equivariance is preferable for not prematurely discarding useful information and for better modeling performance, but invariance alone is sufficient for the conservation laws. Additionally, since we know a priori that the spring problem has Euclidean coordinates, we need not model the kinetic energy $K(\operatorname{\mathbf{p}},m)=\sum_{j=1}^{n}\|\operatorname{\mathbf{p}}_{j}\|^{2}/m_{j}$ and instead focus on modeling the potential $V(\operatorname{\mathbf{q}},k)$. We observe that this additional inductive bias of Euclidean coordinates improves model performance. Table 6 shows the invariance and equivariance properties of the relevant models and baselines. For Noether conservation, we need both to model the Hamiltonian and have the symmetry property.

Dataset Generation: To generate the spring dynamics datasets we generated $D$ systems each with $N=6$ particles connected by springs. The system parameters, mass and spring constant, are set by sampling $\{m_{1}^{(i)},\dots m_{6}^{(i)},k_{1}^{(i)},\dots,k_{6}^{(i)}\}_{i=1}^{N}$, $m_{j}^{(i)}\sim\mathcal{U}(0.1,3.1)$, $k_{j}^{(i)}\sim\mathcal{U}(0,5)$. Following Sanchez-Gonzalez et al. , we set the spring constants as $k_{ij}=k_{i}k_{j}$. For each system $i$, the position and momentum of body $j$ were distributed as $\operatorname{\mathbf{q}}_{j}^{(i)}\sim\mathcal{N}(0,0.16I)$, $\operatorname{\mathbf{p}}_{j}^{(i)}\sim\mathcal{N}(0,0.36I)$. Using the analytic form of the Hamiltonian for the spring problem, $\mathcal{H}(\operatorname{\mathbf{q}},\operatorname{\mathbf{p}})=K(\operatorname{\mathbf{p}},m)+V(\operatorname{\mathbf{q}},k)$, we use the RK4 numerical integration scheme to generate $5$ second ground truth trajectories broken up into 500 evaluation timesteps. We use a fixed step size scheme for RK4 chosen automatically (as implemented in Chen et al. ) with a relative tolerance of 1e-8 in double precision arithmetic. We then randomly selected a single segment for each trajectory, consisting of an initial state $\operatorname{\mathbf{z}}_{t}$ and $\tau=4$ transition states: $(\operatorname{\mathbf{z}}_{t+1}^{(i)},\dots,\operatorname{\mathbf{z}}_{t+\tau}^{(i)})$.

Training: All models were trained in single precision arithmetic (double precision did not make any appreciable difference) with an integrator tolerance of 1e-4. We use a cosine decay for the learning rate schedule and perform early stopping over the validation MSE. We trained with a minibatch size of $200$ and for $100$ epochs each using the Adam optimizer [Kingma and Ba, 2014] without batch normalization. With 3k training examples, the HLieConv model takes about 20 minutes to train on one 1080Ti.

For the examination of performance over the range of dataset sizes in 8, we cap the validation set to the size of the training set to make the setting more realistic, and we also scale the number of training epochs up as the size of the dataset shrinks (epochs $=100(\sqrt{10^{3}/D})$) which we found to be sufficient to fit the training set. For $D\leq 200$ we use the full dataset in each minibatch.

Hyperparameter tuning: Model hyperparameters were tuned by grid search over channel width, number of layers, and learning rate. The models were tuned with training, validation, and test datasets consisting of 3000, 2000, and 2000 trajectory segments respectively.

C.5 Details for Image and Molecular Experiments

RotMNIST Hyperparameters: For RotMNIST we train each model for 500 epochs using the Adam optimizer with learning rate 3e-3 and batch size 25. The first linear layer maps the 1-channel grayscale input to $k=128$ channels, and the number of channels in the bottleneck blocks follow the scaling law from Appendix C.3 as the group elements are downsampled. We use 6 bottleneck blocks, and the total downsampling factor $S=1/10$ is split geometrically between the blocks as $s=(1/10)^{1/6}$ per block. The initial radius $r$ of the local neighborhoods in the first layer is set so as to include 1/15 of the total number of elements in each neighborhood and is scaled accordingly. The subsampled neighborhood used to compute the Monte Carlo convolution estimator uses $p=25$ elements. The models take less than 12 hours to train on a 1080Ti.

QM9 Hyperparameters: For the QM9 molecular data, we use the featurization from Anderson et al. , where the input features $f_{i}$ are determined by the atom type (C,H,N,O,F) and the atomic charge. The coordinates $x_{i}$ are simply the raw atomic coordinates measured in angstroms. A separate model is trained for each prediction task, all using the same hyperparameters and early stopping on the validation MAE. We use the same train, validation, test split as Anderson et al. , with 100k molecules for train, 10% for test and the remaining for validation. Like with the other experiments, we use a cosine learning rate decay schedule. Each model is trained using the Adam optimizer for $1000$ epochs with a learning rate of 3e-3 and batch size of 100. We use SO($3$) data augmentation, 6 bottleneck blocks, each with $k=1536$ channels. The radius of the local neighborhood is set to $r=\infty$ to include all elements. The model takes about 48 hours to train on a single 1080Ti.

C.6 Local Neighborhood Visualizations