Representing and Learning
Functions Invariant Under Crystallographic Groups
Ryan P Adams and Peter Orbanz

Crystallographic groups describe the symmetries of
crystals and other repetitive structures encountered in nature and the sciences.
These groups include the wallpaper and space groups.
We derive linear and nonlinear representations of functions that are (1) smooth and
(2) invariant under such a group.
The linear representation generalizes the Fourier basis to crystallographically invariant
basis functions. We show that such a basis exists for each crystallographic group,
that it is orthonormal in the relevant $\L_2$ space, and recover the standard Fourier
basis as a special case for pure shift groups. The nonlinear representation
embeds the orbit space of the group into a finite-dimensional Euclidean space.
We show that such an embedding exists for every crystallographic group, and that it factors functions
through a generalization of a manifold called an orbifold.
We describe algorithms that, given a standardized description of the group,
compute the Fourier basis and an embedding map. As examples,
we construct crystallographically invariant neural networks, kernel machines,
and Gaussian processes.

§ INTRODUCTION

Among the many forms of symmetry observed in nature, those that arise from repetitive spatial patterns
are particularly important. These are described by sets of transformations
of Euclidean space called crystallographic groups [Vinberg and Shvartsman, 1993, Thurston, 1997].
For example, consider a problem in materials science, where atoms are arranged in a crystal lattice.
The symmetries of the lattice are then characterized by a crystallographic group $\group$.
Symmetry means that, if we apply one of the transformations in $\group$
to move the lattice—say to rotate or shift it—the transformed lattice is indistinguishable from the
untransformed one.
In such a lattice, the Coulomb potential acting on any single electron due to a collection of fixed nuclei
does not change under any of the transformations in $\group$ [Bell, 1954, Lax, 2001, Johnston, 1960].
If we think of the potential field as a function
on $\mathbb{R}^3$, this is an example of a $\group$-invariant function, i.e., a function whose values do
not change if its arguments are transformed by elements of the group.
When solving the resulting Schrödinger equation for single particle states, members of the group commute with the Hamiltonian, and quantum observables are again $\group$-invariant [Bell, 1954, Killingbeck, 1970, Johnston, 1960, Slater, 1965, Lax, 2001, Dresselhaus et al., 2007].
A different example are ornamental tilings on the walls of the Alhambra, which, when regarded
as functions on $\mathbb{R}^2$, are invariant under two-dimensional crystallographic groups
[Speiser, 1956].
The purpose of this work is to construct smooth invariant functions for any given crystallographic
group in any dimension.

For finite groups, invariant functions can be constructed easily by summing over all group elements;
for compact infinite groups, the sum can be replaced by an integral. This and related ideas have
received considerable attention in machine learning
Such summations are not possible for crystallographic groups, which are neither finite nor compact,
but their specific algebraic and geometric properties
allow us to approach the problem in a different manner.
We postpone a detailed literature review to <ref>,
and use the remainder of this section to give a non-technical sketch of our results.

§.§ A non-technical overview

This section sketches our results in a purely heuristic way; proper
definitions follow in <ref>.

Crystallographic symmetry.
Crystallographic groups are groups that tile a Euclidean space $\mathbb{R}^n$ with a convex shape.
Suppose we place a convex polytope $\Pi$ in the space $\mathbb{R}^n$, say a square or a rectangle in the plane.
Now make a copy of $\Pi$, and use a transformation ${\phi:\mathbb{R}^2\rightarrow\mathbb{R}^2}$ to
move this copy to another location.
We require that $\phi$ is an isometry, which means it may shift, rotate or flip $\Pi$, but does not
change its shape or size.
Here are some examples, where the original shape $\Pi$ is marked in red:

at (0,0) ;
at (4.5,0) ;
at (9,0) ;
at (0,-1.65) ;
at (4.5,-1.65) ;
at (9,-1.65) ;

The descriptors , , and  follow the naming standard for groups
developed by crystallographers [Hahn, 2011], and the symbol “F” is inscribed
to clarify which transformations are used.
The transformations in these examples are horizontal and vertical shifts (in ), rotations around the
corners of the rectangle in (), and reflections about its edges ().
Suppose we repeat one of these processes indefinitely so that the copies of $\Pi$ cover the entire plane
and overlap only on their boundaries. That requires a countably infinite number of
transformations, one per copy. Collect these into a set $\group$.
If this set forms a group, this group is called crystallographic.
Such groups describe all possible symmetries of crystals, and have been thoroughly studied in
crystallography. For each dimension $n$, there is—up to a natural notion of isomorphy that we explain in <ref>—only
a finite number of crystallographic groups: Two on $\mathbb{R}$, 17 on $\mathbb{R}^2$,
230 on $\mathbb{R}^3$, and so forth.
Those on $\mathbb{R}^2$ are also known as wallpaper groups, and those on $\mathbb{R}^3$ as space groups.

The objects of interest.
A function $f$ is invariant under  $\group$ if it satisfies

\begin{equation*}
f(\phi x)\;=\;f(x)\qquad\text{ for all }\phi\in\group\text{ and all }x\in\mathbb{R}^n\;.
\end{equation*}

A simple way to construct such a function is to start with a tiling as above, define a function
on $\Pi$, and then replicate it on every copy of $\Pi$. Here are two examples on $\mathbb{R}^2$,
corresponding to (ii) and (iii) above, and an example on $\mathbb{R}^3$:

at (0,0) ;
at (5,0) ;
at (10,0) ;
at (0,-2) ;
at (5,-2) ;
at (11.5,-2) ;

However, as the examples illustrate,
functions obtained this way are typically not continuous.
Our goal is to construct smooth invariant functions, such as these:

at (0,0) ;
at (5,0) ;
at (10,0) ;
at (0,-2) ;
at (5,-2) ;
at (11.5,-2) ;

We identify two representations of such functions, one linear and one nonlinear.
Working with either representation algorithmically requires a data structure representing the invariance
constraint. We construct such a structure, which we call an orbit graph,
in <ref>. This graph is constructed from a description of the group
(which can be encoded as a finite set of matrices) and of $\Pi$ (a finite set of vectors).

Linear representations: Invariant Fourier transforms.
We are primarily interested in two and three dimensions, but a one-dimensional example
is a good place to start: In one dimension, a convex polytope is always an interval, say ${\Pi=[0,1]}$. If we choose
$\group$ as the group $\mathbb{Z}$ of all shifts of integer length, it tiles the line $\mathbb{R}$ with $\Pi$.
In this case, an invariant function is simply a periodic function with period $1$. Smooth periodic
functions can be represented as a Fourier series,

\begin{equation*}
\end{equation*}

for sequences of scalar coefficients $a_i$ and $b_i$.
Note each sine and cosine on the right is $\group$-invariant and infinitely often differentiable.
Now suppose we abstract from the specific form of these sines and cosines, and only
regard them as $\group$-invariant functions that are very smooth.
The series representation above then has the general form

\begin{equation*}
\end{equation*}

where the $\ef_i$ are smooth, $\group$-invariant functions that depend only on $\group$ and $\Pi$,
and the $c_i$ are scalar coefficients that depend on $f$. (In the Fourier series, $\ef_i$ is
a cosine for odd and a sine for even indices.) In <ref>, we obtain generalizations of
this representation to crystallographically invariant functions. To do so,
we observe that the Fourier basis can be derived as the set of eigenfunctions of the Laplace operator:
The sine and cosine functions above are precisely those functions ${\ef:\mathbb{R}\rightarrow\mathbb{R}}$ that solve

\begin{align*}
-\Delta \ef&\;=\;\lambda\ef \\
\text{ subject to }\quad\qquad\ef&\text{ is periodic with period }1
\end{align*}

for some ${\lambda\geq 0}$. (The negative sign is chosen to make the eigenvalues $\lambda$ non-negative.)
The periodicity constraint is equivalent to saying that $\ef$ is invariant
under the shift group ${\group=\mathbb{Z}}$. The corresponding problem for a general crystallographic group $\group$ on $\mathbb{R}^n$ is hence

\begin{equation}
\label{eigenproblem:intro}
\begin{split}
-\Delta \ef\;&=\;\lambda\ef \\
\text{ subject to }\quad\qquad\ef\;&=\;\ef\circ\phi \text{ for all }\phi\in\group\;.
\end{split}
\end{equation}

<ref> shows that this problem has solutions for any dimension $n$, convex polytope ${\Pi\subset\mathbb{R}^n}$, and crystallographic
group $\group$ that tiles $\mathbb{R}^n$ with $\Pi$. As in the Fourier case, the solution functions ${\ef_1,\ef_2,\ldots}$ are very smooth.

If we choose ${\Pi\subset\mathbb{R}^2}$ as the square ${[0,1]^2}$ and $\group$ as the group $\mathbb{Z}^2$ of discrete horizontal and vertical shifts—that is, the two-dimensional analogue of the example above—we recover the
two-dimensional Fourier transform. The function $\ef_0$ is constant; the functions ${\ef_1,\ldots,\ef_5}$ are shown in
If the group also contains other transformations, the basis looks less familiar.
These are the basis functions ${\ef_1,\ldots,\ef_5}$ for a group () containing shifts and rotations of order three:

The same idea applies in any finite dimension $n$. For ${n=3}$, the $\ef_i$ can be visualized
as contour plots. For instance, the first five non-constant basis elements for a specific three-dimensional
group, designated  by crystallographers, look like this:

Our results show that any continuous invariant function can be represented
by a series expansion in functions $\ef_i$. As for the Fourier transform, the functions form a orthonormal
basis of the relevant $\L_2$ space. The functions $\ef_i$ can hence be seen as a generalization of the
Fourier transform from pure shift groups to crystallographic groups. All of this is made precise in <ref>.

Nonlinear representations: Factoring through an orbifold.
The second representation, in <ref>, generalizes an idea of David MacKay [MacKay, 1998], who constructs periodic
functions on the line as follows: Start with a continuous function ${h:\mathbb{R}^2\rightarrow\mathbb{R}}$.
Choose a circle of circumference $1$ in $\mathbb{R}^2$, and restrict $h$ to the circle. The restriction is still
continuous. Now “cut and unfold the circle with $h$ on it” to obtain a function on the unit interval.
Since this function takes the same value at both interval boundaries, replicating it by shifts of integer length
defines a function on $\mathbb{R}$ that is periodic and continuous:

at (0,0) ;
at (3.5,0) ;
at (9,0) ;
at (0,-2) function $h$ on $\mathbb{R}^2$;
at (3.5,-2) restrict $h$ to circle;
at (9,-2) unfold circle and replicate;

More formally, MacKay's approach constructs a function
${\rho:\mathbb{R}\rightarrow\text{circle}\subset\mathbb{R}^2}$ such that

\begin{equation*}
f\text{ is continuous and periodic on $\mathbb{R}$ }
\quad\Leftrightarrow\quad
f\;=\;h\circ\rho\quad\text{ for some continuous }h:\mathbb{R}^2\rightarrow\mathbb{R}\;.
\end{equation*}

We show how to generalize this  construction to any finite dimension $n$, any crystallographic
group $\group$ on $\mathbb{R}^n$, and any convex polytope with which $\group$ tiles the space:
For each $\group$ and $\Pi$, there is a continuous, surjective map

\begin{equation}
\label{intro:embedding:map}
\rho:\mathbb{R}^n\rightarrow\Omega\qquad\text{ for some finite }N\geq n\text{ and a compact set }\Omega\subset\mathbb{R}^N
\end{equation}

such that

\begin{equation*}
f\text{ is continuous and invariant }
\quad\Leftrightarrow\quad
f\;=\;h\circ\rho\quad\text{ for some continuous }h:\mathbb{R}^N\rightarrow\mathbb{R}\;.
\end{equation*}

This is <ref>.
<ref> shows how to compute a representation of $\rho$ using multidimensional scaling.

The set $\Omega$ can be thought of as an $n$-dimensional surface in a higher-dimensional space $\mathbb{R}^N$.
If $\group$ contains only shifts, this surface is completely smooth, and hence a manifold. That is the
case in MacKay's construction, where $\Omega$ is the circle, and
the group  on $\mathbb{R}^2$, for which $\Omega$ is the torus shown on the left:

at (0,0) ;
at (6,0) ;

For most crystallographic groups, $\Omega$ is not a manifold, but
rather a more general object called an orbifold. The precise definition (see <ref>) is somewhat technical, but
loosely speaking, an orbifold is a surface that resembles a manifold almost everywhere, except at a small number of
points at which it is not smooth.
That is illustrated by the orbifold on the right, which represents a group containing rotations, and has several “sharp corners”.

Applications I: Neural networks.
We can now define $\group$-invariant models by factoring through $\rho$.
To define an invariant
neural network, for example, start with a continuous neural network ${h_{\theta}:\mathbb{R}^N\rightarrow Y}$ with weight vector $\theta$ and some
output space $Y$. Then ${\rho\circ h_{\theta}}$ is a continuous and invariant neural network ${\mathbb{R}^n\rightarrow Y}$.
Here are examples for three groups (, , and ) on $\mathbb{R}^2$, with three hidden layers and randomly generated weights:

at (0,0) ;
at (5,0) ;
at (10,0) ;

Applications II: Invariant kernels.
We can similarly define $\group$-invariant reproducing kernels on $\mathbb{R}^n$, by starting with a
kernel $\hat{\kappa}$ on $\mathbb{R}^N$ and defining a function on $\mathbb{R}^n$ as

\begin{equation*}
\kappa(x,y)\;=\;\hat{\kappa}\circ(\rho\otimes\rho)(x,y)\;=\;\hat{\kappa}(\rho(x),\rho(y))\;.
\end{equation*}

This function is again a kernel. In <ref>, we
show that its reproducing kernel Hilbert space consists of continuous $\group$-invariant
functions on $\mathbb{R}^n$.
We also show that, even though $\mathbb{R}^n$ is not compact, $\kappa$ behaves essentially like a kernel on a
compact domain (<ref>). In particular,
it satisfies a Mercer representation and a compact embedding property, both of which usually require
compactness. This behavior is specific to kernels invariant under crystallographic groups, and
does not extend to more general groups of isometries on $\mathbb{R}^n$.

Applications III: Invariant Gaussian processes.
There are two ways in which a Gaussian process (GP) can be invariant under a group: A GP is a distribution on
functions, and we can either ask for each function it generates to be invariant, or only require that its distribution is
invariant (see <ref> for definitions). The former implies the latter. Both types of processes can be constructed
by factoring through an orbifold:
Suppose we start with a kernel $\hat{\kappa}$ (a covariance function) and a real-valued function $\hat{\mu}$ (the mean function),
both defined on $\mathbb{R}^N$. If we then generate a random function $F$ on $\mathbb{R}^n$ as

\begin{equation*}
F\;:=\;H\circ\rho\qquad\text{ where }\qquad H\,\sim\,\text{GP}(\hat{\mu},\hat{\kappa})\;,
\end{equation*}

the function $F$ is $\group$-invariant with probability 1.
The following are examples of such random functions, rendered as contour plots with non-smooth colormaps.

at (0,0) ;
at (3.5,0) ;
at (7,0) ;
at (10.5,0) ;

If we instead generate $F$ as

\begin{equation*}
F\,\sim\,\text{GP}(\mu,\kappa)\qquad\text{ where }\qquad \mu\;:=\;\hat{\mu}\circ\rho\;\text{ and }\;\kappa\;:=\;\hat{\kappa}\circ(\rho\otimes\rho)\;,
\end{equation*}

the distribution of $F$ is $\group$-invariant. See <ref>.

Properties of the Laplace operator.
<ref> studies differentials and Laplacians of crystallographically invariant functions
The results are then used in the proof of the Fourier representation.
Consider a vector field $F$, i.e., a function ${F:\mathbb{R}^n\rightarrow\mathbb{R}^n}$. An example
of such a vector field is the gradient ${F=\nabla f}$. <ref> shows that
the gradient transforms under elements $\phi$ of $\mathbb{G}$ as

\begin{equation*}
\nabla f(\phi x)\;=\;(\text{linear part of }\phi)\cdot\nabla f
\quad\text{ or abstractly }\quad
F\circ\phi\;=\;(\text{linear part of }\phi)\circ F\;.
\end{equation*}

<ref> shows that, for any vector field $F$ that transforms in this way,
the total flux through the boundary of the polytope $\Pi$ vanishes,

\begin{equation*}
\mint_{\partial\Pi}F(x)^{\trans}(\text{normal vector of $\partial\Pi$ at $x$})dx\;=\;0\;.
\end{equation*}

We can combine this fact with a result from the theory of partial differential equations, the
so-called Green identity, which decomposes the Laplacian on functions on $\Pi$ as

\begin{equation*}
-\Delta f\;=\;\text{self-adjoint component on interior of }\Pi\;-\;\text{correction term on }\partial\Pi\;.
\end{equation*}

<ref> makes the statement precise. Using the fact that the flux vanishes,
we can show that the correction term on $\partial\Pi$ vanishes, and from that deduce that the
Laplace operator on invariant functions is self-adjoint (<ref>).
That allows us to draw on
results from the spectral theory of self-adjoint operators to solve (<ref>).

Background and reference results.
Since our methods draw on a number of different fields,
the appendix provides additional background on groups of isometries (App. <ref>),
functional analysis (App. <ref>), and orbifolds (App. <ref>),
and spectral theory (App. <ref>).

§ PRELIMINARIES: CRYSTALLOGRAPHIC GROUPS

Throughout, we consider a Euclidean space $\mathbb{R}^n$, and write $d_n$ for Euclidean distance in $n$ dimensions.
Euclidean volume (that is, Lebesgue measure on $\mathbb{R}^n$) is denoted $\vol_n$.
As we work with both sets and their boundaries, we must carefully distinguish
dimensions: The span of a set ${A\subset\mathbb{R}^n}$ is the smallest affine subspace that contains it.
We define the dimension and relative interior of $A$ as

\begin{equation*}
\dim A\,:=\,\dim\text{span}\,A
\quad\text{ and }\quad
A^{\circ}\,:=\,\text{ largest subset of }A\text{ that is open in }\text{span}\,A\;.
\end{equation*}

The boundary of $A$ is the set ${\partial A:=A\setminus A^{\circ}}$.
If $A$ has dimension ${k<n}$, then ${\vol_k(A)}$ denotes Euclidean volume in $\text{span}\,A$.
For example: If ${A\subset\mathbb{R}^3}$ is a closed line segment, then ${\dim A=1}$, and
${\vol_1(A)}$ is the length of the line segment, whereas ${\vol_3(A)=\vol_2(A)=0}$.
Taking the relative interior $A^\circ$ removes
the two endpoints, whereas interior of $A$ in $\mathbb{R}^3$ is the empty set.
(No such distinction is required for the closure $\bar{A}$,
since $A$ is closed in $\text{span}\,A$ if and only if it is closed in $\mathbb{R}^n$.)

§.§ Defining crystallographic groups

Consider a group $\group$
of isometries of $\mathbb{R}^n$. (See <ref> for a
brief review of definitions.) Every isometry $\phi$ of $\mathbb{R}^n$ is of the form

\begin{equation}
\label{eq:euclidean:isometry}
\phi x\;=\;A_\phi x+b_\phi
\qquad\text{ for some orthogonal $n\times n$ matrix }A_\phi\text{ and some }b_\phi\in\mathbb{R}^n\;.
\end{equation}

Let ${M\subset\mathbb{R}^n}$ be a set. We say that $\group$
tiles the space $\mathbb{R}^n$ with $M$
if the image sets ${\phi M}$ completely cover the space so that only their boundaries overlap:

\begin{equation*}
\medcup\nolimits_{\phi\in\group}\phi M\;=\;\mathbb{R}^n
\qquad\text{ and }\qquad
\phi M\cap\psi M\;\subset\;\partial(\phi M)
\quad\text{ whenever }\phi\neq\psi\;.
\end{equation*}

Each set $\phi M$ is a tile,
and the collection ${\group M:=\braces{\phi M|\phi\in\group}}$ is a tiling of $\mathbb{R}^n$.

By a convex polytope, we mean the convex hull of a finite set of points
[Ziegler, 1995]. Let ${\cF\subset\mathbb{R}^n}$ be an $n$-dimensional convex polytope.
The boundary $\partial\cF$
consists of a finite number of ${(n-1)}$-dimensional convex polytopes, called the facets
of $\cF$. Thus, if $\group$ tiles $\mathbb{R}^n$ with $\Pi$, only points on
facets are contained in more than one tile.

A crystallographic group is a group of isometries that
tiles $\mathbb{R}^n$ with an $n$-dimensional convex polytope $\cF$.

The polytope $\cF$ is then also called a fundamental region (in geometry) or an
asymmetric unit (in materials science) for $\group$.
This definition of crystallographic groups
differs from those given in the literature, but we clarify in
<ref> that it is equivalent.

Tiling behavior of the 17 crystallographic groups on $\mathbb{R}^n$ (the wallpaper groups).

§.§ Basic properties

Some properties of $\group$ can be read right off the definition:
Since $\group$ tiles the entire space with a set $\cF$ of finite diameter, we must have ${|\group|=\infty}$.
Since $\Pi$ is $n$-dimensional and convex, it contains an open metric ball of positive radius.
Each tile contains a copy of this ball, and these copies do not overlap. It follows that

\begin{equation}
\label{discrete}
d(\phi(x),\psi(x))\;>\;\varepsilon\qquad\text{ for all distinct }\phi,\psi\in\group\text{ and all }x\in\Pi^{\circ}\;.
\end{equation}

A group of isometries that satisfies (<ref>)
for some ${\varepsilon>0}$ is called discrete, in contrast to groups which contain, e.g., continuous rotations.
Discreteness implies $\group$ is countable, but not all countable groups of isometries are discrete
(the group $\mathbb{Q}^n$ of rational-valued shifts is a non-example).
In summary, every crystallographic group is
an infinite, discrete (and hence countable) subgroup of the Euclidean group on $\mathbb{R}^n$.

Suppose we choose one of the tilings in <ref>,
and rotate or shift the entire plane with the tiling on it.
Informally speaking, that changes the tiling, but not the tiling mechanism, and it is natural to consider
the two tilings isomorphic. More formally,
two crystallographic groups $\group$ and $\group'$ are isomorphic if there is
an orientation-preserving, invertible, and affine (but not necessarily isometric) map
${\gamma:\mathbb{R}^n\rightarrow\mathbb{R}^n}$ such that
${\group'=\gamma\group}$, where ${\gamma\group:=\braces{\gamma\phi\gamma^{-1}\,|\,\phi\in\group}}$.
Up to isomorphy, there are only finitely many crystallographic groups on $\mathbb{R}^n$ for each ${n\in\mathbb{N}}$.
Specifically, there 17 such groups for ${n=2}$, and 230 for ${n=3}$.

§ PRELIMINARIES: INVARIANT FUNCTIONS

A function ${f:\mathbb{R}^n\rightarrow\xspace}$, with values in some set $\xspace$, is
$\phi$-invariant if it satisfies

\begin{equation*}
f(\phi x)=f(x)\quad\text{ for all }x\in\mathbb{R}^n\qquad\text{ or in short }\quad f\circ\phi=f\;.
\end{equation*}

It is $\group$-invariant if it is $\phi$-invariant
for all ${\phi\in\group}$.
We are specifically interested in $\group$-invariant functions that
are continuous, and write

\begin{equation*}
\C(M):=\braces{f:\mathbb{R}^n\rightarrow\mathbb{R}\,|\,\text{$f$ continuous}}
\quad\text{ and }\quad
\CG:=\braces{f\in\C(\mathbb{R}^n)\,|\,
f\text{ is $\group$-invariant}}\;.
\end{equation*}

More generally, a function
${f:(\mathbb{R}^n)^k\rightarrow\xspace}$ is $\group$-invariant in each argument if

\begin{equation}
\label{eq:invariant:k:arguments}
f(\phi_1x_1,\ldots,\phi_kx_k)\;=\;f(x_1,\ldots,x_k)\quad\text{ for all }\phi_1,\ldots\phi_k\in\group\text{ and }x_1,\ldots,x_k\in\mathbb{R}^n\;.
\end{equation}

§.§ Tiling with functions

To construct a $\group$-invariant function, we may start with a function $h$ on $\Pi$
and “replicate it by tiling”. For that to be possible,
$h$ must in turn be the restriction of a $\group$-invariant function to
$\Pi$. It must then satisfy ${h(\phi x)=h(x)}$ if both $\phi x$ and $x$ are in $\Pi$.
We hence define the relation

\begin{equation*}
\quad:\Longleftrightarrow\quad
x,y\in\cF\text{ and }y=\phi(x)\text{ for some }\phi\in\group\setminus\braces{\Id}\;.
\end{equation*}

We note immediately that ${x\sim y}$ implies each point is also contained in
an adjacent tile, so both must be on the boundary $\partial\Pi$ of $\Pi$.
The requirement

\begin{equation}
\label{pbc}
h(x)\;=\;h(y) \qquad\text{ whenever }x\sim y
\end{equation}

is therefore a periodic boundary condition.
If it holds, the function

\begin{equation}
\label{eq:tiled:function}
f(x)\;:=\;h(\phi^{-1} x)\qquad\text{ for }x\in\phi(\cF)\text{ and each }\phi\in\group
\end{equation}

is well-defined on $\mathbb{R}^n$, and is $\group$-invariant.
Conversely, every $\group$-invariant function $f$ can be obtained this way (by choosing $h$ as the
restriction $f|_{\Pi}$).
Informally, (<ref>) says that we stitch together function segments on tiles that are all copies of
$h$, and these segments overlap on the tile boundaries. The boundary condition ensures
that wherever such overlaps occur, the segments have the same value, so that
(<ref>) produces no ambiguities.
The special case of (<ref>) for pure shift groups—where
$A_\phi$ is the identity matrix for all ${\phi\in\group}$—is known as a
Born-von Karman boundary condition (e.g., Ashcroft and Mermin, 1976).

§.§ Orbits and quotients

An alternative way to express invariance is as follows: A function is $\group$-invariant
if and only if it is constant on each set of the form

\begin{equation*}
\group(x)\;:=\;\braces{\phi x \,|\,\phi\in\group}
\qquad\text{ for each }x\in\mathbb{R}^n\;.
\end{equation*}

The set $\group(x)$ is called the orbit of $x$.
We see immediately that each orbit of a crystallographic group is countably infinite, but
locally finite: The definition of discreteness in (<ref>)
implies that every bounded subset of $\mathbb{R}^n$
contains only finitely many points of each orbit.
We also see that each point ${x\in\mathbb{R}^n}$ is in one and only one orbit, which means
the orbits form a partition of $\mathbb{R}^n$. The assignment ${x\mapsto\group(x)}$ is
hence a well-defined map

\begin{equation*}
\quad\text{ with image }\quad
\mathbb{R}^n/\group\;:=\;q(\mathbb{R}^n)\;=\;\braces{\group(x)\,|\,x\in\mathbb{R}^n}\;.
\end{equation*}

The orbit set $\mathbb{R}^n/\group$ is also called the quotient set or just the quotient of
$\group$, and $q$ is called the quotient map (e.g., Bonahon, 2009).
Since the orbits are mutually disjoint, we can informally think of $q$ as collapsing each
orbit into a single point, and $\mathbb{R}^n/\group$ is the set of such points.

Left: A point $x$ and points on its orbit $\group(x)$ in a shift tiling of $\mathbb{R}^2$.
Right: The infimum in the definition of $d_\group(\group(x),\group(y))$ may be attained by the
Euclidean distance $d_n(x,y)$ between a point $y$
in $\Pi$ and a point $\phi x$ in a different tile $\phi\Pi$.

at (0,0) ;
at (7,0) ;

Quotient spaces are abstract but useful tools for expressing invariance properties:
For any function ${f:\mathbb{R}^n\rightarrow\mathbb{R}}$, we have

\begin{equation}
\label{eq:bijection:f:quotient}
f\text{ is $\group$-invariant }
\quad\Longleftrightarrow\quad
f=\hat{f}\circ q
\quad\text{ for some function }\hat{f}:\mathbb{R}^n/\group\rightarrow\mathbb{R}\;,
\end{equation}

since each point of ${\mathbb{R}^n/\group}$ represents an orbit and $f$ is invariant iff it is constant on orbits.
We can also use the quotient to express continuity, by equipping it with a topology that satisfies

\begin{equation}
\label{eq:bijection:f:quotient:continuous}
\quad\Longleftrightarrow\quad
f=\hat{h}\circ q
\quad\text{ for some continuous }\hat{h}:\mathbb{R}^n/\group\rightarrow\mathbb{R}\;.
\end{equation}

There is exactly one such topology, called the quotient topology in the literature. Its definition
can be made more concrete by metrizing it:
[see Bonahon, 2009, Theorem 7.7]

If $\group$ is crystallographic, the function

\begin{equation*}
d_\group(\omega_1,\omega_2)\,:=\,\inf\braces{d(x,y)\,|\,{x\in \omega_1,y\in \omega_2}}
\quad\text{ for }\omega_1,\omega_2\in\mathbb{R}^n/\group
\end{equation*}

is a valid metric on ${\mathbb{R}^n/\group}$, and it metrizes the quotient
topology. A subset ${U\subset\mathbb{R}^n/\group}$ is open if and only
if its preimage ${q^{-1}U}$ is open in $\mathbb{R}^n$.

Since $\group$ is discrete, the infimum in $d_\group$ is a minimum.
The distance of two orbits (considered as points in $\mathbb{R}^n/\group$)
is hence the shortest Euclidean distance between points in these orbits (considered as sets
in $\mathbb{R}^n$), see <ref> (right).
If $x$ and $y$ are points in the polytope $\Pi$, we have

\begin{equation*}
x\sim y\quad\Longleftrightarrow\quad d_{\group}(\group(x),\group(y))\;=\;0\;.
\end{equation*}

Informally speaking, $d_\group$ implements the periodic boundary condition (<ref>).
The metric space $(\mathbb{R}^n/\group,d_\group)$ is also called the quotient space or
orbit space of $\group$. A very important property of crystallographic groups is that they
have compact quotient spaces:[<cit.>]

If a discrete group $\group$ of isometries tiles $\mathbb{R}^n$ with a set $M$,
the quotient space $(\mathbb{R}^n/\group,d_\group)$ is homeomorphic to the quotient space $M/\group$.
If $\group$ is crystallographic and tiles
with a convex polytope, then $(\mathbb{R}^n/\group,d_\group)$ is compact.

§.§ Transversals and projections

Since orbit spaces are abstract objects, we can only work with them implicitly.
One way to do so is by representing each orbit by one of its points in $\mathbb{R}^n$.
A subset of $\mathbb{R}^n$ that contains exactly one point of each orbit is called a
transversal. In general, transversals can be exceedingly complex sets [Becker and Kechris, 1996],
but crystallographic
groups always have simple transversals. <ref> in the next section
constructs a transversal explicitly.
In the following, we will always write $\tilde{\Pi}$ to mean

\begin{equation*}
\tilde{\Pi}\;:=\;
\text{ a transversal contained in $\Pi$ computed by \cref{algorithm:transversal}.}
\end{equation*}

Given such a transversal, we can define the projector ${p:\mathbb{R}^n\rightarrow\Pi}$ as

\begin{equation}
\label{eq:projector}
p(x)\;:=\;\text{ the unique element of }\group(x)\cap\widetilde{\Pi}\;.
\end{equation}

If we think of each point in $\tilde{\Pi}$ as a concrete representative of an element of ${\mathbb{R}^n/\group}$,
then $p$ is similarly a concrete representation of the quotient map $q$, and we can translate
the identities above accordingly:
The projector is by definition $\group$-invariant,
since we can write $f$ in (<ref>) as ${f=h\circ p}$. That shows

\begin{equation}
\label{eq:bijection:Rn:Pi}
f:\mathbb{R}^n\rightarrow\mathbb{R}\text{ is $\group$-invariant }
\;\;\Longleftrightarrow\;\;
f\;=\;h\circ p \text{ for some }h:\cF\rightarrow\mathbb{R}\text{ satisfying }\eqref{pbc}\;.
\end{equation}

Although $p$ is not continuous as a function ${\mathbb{R}^n\rightarrow\Pi}$, continuity only fails at the boundary, and
$p$ behaves like a continuous function when composed with $h$:

Let ${h:\Pi\rightarrow\xspace}$ be a continuous function with values in a topological space $\xspace$.
If $h$ satisfies (<ref>), then ${h\circ p}$ is a continuous $\group$-invariant function
${\mathbb{R}^n\rightarrow\xspace}$. It follows that

\begin{equation}
\label{eq:projector:composition}
\;\;\Longleftrightarrow\;\;
f\;=\;h\circ p \quad\text{ for some continuous }h:\cF\rightarrow\mathbb{R}\text{ satisfying }\eqref{pbc}\;.
\end{equation}

Since $p$ exists for any choice of $\group$ and $\Pi$, and since it can be evaluated algorithmically,
we have hence reduced the problem of constructing continuous invariant functions to the problem
of finding functions that satisfy the periodic boundary condition (<ref>).

§ TAKING QUOTIENTS ALGORITHMICALLY: ORBIT GRAPHS

at (0,0) ;
at (4,0) ;
at (8,0) ;
at (12,0) ;

at (0,0) $\Pi$ within the tiling;
at (4,0) construct net;
at (8,0) apply elements of $\mathcal{A}_\Pi$;
at (12,0) orbit graph;

at (0,0) ;
at (4,0) ;
at (8,0) ;
at (12,0) ;

Orbit graph construction for two plane groups:  and . Starting with the fundamental region (left), place points to construct an $\epsilon$-net (middle left),
apply local group operations to these points (middle right), then add edges which may include vertices outside the fundamental region (right).

To work with invariant functions computationally, we
must approximate the quotient metric. We do so using a data structure that we call an orbit graph,
in which two points are connected if their orbits are close to each other. More formally, any undirected graph
is a metric space when equipped with path length as distance. The metric space defined by the graph
$\mathcal{G}$ below discretizes the metric space ${(\mathbb{R}^n/\group,d_\group)}$.
To define $\mathcal{G}$, fix constants ${\varepsilon,\delta>0}$.
A finite set $\Gamma$ is an $\varepsilon$-net in $\Pi$ if each point lies within distance $\varepsilon$ of $\Gamma$,

\begin{equation*}
\Gamma\text{ is an $\varepsilon$-net }
\quad:\Leftrightarrow\quad
\text{ for each }x\in\Pi\text{ there exists }z\in\Gamma\text{ such that }d_n(x,z)\,<\,\varepsilon\;,
\end{equation*}

see e.g., Cooper et al., 2000.
If $\Gamma$ is an $\varepsilon$-net in $\Pi$, we call the graph

\begin{equation*}
\mathcal{G}=\mathcal{G}(\varepsilon,\delta)=(\Gamma,E)
\quad\text{ where }\quad E:=\braces{(x,y)\,|\,d_\group(\group(x),\group(y)<\delta}
\end{equation*}

an orbit graph for $\group$ and $\Pi$.

§.§ Computing orbit graphs

Algorithmically, an orbit graph can be constructed as follows: Constructing an $\varepsilon$-net is a standard problem in
computational geometry and can be solved efficiently (e.g., Haussler and Welzl, 1987).
Having done so, the problem we have to solve is:

\begin{equation*}
\text{Given points }x,y\in\Pi\text{, determine }d_\group(\group(x),\group(y))\;.
\end{equation*}

Since $\Pi$ is a polytope, its diameter

\begin{equation*}
\text{diam}(\Pi)\;:=\;\max\braces{d_n(x,y)\,|\,x,y\in\Pi}\;<\;\infty
\end{equation*}

can also be evaluated computationally. By definition of $d_\group$, we have

\begin{equation*}
\;=\;
\min_{\phi\in\group}d_n(x,\phi y)
\;\leq\;
\;\leq\;
\text{diam}(\Pi)\;.
\end{equation*}

That shows the minimum is always attained for a point $\phi y$ on a tile $\phi\Pi$ that lies within
distance $\text{diam}(\Pi)$ of $x$. The set of transformations that specify these tiles is

\begin{equation*}
\mathcal{A}_\Pi
\;=\;
\braces{\phi\in\group \,|\,d_n(x,\phi z)\leq\text{diam}(\Pi)\text{ for some }z\in\Pi}\;.
\end{equation*}

This set is always finite, since $\group$ is discrete and the ball of radius $\text{diam}(\Pi)$ is compact.
We can hence evaluate the quotient metric as

\begin{equation*}
\;=\;
\min\braces{d_n(x,\phi y)\,|\,\phi\in\mathcal{A}_\Pi}\,,
\end{equation*}

which reduces the construction of $E$ to a finite search problem. In summary:
[Constructing the orbit graph]

1.)     Construct the $\varepsilon$-net $\Gamma$.

2.)     Find local group elements $\mathcal{A}_\Pi$.

3.)     For each pair $x,y\in\Gamma$, find $d_\group(\group(x),\group(y))=\min\braces{d_n(x,\phi y)\,|\,\phi\in\mathcal{A}_\Pi}$.

4.)     Add an edge between $x$ and $y$ if $d_\group(\group(x),\group(y)) < \delta$.

The construction is illustrated in <ref>.

§.§ Computing a transversal

Recall that the faces of a polytope are its vertices, edges, and so forth; the
facets are the ${(n-1)}$-dimensional faces. The polytope itself is also a face, of dimension $n$.
See [Ziegler, 1995] for a precise definition. Given $\Pi$ and $\group$, we will call two
faces $S$ and $S'$ $\group$-equivalent if ${S'=\phi S}$ for some ${\phi\in\group}$.
Thus, if ${S=\Pi}$, its equivalence class is $\braces{\Pi}$.
If $S$ is a facet, it is equivalent to at most one distinct facet, so its
equivalence class has one or two elements. The equivalence classes
of lower-dimensional faces may be larger—if $\group$ is  and $\Pi$ a square,
for example, all four vertices of $\Pi$ are $\group$-equivalent.
[Constructing a transversal]

1)     Start with an exact tiling. Enumerate all faces of $\Pi$.

2)     Sort faces into $\group$-equivalence classes.

3)     Select one face from each class and take its relative interior.

4)     Output the union $\tilde{\Pi}$ of these relative interiors.

The set $\tilde{\Pi}$ constructed by <ref> is a transversal.

The relative interiors of the faces of a convex polytope are mutually disjoint and
their union is $\Pi$, so each point ${x\in\Pi}$ is on exactly one such relative interior.
Let $S$ be the face with ${x\in S^{\circ}}$, and consider any ${\phi\in\group}$.
Since the tiling is exact, ${\phi S}$ is either a face of $\Pi$ or ${\phi S\cap\Pi=\varnothing}$.
If ${\phi x\in\Pi}$, the intersection cannot be empty, so $\phi S$ is a face and hence
$\group$-equivalent to $S$.
It follows that the interior of a face of $\Pi$ intersect the orbit $\group x$ if and only
if it is in the equivalence class of $S$.
Since we select exactly one element of this class, exactly one point of $\group x$
is contained in $\tilde{\Pi}$.

§.§ Computing the projector

Since $\group$ is crystallographic, it contains shifts in $n$ linearly independent directions, and these shifts hence specify
a coordinate system of $\mathbb{R}^n$. More precisely: There are $n$ elements ${\phi_1,\ldots,\phi_n}$ of
$\group$ that (1) are pure shifts (satisfy ${A_{\phi_i}=\mathbb{I}}$), (2) are linearly independent, and (3) are the shortest
such elements (in terms of the Euclidean norm of $b_{\phi_i}$). Up to a sign, each of these elements is uniquely determined.
We refer to the vectors ${\phi_1,\ldots,\phi_n}$ as the shift coordinate system of $\group$.
[Computing the projector]

1.)     Perform a basis change from the shift coordinates to the standard basis of $\mathbb{R}^n$.

2.)     Set ${\tilde{x}=(x_1\mod 1,\ldots,x_n\mod 1)}$.

3.)     Find ${\phi\in\group}$ such that ${\phi\tilde{x}\in\tilde{\Pi}}$.

4.)     Apply the reverse change of basis from standard to shift coordinates.

§ LINEAR REPRESENTATION: INVARIANT FOURIER TRANSFORMS

In this section, we obtain a basis representation for invariant functions: given a
crystallographic group $\group$, we construct a sequence of $\group$-invariant functions
${\ef_1,\ef_2,\ldots}$ on $\mathbb{R}^n$ such that any $\group$-invariant continuous function
can be represented as a (possibly infinite) linear combination ${\sum_{i\in\mathbb{N}}c_i\ef_i}$.
If $\group$ is generated by $n$ orthogonal shifts, the functions $e_i$
are an $n$-dimensional Fourier basis.
<ref> below obtains an analogous basis for each crystallographic group $\group$.

§.§ Representation theorem

For any open set ${M\subseteq\mathbb{R}^n}$, we define the Laplace operator on twice differentiable
functions ${h:M\rightarrow\mathbb{R}}$ as

\begin{equation*}
\Delta h
\;:=\;
\mfrac{\partial^2 h}{\partial x_1^2}+\ldots+\mfrac{\partial^2 h}{\partial x_n^2}
\;=\;
\nabla^{\trans}(\nabla h)\;.
\end{equation*}

Now consider specifically functions ${e:\mathbb{R}^n\rightarrow\mathbb{R}}$.
Fix some ${\lambda\in\mathbb{R}}$, and consider the constrained partial differential equation

\begin{equation}
\label{sturm:liouville}
\begin{aligned}
-\Delta e&\;=\;\lambda e\quad&&\text{ on }\mathbb{R}^n\\
\text{subject to }\qquad e&\;=\;e\circ\phi\quad&&\text{ for }\phi\in\group\;.
\end{aligned}
\end{equation}

Clearly, there is always a trivial solution, namely the constant function ${e=0}$.
If (<ref>) has a non-trivial solution $e$, we call this $e$ a $\group$-eigenfunction and
$\lambda$ a $\group$-eigenvalue
of the linear operator $-\Delta$. Denote the set of solutions by

\begin{equation*}
\mathcal{V}(\lambda)
\;:=\;
\braces{e:\mathbb{R}^n\rightarrow\mathbb{R}\,|\,e\text{ satisfies }\eqref{sturm:liouville}}\;.
\end{equation*}

Since $0$ is a solution, and any linear combination of solutions is again a solution,
${\mathcal{V}(\lambda)}$ is a vector space,
called the eigenspace of $\lambda$. Its dimension

\begin{equation*}
\end{equation*}

is the multiplicity of $\lambda$.

Let $\group$ be a crystallographic group that tiles $\mathbb{R}^n$ with
a convex polytope $\Pi$.
Then the constrained problem (<ref>) has solutions for countably many
distinct values ${\lambda_1,\lambda_2,\ldots}$ of $\lambda$, and these values satisfy

\begin{equation*}
\qquad\text{ and }\qquad
\lambda_i\;\xrightarrow{i\rightarrow\infty}\;\infty\;.
\end{equation*}

Every solution function $e$ is infinitely often differentiable.
There is a sequence ${e_1,e_2,\ldots}$ of solutions whose restrictions
${e_1|_\Pi,e_2|_\Pi,\ldots}$ to $\Pi$ form an orthonormal basis of the the space $\L_2(\Pi)$,
and satisfy

\begin{equation*}
\big|\braces{\,j\in\mathbb{N}\,|\,e_j\in\mathcal{V}(\lambda_i)}\big|\;=\;k(\lambda_i)
\qquad\text{ for each }i\in\mathbb{N}\;.
\end{equation*}

A function ${f:\mathbb{R}^n\rightarrow\mathbb{R}}$ is $\group$-invariant and continuous if and only if

\begin{equation*}
\qquad\text{ for some sequence }c_1,c_2,\ldots\in\mathbb{R}\;,
\end{equation*}

where the series converges in the supremum norm.

See <ref>.

The space ${\L_2(\mathbb{R}^n)}$ contains no non-trivial $\group$-invariant functions, since
for every ${f\in\C_\group}$

\begin{equation*}
\|f\|_{\L_2(\mathbb{R}^n)}\;=\;\tsum_{\phi\in\group}\|\,f|_{\phi\Pi}\,\|_{\L_2(\phi\Pi)}\;=\;
\begin{cases}
0 & f=0\text{ almost everywhere}\\
\infty &\text{otherwise}
\end{cases}\;.
\end{equation*}

On the other hand, the restriction $f|_\Pi$ is in $\L_2(\Pi)$, and completely determines $f$.
That makes $\L_2(\Pi)$ the natural $\L_2$-space in the context of crystallographic invariance,
and is the reason why the restrictions $e_i|_\Pi$ are used in the theorem.
Since $\L_2(\phi\Pi)$ is isometric to $\L_2(\Pi)$ for all ${\phi\in\group}$, it does not matter which tile we restrict to.

at (0,0) ;
at (3,0) ;
at (6,0) ;
at (9,0) ;
at (12,0) ;
at (0,-1.7) 39.48;
at (3,-1.7) 39.48;
at (6,-1.7) 39.48;
at (9,-1.7) 39.48;
at (12,-1.7) 78.95;

The first five non-constant basis functions ${e_2,\ldots,e_{6}}$
in <ref>, with their eigenvalues $\lambda_i$, for the group .
The eigenbasis in this case is precisely the standard Fourier basis on $\mathbb{R}^2$.

§.§ Relationship to Fourier series

at (0,0) ;
at (3,0) ;
at (6,0) ;
at (9,0) ;
at (12,0) ;
at (0,-1.7) 52.62;
at (3,-1.7) 157.87;
at (6,-1.7) 210.87;
at (9,-1.7) 368.02;
at (12,-1.7) 368.34;

at (0,0) ;
at (3,0) ;
at (6,0) ;
at (9,0) ;
at (12,0) ;
at (0,-1.7) 473.55;
at (3,-1.7) 631.41;
at (6,-1.7) 683.39;
at (9,-1.7) 684.06;
at (12,-1.7) 841.80;

The first ten non-constant basis functions ${e_2,\ldots,e_{11}}$
in <ref>, with their approximate eigenvalues $\lambda_i$, for the group .
In this case, the mulitplicities $k(\lambda_i)$ are not $2n=4$ as for the standard Fourier transform.

The standard Fourier bases for periodic functions on $\mathbb{R}^n$
can be obtained as the special cases of <ref>
for shift groups: Fix some edge width ${c>0}$, and choose $\Pi$ and $\mathbb{G}$ as

\begin{equation*}
\Pi\;=\;[0,c]^n
\qquad\text{ and }\qquad
\group=\bigl\{x\mapsto x+c(i_1,\ldots,i_n)^{\trans}\,\big|\,i_1,\ldots,i_n\in\mathbb{Z}\bigr\}\;.
\end{equation*}

For these groups, all eigenvalue multiplicities are ${k(\lambda_i)=2n}$ for each ${i\in\mathbb{N}}$.
For ${n=2}$, the group $\group$ is  (see <ref>). Its eigenfunctions
are shown in <ref>.

To clarify the relationship in more detail, consider the case
${n=1}$: Since $\Delta$ is a second derivative, the functions
${e(x)=\cos(\nu x)}$ and ${e(x)=\sin(\nu x)}$ satisfy

\begin{equation*}
\Delta e(x)\;=\;-\nu^2e(x)\qquad\text{ for each }\nu\geq 0\;,
\end{equation*}

and are hence eigenfunctions of ${-\Delta}$ with eigenvalue ${\lambda=\nu^2}$.
For this choice of $\Pi$ and $\group$, the invariance constraint
in (<ref>) holds iff ${e(x)=e(x+c)}$ for every ${x\in\mathbb{R}}$.
That is true iff

\begin{equation*}
\nu(x+c)\;=\;\nu x+2\pi(i-1)\quad\text{ for some $i\in\mathbb{N}$,
and hence }\quad
\lambda_i\;=\;\nu^2\;=\;\Bigl(\mfrac{2\pi(i-1)}{c}\Bigr)^2\;.
\end{equation*}

The eigenspaces are therefore the two-dimensional vector spaces

\begin{equation*}
\mathcal{V}(\lambda_i)\;=\;\text{span}\braces{\sin(\sqrt{\lambda_i}\argdot),\cos\sqrt{\lambda_i}\argdot)}
\quad\text{ with }\quad
\qquad\text{ for all }i\in\mathbb{N}\;.
\end{equation*}

Any continuous function $f$ that is $\group$-invariant (or, equivalently, $c$-periodic) can be expanded as

\begin{equation}
\label{eq:fourier:series:1d}
\end{equation}

In the notation of <ref>, the coefficients are
${c_{2i}=a_i}$ and ${c_{2i+1}=b_i}$, and

\begin{equation*}
\quad\text{ and }\quad
\;.
\end{equation*}

Note that the unconstrained equation has solutions for all $\lambda$ in the uncountable set ${[0,\infty)}$.
The invariance constraint limits possible values to the countable set ${\lambda_1,\lambda_2,\ldots}$.
If $f$ was continuous but not invariant, the expansion (<ref>) would hence require an integral
on the right. Since $f$ is invariant, a series suffices.

Fourier series, in particular in one dimension, are often written using complex-valued functions as

\begin{equation*}
f(x)\;=\;\msum_{i\in\mathbb{N}}\gamma_i\exp(J\lambda_i x)\qquad\text{ where }
\gamma_i\in\mathbb{C}\text{ and }J:=\sqrt{-1}\;.
\end{equation*}

Since Euler's formula ${\exp(Jx)=\cos(x)+J\sin(x)}$ shows

\begin{equation*}
\gamma_i\exp(J\sqrt{\lambda_i})
\;=\;
\quad\text{ for }(a_i-Jb_i)=\gamma_i\;,
\end{equation*}

that is equivalent to (<ref>). The complex plane $\mathbb{C}$ is not inherent to the Fourier
representation, but rather a convenient way
to parameterize the two-dimensional eigenspace $\mathcal{V}(\lambda_i)$.
For general crystallographic groups, the complex representation is less useful,
since the multiplicities $k(\lambda_i)$ may not be even, as can be seen in <ref>.

at (0,0) ;
at (3,0) ;
at (6,0) ;
at (9,0) ;
at (12,0) ;
at (0,-1.7) 78.94;
at (3,-1.7) 158.46;
at (6,-1.7) 237.30;
at (9,-1.7) 238.46;
at (12,-1.7) 323.08;

at (0,0) ;
at (3,0) ;
at (6,0) ;
at (9,0) ;
at (12,0) ;
at (0,-1.7) 395.37;
at (3,-1.7) 438.02;
at (6,-1.7) 489.34;
at (9,-1.7) 490.42;
at (12,-1.7) 576.43;

The first ten non-constant basis functions ${e_2,\ldots,e_{11}}$
in <ref>, with their approximate eigenvalues $\lambda_i$, for the group  on $\mathbb{R}^3$.

§.§ Spectral algorithms

The eigenfunctions in <ref> can be
approximated by eigenvectors of a suitable graph Laplacian of the orbit graph as follows.
We first compute an orbit graph ${\mathcal{G}=(\Gamma,E)}$ as described in
<ref>. We weight each edge $(x,y)$ of the graph by

\begin{align}
w(x,y) &=
\begin{cases}
\exp\Bigl(-\frac{d^2_\group(\group(x), \group(y))}{2\epsilon^2}\Bigr) &
\text{ if }(x,y)\in E\\
0 &\text{ otherwise}
\end{cases}
\;.
\label{eqn:edge-weights}
\end{align}

The normalized Laplacian of the weighted graph is

\begin{align}
L = \mathbb{I} - D^{-1}W \qquad\text{ where }W_{xy}=w(x,y)\,,
\label{eqn:graph-laplacian}
\end{align}

and ${D}$ is the diagonal matrix containing the sum of each row of $W$.
See e.g., Chung, 1997 for more on the matrix $L$.
Our estimates of the eigenvalues and -functions of $\Delta$ are the eigenvalues and eigenvectors of $L$,

\begin{equation*}
\hat{\lambda}_i\;:=\;\text{$i$th eigenvalue of }L
\quad\text{ and }\quad
\hat{e}_i\;:=\;\text{$i$th eigenvector of }L\,.
\end{equation*}

These approximate the spectrum of $\Delta$ in the sense that

\begin{equation*}
\lambda_i
\;\approx\;
\mfrac{2}{\epsilon}\,\hat{\lambda}_i
\qquad\text{ and }\qquad
\;\approx\;
\mfrac{2}{\epsilon}\,(\hat{e}_{i})_x
\qquad\text{ for }x\in\Gamma\;,
\end{equation*}

see Singer, 2006. Once an eigenvector $\hat{e}_i$ is computed, values of
$\ef_i$ at points ${x\not\in\Gamma}$ can be estimated using standard interpolation methods.
[Computing Fourier basis]

1.)     Construct the orbit graph $(\Gamma,E)$.

2.)     Compute the normalized Laplacian matrix $L$ according to (<ref>).

3.)     Compute eigenvectors $\hat{\ef}_i$ and eigenvalues $\hat{\lambda}_i$ of $L$.

4.)     Return eigenvalues and interpolated eigenfunctions.

Alternatively, the basis can be computed using a Galerkin approach, which
is described in <ref>. The functions in Figures <ref>, <ref>
and <ref> are computed using the Galerkin method.

The orbit graph automatically enforces the boundary condition (<ref>), since it measures distance in terms of $d_\group$.
The exception are group elements that are reflections, since these imply an additional property that the graph does not resolve:
If $\phi$ is a reflection over a facet $S$, $x$ a point on $S$ (and hence $\phi x=x$), and $f$ a $\phi$-invariant smooth function,
we must have ${\nabla f(x)=-\nabla f(\phi x)}$, and hence $\nabla f=0$ on $S$. In the parlance of PDEs, this is a Neumann boundary condition,
and can be enforced in several ways:

1) For each point ${x_j\in\Gamma}$ that is on $S$, add a point $x_{j'}$ to $\Gamma$ and
the edge $(x_j,x_j')$ to $E$. Then constrain each eigenvector $e_i$ in <ref> to satisfy ${e_{ij}=e_{ij'}}$.
This approach is common in spectral graph theory (e.g., Chung, 1997).

2) Alternatively, one may symmetrize the orbit graph: For vertext $x_j$ that is close to $S$, add its reflection ${x_{j'}:=\phi(x_j)}$ to
$\Gamma$. Now construct the edge set according to $d_\group$ using the augmented vertex set, and again constrain eigenvectors to
satisfy ${e_{ij}=e_{ij'}}$.
Either constrained eigenvalue problem can be solved using techniques of Golub, 1973.

§ NONLINEAR REPRESENTATION: FACTORING THROUGH AN ORBIFOLD

at (.1,0) ;
at (3.5,0) ;
at (6.5,0) ;
at (9.5,0) ;

at (.1,0) ;
at (3.5,0) ;
at (6.5,0) ;
at (9.5,0) ;

Visualizations of four-dimensional orbifold embeddings for wallpaper groups  (top) and  (bottom).  Compare the  embedding to the three-dimensional visualization in Figure <ref>.

We now generalize MacKay's construction, as sketched in the introduction, from shifts to
crystallographic groups. The construction defines a map

\begin{equation*}
\rho|_\Pi:\,\Pi\rightarrow\Omega\subset\mathbb{R}^N
\quad\text{ which in turn defines }\quad
\rho:=\rho_\Pi\circ p:\mathbb{R}^n\rightarrow\Omega\;.
\end{equation*}

In MacKay's case, $\Pi$ is an interval and ${\Omega\in\mathbb{R}^2}$ a circle.
The circle can be obtained from $\Pi$ by “gluing” the ends of the interval to each other.
To generalize this idea, we proceed as follows: Starting with the polytope $\Pi$,
we find any pair of points $x$ and $y$ on the same orbit of $\group$, and “bend” $\Pi$
so that we can glue $x$ to $y$. That results in a surface $\Omega$ in $\mathbb{R}^N$, where ${N\geq n}$
since we have bent $\Pi$.
If we denote the point on $\Omega$ that corresponds to ${x\in\Pi}$ by
${\rho|_\Pi(x)}$, we obtain the maps above. We first show how to implement this construction numerically,
and then consider its mathematical properties.
In mathematical terms, the surface $\Omega$ is an orbifold, a concept that generalizes
the notion of a manifold. The term $\group$-orbifold is made precise in
in <ref>, but can be read throughout this section as a surface in $\mathbb{R}^N$
that is “smooth almost everywhere”.

§.§ Gluing algorithms

The gluing algorithm constructs numerical approximations $\widehat{\rho}$ of $\rho$ and
$\widehat{\Omega}$ of $\Omega$. Here, $\widehat{\Omega}$ is a surface in $\widehat{N}$ dimensions,
where (as we explain below) ${\widehat{N}}$ may be larger than $N$.
As in the linear formulation of <ref>, we start with the orbit graph ${\mathcal{G}=(\Gamma,E)}$,
but in this case weight the edges to obtain a weighted graph

\begin{equation*}
\mathcal{G}_w\;=\;(\Gamma,E_w)
\qquad\text{ with }\qquad
\text{weight}(x,y):=\begin{cases}d_\group(x,y) & \text{if }(x,y)\in E\\
0 & \text{if }(x,y)\not\in E
\end{cases}\;.
\end{equation*}

The weighted graph provides approximate distances in quotient space.
The surface $\widehat{\Omega}$ is constructed from this graph by multidimensional scaling
(MDS) [Kruskal, 1964].
MDS proceeds as follows: Let $R$ be the matrix of squared geodesic distances, with entries

\begin{equation*}
R_{ij}\;=\;(\text{weighted path length from $x_i$ to $x_j$ in $\mathcal{G}_w$})^2\;.
\end{equation*}

Let ${0<\delta_1\leq\ldots\leq\delta_{|\Gamma|}}$ be the eigenvalues and ${v_1,\ldots,v_{|\Gamma|}}$ the eigenvectors of the matrix

\begin{align*}
\tilde{R} &= -\frac{1}{2}\Bigl(\mathbb{I} - \mfrac{1}{|\Gamma|}\mathbb{J}\Bigr)\,R\,\Bigl(\mathbb{I} - \mfrac{1}{|\Gamma|}\mathbb{J}\Bigr)
\qquad\text{ where }
\mathbb{I}=\text{diag}(1,\ldots,1)\text{ and }
\mathbb{J}={\scriptsize
\begin{pmatrix} 1 & \cdots & 1\\ \svdots & & \svdots\\ 1 & \cdots &1
\end{pmatrix}}\;.
\end{align*}

The embedding of each point ${x_i}$ in the $\varepsilon$-net $\Gamma$ is
then given by

\begin{equation*}
\widehat{\rho}(x_i)\;:=\;
\left(\begin{matrix}
\sqrt{\delta_{|\Gamma|}}\,v_{|\Gamma|,i}\\[.1em]
\sqrt{\delta_{|\Gamma|-1}}\,v_{|\Gamma|-1,i}\\
\vdots\\
\sqrt{\delta_{\smash{|\Gamma|-\widehat{N}}}}\,v_{|\Gamma|-\widehat{N},i}
\end{matrix}\right)\,.
\end{equation*}

The dimension $\widehat{N}$ is chosen to minimize error in the distances.
From ${\smash{\widehat{\rho}}(x_1),\ldots,\smash{\widehat{\rho}}(x_{|\Gamma|})}$, the surface $\smash{\widehat{\Omega}}$ and the map
$\smash{\widehat{\rho}|_\Pi}$ are obtained by interpolation.

at (0,0) ;
at (3,0) ;
at (6,0) ;
at (9,0) ;
at (12,0) ;

at (0,0) ;
at (3,0) ;
at (6,0) ;
at (9,0) ;
at (12,0) ;

Visualizations of five-dimensional orbifold embeddings for  (top) and  (bottom).

at (0,0)   ;
at (5.5,0) ;
at (11,0)  ;
at (0,2.2)  p1;
at (5.5,2.2)pm;
at (11,2.2) p2gg;

at (0,0)   ;
at (5.5,0) ;
at (11.5,-.7)  ;
at (0,-2.3)   p2mg;
at (5.5,-2.3) p4gm;
at (11,-2.3)  p3m1;

Visualization of several plane group orbifolds in three dimensions, computed as
embeddings of orbit graphs as described in <ref>. The mesh on each orbifold is the image of
the orbit graph.
Note that reflections result in boundaries for , , , and , and that in particular all of the
edges for  are boundaries.  The self-intersections visible for  are an artifact of visualizing in three dimensions—for a sufficiently high dimension, the surface does not self-intersect.

[Gluing with multidimensional scaling]

1.)     Construct the weighted orbit graph $\mathcal{G}_w$.

2.)     Compute the eigenvalues $\delta_i$ and eigenvectors $v_i$ of $\tilde{R}$.

3.)     Compute vertex embeddings ${\smash{\widehat{\rho}}(x_1),\ldots,\smash{\widehat{\rho}}(x_{|\Gamma|})}$.

4.)     Return interpolated vertex embeddings.

Once $\smash{\widehat{\rho}|_\Pi}$ can be computed, we can also compute

$\widehat{\rho}:=\widehat{\rho}\,|_\Pi\circ p$,

since the projector $p$ can be evaluated using <ref>.

The procedure satisfies two desiderata for constructing the orbifold map:
1) facets to be glued will be brought together, and 2) distances between interior points in $\Pi$ will be approximately preserved.
The embedding is unique up to isometric transformations.
The embedding step is similar to the Isomap [Tenenbaum et al., 2000] algorithm, but
unlike Isomap embeds into a higher-dimensional space rather than a lower-dimensional one.

§.§ Example: Invariant neural networks

Given $\group$ and $\Pi$, compute $\widehat{\rho}$ and $\widehat{\Omega}$ using <ref>.
Choose a neural network

\begin{equation*}
\quad\text{ with parameter vector }\theta\text{ and set }\quad
\end{equation*}

Then $f_\theta$ is a real-valued neural network on $\mathbb{R}^n$.
<ref> shows examples of $f_\theta$ for ${n=3}$, where
$h_{\theta}$ has three hidden layers of ten units each, with
rectified linear (relu) activations, although the input dimension $\widehat{N}$ may vary according to the choice of
$\group$ and $\Pi$. The parameter vector is generated at random.

Since most ways of performing interpolation in the construction of $\smash{\widehat{\rho}}$ are
amenable to automatic differentiation tools, this representation is easy to incorporate into machine learning pipelines.
Moreover, universality results for neural networks (e.g., Hornik et al., 1989) carry over:
If a class of neural networks $h_\theta$ approximates to arbitrary precision in ${\C(\mathbb{R}^N)}$, the
the resulting functions $f_\theta$ approximate to arbitrary precision in $\C_\group$ (though the approximation rate may change under composition with $\rho$). See <ref>.

at (0,0) ;
at (3.8,0) ;
at (7.7,0) ;
at (11.6,0) ;
at (0,-2.3)  ;
at (3.8,-2.3);
at (7.7,-2.3) ;
at (11.6,-2.3) ;

Four random neural networks for four different space groups, each with three hidden layers of ten units each and relu activations.

Dirichlet domains. Left: The set $R_{\phi}(x)$, shown in gray,
is the closed half space of all points at least as
close to $x$ as to $\phi x$. (In machine learning terms, the boundary of $R_\phi(x)$
orthogonally intersects the connecting line between $x$ and $\phi x$ at its center, and is hence
the support vector classifier with two support vectors $x$ and $\phi x$.) Right: A Dirichlet domain
$\mathbf{D}(x)$ defined by the group generated by a vertical shift $\phi$ and a diagonal shift $\psi$ in the plane.

[gray!40!white,fill] (-2.5,-1.5) rectangle (2.5,0);
[circle,scale=.3,fill,label=right:$x$] (a) at (0,-1) ;
[circle,scale=.3,fill,label=right:$\phi x$] (b) at (0,1) ;
[dotted] (a)–(b);
[thick] (-2.5,0)–(2.5,0);
at (-1.8,-.8) $R_\phi(x)$;

[circle,scale=.3,fill,label=below:$x$] (a) at (0,0) ;
(0:) in 60,120,...,360   – (:) ;

[inner sep=1pt,circle] at (:) ;
[circle,scale=.3,fill,label=right:$\phi x$] at ($(a.30)+(30:1.732*\R)$) ;
[circle,scale=.3,fill,label=right:$\psi x$] at ($(a.90)+(90:1.732*\R)$) ;
[circle,scale=.3,fill,label=left:$\phi^{-1}\psi x$] at ($(a.150)+(150:1.732*\R)$) ;
[circle,scale=.3,fill,label=left:$\phi^{-1} x$] at ($(a.210)+(210:1.732*\R)$) ;
[circle,scale=.3,fill,label=left:$\psi^{-1} x$] at ($(a.270)+(270:1.732*\R)$) ;
[circle,scale=.3,fill,label=right:$\phi\psi^{-1} x$] at ($(a.330)+(330:1.732*\R)$) ;

§.§ Exact tilings

Although the properties of general orbifolds constitute one of the more demanding problems
of modern mathematics, orbifolds of crystallographic groups are particularly well-behaved, and are well-understood.
That we can draw directly on this theory is due to the fact that it uses a notion of gluing very similar to
that employed by our algorithms as a proof technique [Bonahon, 2009, Ratcliffe, 2006].
The two notions align under an additional condition:
A convex polytope $\Pi$ is exact for $\group$ if $\group$ tiles with $\Pi$, and if each
face $S$ of $\Pi$ can be represented as

\begin{equation*}
S\;=\;\Pi\cap\phi\Pi\qquad\text{ for some }\phi\in\group\;.
\end{equation*}

Not every $\Pi$ with which $\group$ tiles is exact—in <ref>, for example, the polytopes shown
for  and  are not exact, though all others are. However, given $\Pi$ and $\group$, we can always
construct an exact surrogate as follows:
Choose any point ${x\in\mathbb{R}}$ that is not a fixed point for any ${\phi\in\group\setminus\braces{\Id}}$.
If $\group$ is crystallographic, that is true for every point in the interior of $\Pi$. For each ${\phi\in\group}$, the set

\begin{equation*}
\;:=\;
\braces{y\in\mathbb{R}^n\,|\,d_n(y,x)\leq d_n(y,\phi x)}\;,
\end{equation*}

is a half-space in $\mathbb{R}^n$ (see <ref>/left). The intersection

\begin{equation*}
\mathbf{D}(x)
\;:=\;
\medcap\nolimits_{\phi\in\group} R_\phi(x)
\;=\;
\medcap\nolimits_{\phi\,|\,\phi\Pi\cap\Pi\neq\varnothing} R_\phi(x)
\end{equation*}

of these half-spaces is called a Dirichlet domain for $\group$ (<ref>/right).
If $\group$ is crystallographic,
$\mathbf{D}(x)$ is an exact convex polytope for $\group$.

For illustration, consider the group : We start with a rectangle $\Pi$.
The group is generated by two glide reflections $\phi$ and $\psi$, each of which shifts
$\Pi$ horizontally and then reflects it about one of its long edges (<ref>/left).
Exactness fails because the set ${\Pi\cap\phi\Pi}$, marked in black, is not a complete edge of $\Pi$.
A Dirichlet domain for this tiling differs significantly from $\Pi$ (<ref>/right).
Although substituting $\mathbf{D}(x)$ for $\Pi$ changes the look of the tiling,
it does not change the group—that is, we still work with the same set of transformations (rather than
another group in the same isomorphism class), and the axes of reflections are still
defined by the faces of $\Pi$ rather than those of $\mathbf{D}(x)$.

The Dirichlet domain of a tiling may differ from the convex polytope defining the tiling.
Left: A group generated by two glide reflections $\phi$ and $\psi$ that
tiles with a rectangle $\Pi$. Both glides shift the rectangle
and then reflect it about one of its long edges. The set ${\Pi\cap\phi\Pi}$, marked black,
is not a face of $\Pi$. Right: The Dirichlet domain $\mathbf{D}(x)$. When $\group$ tiles with
$\mathbf{D}(x)$, the axes of reflection are still defined by the edges of $\Pi$.

[gray] (0,1.5) rectangle (2,.5);
[gray] (0,-1.5) rectangle (2,-.5);
[gray] (0,-1.5) rectangle (-2,-.5);
[gray] (0,1.5) rectangle (-2,.5);
[gray] (-1,-.5) rectangle (-3,.5);
[gray] (1,-.5) rectangle (3,.5);
[gray] (1,-.5) rectangle (-1,.5);
[very thick] (0,.5)–(1,.5);
[circle,scale=.3,fill,label=left:$x$] (a) at (.5,0) ;
[circle,scale=.3,fill,label=left:$\phi x$] at ($(a)+(1,1)$) ;
[circle,scale=.3,fill,label=left:$\psi x$] at ($(a)+(1,-1)$) ;
[circle,scale=.3,fill,label=left:$\psi^{-1} x$] at ($(a)+(-1,1)$) ;
[circle,scale=.3,fill,label=left:$\phi^{-1} x$] at ($(a)+(-1,-1)$) ;

[circle,scale=.3,fill,label=below:$x$] (a) at (.5,0) ;
[circle,scale=.3,fill,label=right:$\phi x$] at ($(a)+(1,1)$) ;
[circle,scale=.3,fill,label=right:$\psi x$] at ($(a)+(1,-1)$) ;
[circle,scale=.3,fill,label=left:$\psi^{-1} x$] at ($(a)+(-1,1)$) ;
[circle,scale=.3,fill,label=left:$\phi^{-1} x$] at ($(a)+(-1,-1)$) ;

§.§ Properties of embeddings

<ref> can be interpreted as computing a numerical approximation $\widehat{\rho}$
to a “true” embedding map $\rho$, namely the map in (<ref>) in the introduction.
Our main result on the nonlinear representation, <ref> below,
shows that this map indeed exists for every crystallographic group, and describes some of
its properties. The proof of the theorem shows that $\rho$ and the set $\Omega$ can be constructed
by the following abstract gluing algorithm.

Abstract gluing construction.

1.)     Glue: Identify each ${x\in\partial\Pi}$ with the unique point ${y\in\partial\Pi}$ satisfying ${x\sim y}$.

2.)     Equip the glued set $M$ with metric $d_\group$.

3.)     Embed the metric space $(M,d_\group)$ as a subset ${\Omega\subset\mathbb{R}^N}$ for some ${N\in\mathbb{N}}$.

4.)     For each ${x\in\Pi}$, define $\rho|_{\Pi}(x)$ as the representative of $x$ on $\Omega$.

5.)     Set ${\rho:=\rho|_\Pi\circ p}$.

Since $\Pi$ contains at least one point of each orbit, and the gluing step identifies all points identifies all points
on the same orbit with each other, the glued set $M$ can be regarded as the quotient set ${\mathbb{R}^n/\group}$.
Recall that an embedding is a map ${M\rightarrow\Omega\subset\mathbb{R}^N}$ that is a homeomorphism
(a continuous bijection with continuous inverse) of the metric spaces $(M,d_\group)$ and $(\Omega,d_N)$.

The state the theorem, we need one additional bit of terminology: The stabilizer of $x$ in $\group$ is the set of all
$\phi$ that leave $x$ invariant,

\begin{equation*}
\Stab(x)\;:=\;\braces{\phi\in\group\,|\,\phi x=x}\;
\end{equation*}

see Vinberg and Shvartsman, 1993, Ratcliffe, 2006, Bonahon, 2009.
We explain the role of the stabilizer in more detail in the next subsection.

Let $\group$ be a crystallographic group that tiles $\mathbb{R}^n$ with an exact convex
polytope $\Pi$. Then the set $M$ constructed by
gluing is a compact $\group$-orbifold that
is isometric to $\mathbb{R}^n/\group$. This orbifold can be
embedded into $\mathbb{R}^N$ for some

\begin{equation*}
\end{equation*}

that is, there is compact subset ${\Omega\subset\mathbb{R}^N}$ such that the
metric space ${(\Omega,d_N)}$ is homeomorphic to ${(\mathbb{R}^n/\group,d_\group)}$.
In particular, every point ${x\in\Pi}$ is represented by one and only one point $\rho_{\Pi}(x)$.
We can hence define a map

\begin{equation*}
\rho:\mathbb{R}^n\rightarrow\Omega\subset\mathbb{R}^N
\qquad
\text{ as }
\qquad
\rho(x)\;:=\;\rho_{\Pi}(p(x))\;.
\end{equation*}

The map $\rho$ is continuous, surjective, and $\group$-invariant.
A function ${f:\mathbb{R}^n\rightarrow Y}$, with values
in some topological space $Y$, is $\group$-invariant and continuous if and only if

\begin{equation*}
f=h\circ\rho\qquad\text{ for some continuous }h:\mathbb{R}^N\rightarrow Y\;.
\end{equation*}

$\Omega$ is smooth almost everywhere, in the sense that

\begin{equation*}
\vol_n\braces{x\in\Pi\,|\,\Omega\cap B_{\varepsilon}(\rho(x))\text{ is not a manifold for any }\varepsilon>0}\;=\;0
\end{equation*}

where $B_\varepsilon(z)$ denotes the open Euclidean metric ball of radius $\varepsilon$ centered at ${z\in\mathbb{R}^N}$.

See <ref>.

(a) Note carefully what the theorem does and does not show about the
embedding algorithm in <ref>: It does say that the
glued set constructed by the algorithm discretizes an orbifold, and that an $N$-dimensional
embedding of this orbifold exists. It does
not show that the embedding computed by MDS matches this dimension—indeed, since
MDS attempts to construct an embedding that is also isometric (rather than just homeomorphic), we
must in general expect the MDS embedding dimension to be larger, and we have at present no proof
that an isometric embedding always exists.

(b) If the tiling defined by $\group$ and $\Pi$ is not exact,
we can nonetheless define an embedding $\rho$ that represents continuous functions that are invariant
functions with respect to this tiling: Construct a Dirichlet domain $\mathbf{D}$, and then construct $\rho$ by applying the gluing
algorithm to $\mathbf{D}$. Functions constructed as ${h\circ\rho}$ are then invariant for the tiling $(\group,\Pi)$.

We have now seen different representations of continuous $\group$-invariant functions on $\mathbb{R}^n$, respectively
by continuous functions on $\Pi$, on the abstract space $\mathbb{R}^n/\group$, and on $\Omega$. On $\Pi$, we must explicitly
impose the periodic boundary condition, so we are using the set

\begin{equation*}
\C_{\text{\rm pbc}}(\Pi)\;:=\;\braces{\hat{f}\in\C(\Pi)\,|\,\hat{f}\text{ satisfies \eqref{pbc}}}\;.
\end{equation*}

In these representations, the projector $p$, the quotient map $q$,
and the embedding map $\rho$ play very similar roles.
We can make that observation more rigorous:

Given a crystallographic group $\group$ that tiles with a convex polytope $\Pi$, consider the maps

\begin{align*}
I_\Pi:\C_{\text{\rm pbc}}(\Pi)&\rightarrow\C_\group
&\text{ and }&&
&\text{ and }&&
\hat{f}&\mapsto\hat{f}\circ p
\hat{g}&\mapsto\hat{g}\circ q
\hat{h}&\mapsto\hat{h}\circ\rho
\end{align*}

where $I_\Omega$ is only defined if $\Pi$ is exact. Equip all spaces with the supremum norm.
Then $I_\Pi$ and $I_{\mathbb{R}^n/\group}$ are isometric isomorphisms, and if $\Pi$
is exact, so is $I_\Omega$. In particular, $\C_\group$ is always a separable Banach space.

By <ref>, (<ref>) and <ref>, all three maps are bijections.
We also have

\begin{equation*}
\|\hat{f}\|_{\sup}
\;=\;
\sup\nolimits_{x\in\Pi}|\hat{f}(x)|
\;=\;
\sup\nolimits_{x\in\mathbb{R}^n}|\hat{f}(p(x))|
\;=\;
\|\hat{f}\circ p\|_{\sup}
\quad\text{ for }\hat{f}\in\C_{\text{\rm pbc}}(\Pi)\;,
\end{equation*}

and the same holds mutatis mutandis on $\mathbb{R}^n/\group$ and $\Omega$, so all
maps are isometries. Since $\Pi$ is compact, $\C(\Pi)$ is separable
<cit.>. The same hence holds for the closed subspace $\C_{\text{\rm pbc}}(\Pi)$,
and by isometry for $\C_\group$.

§.§ Why the glued surface may not be smooth

Fixed points become non-smooth under gluing. Left: A triangular polytope $\Pi$ is transformed by a rotation $\phi$, which rotates by $120^\circ$ around the vertex
$x$ of $\Pi$. Note that ${\phi^2=\phi^{-1}}$ and ${\phi^3=\Id}$.
Middle: The gluing defined by $\phi$ identifies the edge $(x,z)$ and the edge $(x,y)$.
Right: The resulting glued surface is a cone, with $x$ mapped to the tip. Locally around $x$, the cone is not a manifold.

Whether or not the glued surface is smooth depends on whether the transformations in $\group$ leave any points invariant. It is a known fact in geometry (and made precise in the proof of <ref>) that

\begin{align*}
\text{ glued surface is a manifold }
\qquad&\Leftrightarrow\qquad
\phi x\;\neq\;x\quad\text{for all }\phi\in\group\setminus\braces{\Id}\text{ and }x\in\mathbb{R}^n\;.\\
\intertext{That can be phrased in terms of the stabilizer as}
\text{ glued surface is a manifold }
\qquad&\Leftrightarrow\qquad
\Stab(x)\;=\;\braces{\Id}\qquad\text{ for all }x\in\mathbb{R}^n\;.
\end{align*}

It is straightforward to check that $\Stab(x)$ is a group [Vinberg and Shvartsman, 1993].
Since each $\phi$ is an isometry,
and shifts of $\mathbb{R}^n$ have no fixed points, $\phi(x)=x$ can only hold if ${b_\phi=0}$.
Thus, $\Stab(x)$ is always a subset of the point group $\group_o$ (in the terminology of <ref>),
which means it is finite. To illustrate its effect on the surface, consider the following examples.

(a) Recall that MacKay's construction [MacKay, 1998], as sketched in the introduction,
can be translated to crystallographic groups by setting ${\Pi=[0,1]}$ and choosing $\group$ as shifts. In this case,
${\Stab(x)=\braces{\Id}}$ for each ${x\in\mathbb{R}}$, and the glued surface is a circle, which is indeed a manifold.
The two-dimensional analogue is to choose ${\Pi=[0,1]^2}$ and $\group$ as the group  in <ref>,
in which case the glued surface is a torus as shown in <ref>, and hence again a manifold.

(b) Now suppose $\Pi$ is a triangle, $x$ one of its corners, and $\phi$ a $120^\circ$ rotation
around $x$, as illustrated in <ref>. Then ${\Stab(x)=\braces{\Id,\phi,\phi^2=\phi^{-1}}}$, and the glued surface ${\Omega=\rho(\Pi)}$ is a cone with $\rho(x)$ as its tip.
That means $\Omega$ is not a manifold, because no neighborhood of the tip can be mapped isometrically to a neighborhood in $\mathbb{R}^2$.

§ INVARIANT KERNELS

at (-.5,0) ;
at (3.2,0) ;
at (7,0) ;
at (10.5,0) ;

Invariant kernel functions using the kernel in (<ref>),
for groups on $\mathbb{R}^2$ and $\mathbb{R}^3$:  with $\ell=0.1$,  with $\ell=0.01$,
with $\ell=0.01$, and  with $\ell=0.05$ (from left to right).

Throughout this section, ${\kappa:\mathbb{R}^n\times\mathbb{R}^n\rightarrow\mathbb{R}}$ is a kernel, i.e., a
positive definite function, and $\mathbb{H}$ is its reproducing kernel Hilbert space, or RKHS.
<ref> reviews definitions.
We consider kernels that are $\group$-invariant in both arguments in the sense of (<ref>),
that is,

\begin{equation*}
\label{invariant:kernel}
\kappa(\phi x,\psi y)\;=\;\kappa(x,y)\qquad\text{ for all }\phi,\psi\in\group\text{ and all }x,y\in\mathbb{R}^n\;.
\end{equation*}

That is the natural notion of invariance for most purposes, since
such kernels are precisely those that define spaces of $\group$-invariant

If and only if $\kappa$ is $\group$-invariant in each argument,
all functions ${f\in\mathbb{H}}$ are $\group$-invariant.
If $\kappa$ is also continuous, all ${f\in\mathbb{H}}$ are continuous, and hence

<ref> implies that, to define an invariant kernel, we can start with any
kernel $\hat{\kappa}$ on the embedding space $\mathbb{R}^N$, and compose it with the embedding
map $\rho$:

Let $\hat{\kappa}$ be a kernel on ${\mathbb{R}^N}$. Then the function

\begin{equation*}
\kappa(x,y):=\hat{\kappa}(\rho(x),\rho(y))
\quad\text{ or in short }\quad
\kappa=\hat{\kappa}\circ(\rho\otimes\rho)
\end{equation*}

is a kernel on ${\mathbb{R}^n}$ that is $\group$-invariant
in both arguments. If $\hat{\kappa}$ is continuous, so is $\kappa$.

That follows immediately from <ref> and the fact that
the restriction of a kernel to a subset is again a kernel [Steinwart and Christmann, 2008].

Suppose $\hat{\kappa}$ is an radial basis function (RBF) kernel with length scale $\ell$ on $\mathbb{R}^N$,
and hence of the form ${\hat{\kappa}(z,z')=\exp(-\|z-z'\|^2/\ell^2)}$. Then $\kappa$ is simply

\begin{equation*}
\label{eqn:covariance-kernel}
\kappa(x,y)
\;=\;
\hat{\kappa}(\rho(x),\rho(z))
\;=\;
\exp\Bigl(-\frac{||\rho(x)-\rho(y)||^2}{2\ell^2}\Bigr)\;.
\end{equation*}

<ref> illustrates this kernel the two-dimensional groups
and  and the three-dimensional groups  and .

Once we have constructed an invariant kernel, its application to machine learning problems is
straightforward. That becomes obvious if we define
${\Phi(x)\;=\;\kappa(x,\argdot)}$, often called the feature map of $\kappa$
[Steinwart and Christmann, 2008].
Using the definition of the scalar product on $\mathbb{H}$ and the reproducing property (see <ref>),
we then have

\begin{equation*}
\Phi:\mathbb{R}^n\rightarrow\mathbb{H}
\qquad\text{ and }\qquad
\kappa(x,y)\;=\;\sp{\Phi(x),\Phi(y)}_{\mathbb{H}}\;.
\end{equation*}

If $\kappa$ is $\group$-invariant, then $\Phi$ is also $\group$-invariant by construction.
Recall that most kernel methods in machine learning are derived by substituting
a Euclidean scalar product by $\sp{\Phi(x),\Phi(y)}_{\mathbb{H}}$, thereby making a linear
method nonlinear. Using a $\group$-invariant kernel results in a $\group$-invariant
[Invariant SVM]
A support vector machine (SVM) with kernel $\kappa$ is determined by two finite sets of points $\mathcal{X}$ and
$\mathcal{Y}$ in $\mathbb{R}^n$. To train the SVM, one maps these points into $\mathbb{H}$ via $\Phi$, finds the shortest
connecting line between the convex hulls of $\Phi(\mathcal{X})$ and $\Phi(\mathcal{Y})$, and determines
a hyperplane $F$ that is orthogonal to this line and intersects its center—equivalently, in dual formulation, the
unique hyperplane that separates the convex hulls of $\Phi(\mathcal{X})$ and $\Phi(\mathcal{Y})$
and maximizes the $\mathbb{H}$-norm distance to both. The set of points $x$ in $\mathbb{R}^n$ whose
image $\Phi(x)$ lies on $F$ is the decision surface of the SVM in $\mathbb{R}^n$.
The hyperplane can be specified by two functions $g$ (an offset vector) and $h$ (a normal vector) in $\mathbb{H}$:
A function ${f\in\mathbb{H}}$ lies on $F$ if and only if

\begin{equation*}
\sp{f-g,h}_{\mathbb{H}}
\;=\;
\qquad\text{ or equivalently }\qquad
\sp{f,h}_{\mathbb{H}}
\;=\;
\sp{g,h}_{\mathbb{H}}\;.
\end{equation*}

Let $x$ be a point in $\mathbb{R}^n$. If $y$ and $z$ are points with ${g=\Phi(y)}$ and ${h=\Phi(z)}$, then

\begin{equation*}
x\text{ is on decision surface }\qquad\Longleftrightarrow\qquad\kappa(x,z)\;=\;\kappa(y,z)\;.
\end{equation*}

Since invariance of $\kappa$ implies ${\kappa(\phi x,z)=\kappa(x,z)}$, that shows the decision surface is $\group$-invariant.
<ref> shows examples.
In these figures the data were randomly generated with regions assigned labels using a random function generated as in <ref>.
The support vectors are highlighted and illustrate the effects of symmetry constraints: the decision surface can be determined by data observed far away.

Two of the most important results on kernels are Mercer's theorem
and the compact inclusion theorem <cit.>. The latter shows the inclusion map
${\mathbb{H}\hookrightarrow\C}$ is compact, and is used in turn to establish good
statistical properties of kernel methods, such as oracle inequalities and finite covering numbers
[Steinwart and Christmann, 2008]. Both results assume that $\kappa$ has
compact support. If $\kappa$ is invariant under a crystallographic group,
its support is necessarily non-compact, but the next result shows that
versions of both theorems hold nonetheless:

If $\kappa$ is continuous and $\group$-invariant in both arguments, the
inclusion map ${\mathbb{H}\hookrightarrow\C_{\group}}$ is compact.
There exist functions ${f_1,f_2,\ldots\in\mathbb{H}}$ and
scalars ${c_1\geq c_2\geq\ldots>0}$ such that

\begin{equation*}
\kappa(x,y)\;=\;\tsum_{i\in\mathbb{N}}c_if_i(x)f_i(y)\qquad\text{ for all }x,y\in\mathbb{R}^n\;,
\end{equation*}

and the scaled sequence ${(\sqrt{c_i}f_i)}$ is an orthonormal basis of $\mathbb{H}$.
With this basis,

\begin{equation*}
\mathbb{H}\;=\;
\braces{\;f\!=\!\tsum_{i\in\mathbb{N}}a_i\sqrt{c_i}f_i\,|\,
a_1,a_2,\ldots\in\mathbb{R}\text{ with }\tsum_{i}|a_i|^2<\infty}\;,
\end{equation*}

where each series converges in $\mathbb{H}$ and hence (by compactness of inclusion)
also uniformly.

Intuitively, that is the case because every $\group$-invariant kernel is the pullback
of a kernel on $\Omega$, and $\Omega$ is compact.
<ref> shows an application of such a kernel to generate a two-class classifier
with an $\group$-invariant decision surface.

at (0,0) ;
at (5,0) ;
at (10,0) ;

Support vector machine decision surfaces for $\group$-invariant kernels constructed as in
<ref>. The groups used are  (left),  (middle),
and  (right). Support vectors are highlighted with red circles.

§ INVARIANT GAUSSIAN PROCESSES

We now consider the problem of generating random functions
${F:\mathbb{R}^n\rightarrow\mathbb{R}}$ such that each instance of $F$ is
continuous and $\group$-invariant with probability 1. That can be done linearly
using the generalized Fourier representation, by generating the coefficients $c_i$
in <ref> at random. Here, we consider the nonlinear representation
instead: If we set

\begin{equation*}
F\;:=\;H\circ\rho\qquad\text{ for a random continuous function }H:\mathbb{R}^N\rightarrow\mathbb{R}\;,
\end{equation*}

<ref> implies that $F$ is indeed continuous and $\group$-invariant with probability $1$,
and hence a random element of
$\C_{\group}$. Conversely, the result also implies that every random
element of $\C_{\group}$ is of this form, for some random element $H$ of $\C(\mathbb{R}^N)$.

§.§ Almost surely invariant processes

Recall that a random function ${F:M\subseteq\mathbb{R}^n\rightarrow\mathbb{R}}$ is a
Gaussian process if
the joint distribution of the random vector ${(F(x_1),\ldots,F(x_k))}$ is
Gaussian for any finite set of points ${x_1,\ldots,x_k\in M}$.
The mean and covariance function of a Gaussian process are defined as

\begin{equation*}
\mu(x)\;:=\;\mean[F(x)]
\quad\text{ and }\quad
\kappa(x,y)\;:=\;\mean[(F(x)-\mu(x))(F(y)-\mu(y)]
\quad\text{ for }x,y\in M\;.
\end{equation*}

The covariance function is always positive definite, and hence a kernel on
$M$. The distribution of a Gaussian process is completely determined by $\mu$ and $\kappa$,
and conditions for $F$ to satisfy continuity or stronger regularity conditions can
be formulated in terms of $\kappa$. See e.g., Marcus and Rosen, 2006 for more

Let $H$ be a continuous Gaussian process on $\mathbb{R}^N$, with mean $\mu$ and covariance function $\kappa$.
Then ${F:=H\circ\rho}$ is a continuous random function on $\mathbb{R}^n$, and is $\group$-invariant
with probability 1. Consider any finite set of points

\begin{equation*}
x_1,\ldots,x_k\in\mathbb{R}^n\quad\text{ such that }
x_i\not\sim x_j\text{ for all distinct }i,j\leq k\;.
\end{equation*}

Then ${(F(x_1),\ldots,F(x_k))}$ is a Gaussian random vector, with mean and covariance

\begin{equation*}
\mean[F(x_i)]=\mu(\rho(x_i)) \quad\text{ and }\quad
\text{\rm Cov}[F(x_i),F(x_j)]=\kappa(\rho(x_i),\rho(x_j))
\quad\text{ for }i,j\leq k\;.
\end{equation*}

Clearly, $F$ cannot be a Gaussian process on $\mathbb{R}^n$: Since $F$ is invariant,
$F(x)$ completely determines $F(\phi(x))$, so ${(F(x),F(\phi x))}$ cannot be jointly Gaussian.
Put differently, conditioning $F$ on its values on $\Pi$ renders $F$ non-random.
Loosely speaking, the proposition hence says that $F$ is “as Gaussian” as a $\group$-invariant
random function can be. <ref> illustrates random functions generated by such a process.

The construction of MacKay, 1998 described in the introduction
was designed specifically for Gaussian processes, to generate periodic functions at random.
We can now generalize these processes from periodicity to crystallographic invariance:
Given $\group$ and $\Pi$, construct the embedding map ${\rho:\mathbb{R}^n\rightarrow\mathbb{R}^N}$.
Choose $\hat{\kappa}$ as the RBF kernel (<ref>) on $\mathbb{R}^N$, and $\hat{\mu}$ as the constant function $0$ on $\mathbb{R}^N$.
Then generate $F$ as

\begin{equation*}
\qquad\text{ and }\qquad
\end{equation*}

For visualization, draws can be approximated by the randomized feature scheme of
Rahimi and Recht, 2007. <ref> shows examples for $\group$ chosen as
and  on $\mathbb{R}^2$, and for  and  on $\mathbb{R}^3$.

at (0,0) ;
at (3.7,0) ;
at (7.4,0) ;
at (11.1,0) ;

Random invariant functions on $\mathbb{R}^2$ and $\mathbb{R}^3$,
generated by Gaussian processes as described in <ref>.
The groups are, from left to right, , , , and .

§.§ Distributionally invariant processes

Another type of invariance that random functions can satisfy is distributional $\group$-invariance,
which holds if

\begin{equation*}
F\;\equdist\;F\circ\phi\qquad\text{ for all }\phi\in\group\;.
\end{equation*}

Here, $\smash{\equdist}$ denotes equality in distribution.
That is equivalent to requiring that
the distribution $P$ of $F$ satisfies ${P(\phi A)=P(A)}$ for every measurable set $A$.
For crystallographic groups, distributionally invariant Gaussian processes can
be constructed by factoring the parameters, rather than
the random function $F$, through the embedding in <ref>:
Let $\mu$ be a real-valued function and $\kappa$ a kernel on $\mathbb{R}^N$.
If $F$ is the Gaussian process on $\mathbb{R}^n$ with mean ${\mu\circ\rho}$
and covariance function ${\kappa\circ(\rho\otimes\rho)}$, then
$F$ is distributionally $\group$-invariant, i.e. ${F\circ\phi\equdist F}$ for all ${\phi\in\group}$.

Almost sure invariance implies distributional invariance; distributional invariance is
typically a much weaker property.
Frequently encountered examples of distributional invariance are all forms of
stationarity (distributional invariance under shift groups) and of
exchangeability  (permutation groups).

§ THE LAPLACE OPERATOR ON INVARIANT FUNCTIONS

The results in this section describe the behavior of the Laplace operator on
$\group$-invariant functions. All of these are ingredients in the proof of the Fourier representation.
We first describe the transformation behavior of differentials of invariant
functions, in <ref>. Gradients turn out to be invariant under shifts and
equivariant under orthogonal transformations.
Gradient vector fields, and more generally vector fields
with the same transformation behavior as gradients, have a cancellation property—their integral
orthogonal to the tile boundary vanishes (<ref>).
We then define the relevant solution space for the spectral problem, which has Hilbert space
structure (so that we can define orthogonality and self-adjointness) but has smoother elements than
$\L_2$, in <ref>. Once the Laplacian has been properly defined on this space, we can
use the cancellation property to show it is self-adjoint.

§.§ Differentials and gradients of invariant functions

Given a differentiable function ${f:\mathbb{R}^n\rightarrow\mathbb{R}^m}$, denote the differential at $x$ as

\begin{equation*}
\;=\;
\left(\mfrac{\partial f_j}{\partial x_i}\right)_{i\leq n,j\leq m}
\;\in\mathbb{R}^{n\times m}\;.
\end{equation*}

The next result summarizes how invariance of $f$ under a transformation $\phi$ affects $Df$.
Note the order of operations matters: $D(f\circ\phi)$ is the differential of the function $f\circ\phi$, whereas
${(Df)\circ\phi}$ transforms the differential $Df$ of $f$ by $\phi$.

If ${f:\mathbb{R}^n\rightarrow\mathbb{R}^m}$ is
invariant under an isometry ${\phi x=A_\phi x+b_\phi}$ and differentiable, then

\begin{equation}
\label{eq:equivariant:differential}
(Df)(\phi x)\;=\;Df(x)\cdot A_{\phi}^{\trans}\;.
\end{equation}

If in particular ${f:\mathbb{R}^n\rightarrow\mathbb{R}}$, the gradient satisfies

\begin{equation}
\label{eq:equivariant:gradient}
\nabla f(\phi x)\;=\;A_{\phi}\cdot\nabla f(x) \qquad\text{ for all }x\in\mathbb{R}^n\;.
\end{equation}

The Hessian matrix $H_f$ and the Laplacian satisfy

\begin{equation}
\label{eq:invariant:Laplacian}
H_f(\phi x)\;=\;A_{\phi}H_f(x)A_{\phi}^{\trans}
\quad\text{ and }\quad
\Delta f(\phi x)\;=\;\Delta f(x)
\qquad\text{ for all }x\in\mathbb{R}^n\;.
\end{equation}

Since $\phi$ is affine, its differential ${(D\phi)(x)=A_\phi}$ is constant. The chain rule

\begin{equation*}
\;=\;
(Df)\circ\phi\cdot (D\phi)
\;=\;
((Df)\circ\phi)\cdot A_\phi\;.
\end{equation*}

By invariance, ${f\circ\phi}$ and $f$ are the same function, and hence ${D(f\circ\phi)=Df}$.
Substituting into the identity above shows (<ref>), since ${A_\phi^{\trans}=A_\phi^{-1}}$.
For ${m=1}$, the transpose ${D^{\trans}=\nabla}$ is the gradient, and (<ref>) becomes (<ref>). Using (<ref>), the Hessian can be written as

\begin{equation*}
H_f\;=\;D(\nabla f)\;=\;D(A_{\phi}\nabla f\circ\phi^{-1})\;.
\end{equation*}

Another application of the chain rule then shows

\begin{align*}
\;=\;
D(A_{\phi}\nabla f)\circ\phi^{-1}\cdot D\phi^{-1}
\;=\;
A_{\phi}(D\nabla f)\circ\phi^{-1}\cdot A_{\phi^{-1}}
\;=\;
A_{\phi}(H_f\circ \phi^{-1})A_{\phi}^{\trans}\;,
\end{align*}

which is the first statement in (<ref>). Since the Laplacian is the trace of $H_f$, and
the trace in invariant under change of basis, that implies

\begin{equation*}
\Delta f(\phi x)
\;=\;
\text{tr}(H_f(\phi x))
\;=\;
\text{tr}(A_{\phi}H_f(x)A_{\phi}^{\trans})
\;=\;
\text{tr}(H_f(x))
\;=\;
\Delta f(x)\;.\qedhere
\end{equation*}

§.§ Flux through the tile boundary

The next result is the key tool we use to prove self-adjointness of the Laplacian.
We have seen above that the gradient of a $\group$-invariant function transforms
under $\group$ according to (<ref>).
We now abstract from the specific function $\nabla f$,
and consider any vector field ${F:\Pi\rightarrow\mathbb{R}^n}$
that transforms like the gradient on the tile boundary, i.e.

\begin{equation}
\label{periodic:condition:field}
F(y)\;=\;A_{\phi}F(x)\qquad\text{ whenever }y=\phi x\;.
\end{equation}

For a polytope $\Pi$ with facets ${S_1,\ldots,S_k}$, we define the
normal field on the boundary as

\begin{equation*}
\n_{\Pi}:\partial\Pi\rightarrow\mathbb{R}^n
\qquad\text{ given by }\qquad
\n_{\cF}(x)\;:=\;\begin{cases}
\n_i & \text{ if }x\in S_i^{\circ}\\
0 & \text{ otherwise }
\end{cases}
\end{equation*}

where $\n_i$ is the unit normal vector of the facet $S_i$, directed outward with respect to $\Pi$.
In vector analysis, the projection ${F^{\trans}\n_{\Pi}}$ of a vector field onto the direction orthogonal
to $\partial\Pi$ is known as the flux of $F$ through the boundary.

Let $\group$ be a crystallographic group that tiles $\mathbb{R}^n$ with a convex polytope $\Pi$.
If a vector field ${F:\Pi\rightarrow\mathbb{R}^n}$ is integrable on $\partial\Pi$
and satisfies (<ref>), then

\begin{equation*}
\mint_{\partial\cF}F(x)^{\trans}\n_{\cF}(x)\vol_{n-1}(dx)\;=\;0    \;.
\end{equation*}

See <ref>.

§.§ The Sobolev space of invariant functions

The proof of <ref> follows a well-established strategy in spectral theory:
The relevant spectral results hold for self-adjoint operators, and self-adjointness can only
be defined with respect to an inner product. Since the space $\C^2$ on which the Laplace
operator is defined is a Banach space, but has no inner product, one must hence first embed the problem into
a suitable Hilbert space. For the Laplacian, this is generally a first-order Sobolev space;
see <ref> for a review of definitions, and Brezis, 2011, Maz'ya, 2011, McLean, 2000 for more
on spectral theory and the general approach.

In our case, we proceed as follows:
Since invariant functions are completely determined by their values on $\Pi$, we can equivalently
solve the problem on the bounded domain $\Pi$ rather than the unbounded domain $\mathbb{R}^n$.
That gives us access to a number of results specific to bounded domains.
We also observe that the invariance constraint ${e=e\circ\phi}$ is a linear constraint—if two functions
satisfy it, so do their linear combinations—so the feasible set of this constraint is a vector space,
and we can encode the constraint by restriction to a suitable subspace.
We start with the vector space

\begin{equation}
\label{definition:space:H}
\mathcal{H}
\;:=\;
\braces{f|_{\Pi^{\circ}}\,|\,f:\mathbb{R}^n\rightarrow\mathbb{R}\text{ infinitely often differentiable and $\group$-invariant}}\;.
\end{equation}

The elements of $\mathcal{H}$ are hence infinitely often differentiable on $\Pi^{\circ}$, and their
continuous extensions to the closure $\Pi$ satisfy the periodic boundary condition (<ref>).
We then define the Sobolev space of candidate solutions as

\begin{equation*}
\H_{\group}\;:=\;\text{ closure of $\mathcal{H}$ in }\H^1(\Pi^{\circ})\;,
\end{equation*}

equipped with the norm and inner product of ${\H^1(\Pi^\circ)}$.
As a closed subspace of a Hilbert space, it is a Hilbert space.

§.§ The Laplace operator on $\H_{\group}$

We now have to extend $\Delta$ to all elements of $\H_\group$. In general,
a linear operator $\Lambda$ on a closed subspace ${V\subset\H^1(\Pi^\circ)}$ is an
extension of $\Delta$ to $V$ if it satisfies

\begin{equation}
\label{extends}
\Lambda f\;=\;\Delta f\qquad\text{ for all }f\in V\cap\C^2(\Pi^\circ)\;.
\end{equation}

The extended operator is self-adjoint on $V$ if

\begin{equation*}
\sp{\Lambda f,h}_{\H^1}\;=\;\sp{f,\Lambda h}_{\H^1}\qquad\text{ for all }f,h\in V\;.
\end{equation*}

To prove self-adjointness, one decomposes $\Lambda$ as

\begin{equation*}
\sp{-\Delta f,h}_{\L_2}
\,=\,
(\text{integral over $\Gamma^{\circ}$ that is symmetric in $f$ and $h$})
\;-\;
(\text{integral over }\partial\Gamma)\;.
\end{equation*}

This is the Green identity alluded to in the introduction.
To make it precise, we need two quantities: One is the energy form
or energy product

\begin{equation}
\label{energy:form}
a(f,h)\;:=\;\mint_{\Gamma}\nabla f(x)^{\trans}\nabla h(x)\,\vol_n(dx)\;.
\end{equation}

Since it only involves first derivatives, and both appear under the integral, it is
well-defined for any ${f,h\in\H^1(\Pi^{\circ})}$, and is hence a symmetric bilinear form
${a:\H^1\times\H^1\rightarrow\mathbb{R}}$. It is positive definite, since

\begin{equation}
\label{energy:positive:semidefinite}
a(f,h)\;=\;\msum_{i\leq n}\sp{\partial_i f,\partial_i h}_{\L_2}
\quad\text{ and hence }\quad
a(f,f)\;=\;\msum_{i\leq n}\|\partial_if\|_{\L_2}\;\geq\;0\;.
\end{equation}

Substituting the definition of $a$ into that of the $\H^1$ scalar product
in (<ref>) shows that

\begin{equation}
\label{eq:H1:L2:a}
\sp{f,h}_{\H^1}\;=\;\sp{f,h}_{\L_2}\,+\,a(f,h)
\qquad\text{ for all }f,h\in\H^1(\Gamma^\circ)\;.
\end{equation}

The second quantity is the conormal derivative

\begin{equation*}
\partial_{\n}f(x)\;:=\;\nabla f(x)^{\trans}\n_\Pi(x)\;.
\end{equation*}

The precise statement of the decomposition above is then as follows.
[Green's identity]

If the domain $\Pi$ is sufficiently regular—in particular, if $\Pi$ is a convex polytope—then

\begin{equation*}
\sp{-\Lambda f,h}_{\L_2}
\;=\;
\,-\,
\mint_{\partial\Gamma}\partial_{\n}f(x)h(x)\vol_{n-1}(dx)
\qquad\text{ for }f,h\in\H^1(\intF)\;.
\end{equation*}

Informally, this shows that $\Delta$ “behaves self-adjointly” in the interior
of $\Pi$, where derivatives can be computed in all directions around a point.
At points on $\partial\Pi$, the boundary truncates derivatives in some direction,
and that requires a correction term $\partial_{\n}f$.

Let $\group$ be a crystallographic group that tiles $\mathbb{R}^n$ with a convex polytope $\Pi$.
Then $\Delta$ has a unique extension to a linear operator $\Lambda$ on $\H_\group$.
This operator is self-adjoint and continuous on $\H_\group$, and satisfies

\begin{equation}
\label{eq:theorem:laplacian}
\text{\rm(i) }\;\sp{-\Lambda f,h}_{\L_2}\;=\;a(f,h)
\qquad\text{ and }\qquad
\text{\rm(ii) }
\sp{-\Lambda f,f}_{\H^1}\;\geq\;\|f\|^2_{\H^1}-\|f\|^2_{\L_2}
\end{equation}

for all ${f,h\in\H_\group}$.

The proof uses the flux property to show that crystallographic symmetry
makes the boundary term cancel. Since $a$ is symmetric, that makes $\Lambda$ self-adjoint.
In the parlance of elliptic differential equations, (<ref>ii) says
that $\Lambda$ is coercive on $\H_\group$ (see [McLean, 2000]).
See <ref>.

§.§ Linear representations from a nonlinear ansatz

The properties of Laplace operators lead naturally to a class of numerical approximations
known as Galerkin methods (e.g., Braess, 2007). Using the embedding map
$\rho$, we can derive a Galerkin method that can be used to compute the Fourier basis functions
in <ref>—that is, we can use the nonlinear representation approach in the
numerical approximation of the linear representation. The Galerkin method
can be more accurate than the spectral approach in <ref>,
and was used to render Figures <ref>, <ref>
and <ref>.

Galerkin methods posit basis functions
${\chi_1,\ldots,\chi_m}$ and approximate an infinite dimensional function space by the finite-dimensional subspace
In our case, we approximate solutions $e$ of (<ref>) by approximating their
restrictions $e|_\Pi$. We hence need functions ${\chi_i:\Pi\rightarrow\mathbb{R}}$.
W we start with functions ${\tilde{\chi}_i:\mathbb{R}^N\rightarrow\mathbb{R}}$, and set
${\chi_i:=\tilde{\chi}_i\circ\rho}$. We then assume ${e|_{\Pi}}$ of (<ref>) is in the
span, and hence of the form

\begin{equation}
\label{galerkin:ansatz}
e|_{\Pi}\;=\;\msum_{i\leq m}c_i\chi_i\;.
\end{equation}

If $e$ solves the eigenvalue problem (<ref>), $e|_{\Pi}$ satisfies

\begin{align*}
\sp{-\Delta e|_{\Pi},\chi_j}_{\L_2}\;=\;\lambda\sp{e|_{\Pi},\chi_j}_{\L_2}\qquad\text{ for all }j\leq m\;.
\end{align*}

Applying (<ref>) and substituting in (<ref>) shows

\begin{equation*}
\quad\text{ and }\quad
\msum_{i}c_ia(\chi_i,\chi_j)\,=\,\lambda\msum_ic_i\sp{\chi_i,\chi_j}_{\L_2}\,.
\end{equation*}

If we define matrices with entries ${A_{ij}:=a(\chi_i,\chi_j)}$
and ${B_{ij}:=\sp{\chi_i,\chi_j}_{\L_2}}$, that becomes

\begin{equation*}
Ac\;=\;\lambda Bc\qquad\text{ where }c=(c_1,\ldots,c_m)^{\trans}\;.
\end{equation*}

The entries of $A$ and $B$ can be computed with off-the-shelf cubature
methods, and we can then solve for the pair $(\lambda,c)$.

(a) If $\rho$ and the basis functions are implemented with JAX [Bradbury et al., 2018] or a similar
automatic differentiation tool, the gradients in (<ref>) are available, which avoids finite
difference approximation and explicit computation of second derivatives.

(b) Neumann boundary conditions for reflections (see <ref>) can be
enforced using the methods of Golub, 1973.

(c) The basis functions $\tilde{\chi}_i$ can be almost any basis on $\mathbb{R}^N$.
Figures <ref>–<ref> were rendered by placing points ${x_1,x_2,\ldots}$ uniformly on $\Pi$,
and centering radial basis functions at the points ${\rho(x_i)}$ in $\mathbb{R}^N$.

§ RELATED WORK AND ADDITIONAL REFERENCES

In machine learning.
There has been substantial work on group invariance and equivariance in machine learning, with a focus on finite and compact groups.
Most salient has been work on approximate translation invariance and equivariance in convolutional neural networks for images [LeCun et al., 1989, Krizhevsky et al., 2017] and speech [Abdel-Hamid et al., 2014], although this work has not been framed in a group-theoretic way.
To our knowledge the earliest explicit consideration of compact and finite group structure in machine learning was from a Fourier perspective by Kondor, 2008; this was primarily in the context of Hilbert-space formalisms of learning.
The current perspective on compact and finite group equivariance in deep learning arose largely from Cohen and Welling, 2016.
There has been widespread application of machine learning models when group invariance or equivariance is desired, e.g., permutation invariance for sets [Zaheer et al., 2017] and equivariance for neural auction design [Rahme et al., 2021].
In the natural sciences, rotation invariance has been used for astronomy [Dieleman et al., 2015] and $E(3)$ equivariance has proved important for molecular applications [Batzner et al., 2022].
Permutation equivariance of transformer architectures plays a crucial role in large language model [Vaswani et al., 2017].

In crystallography.
Crystallographers have completely described the 17 two-dimensional and 230 three-dimensional
crystallgraphic groups and various tilings they describe,
and tabulated many of their properties [Hahn, 2011].
The emphasis in this work differs somewhat from that in mathematics—in particular,
work in crystallography emphasizes polytopes $\Pi$ that occur in crystal structures
(and which are not necessarily exact in the terminology used in <ref>),
whereas more abstract work in geometry tends to work with Dirichlet domains or other exact tilings.
A long line of work in the context of X-ray crystallography modifies the matrices
that occur in fast Fourier transforms (FFTs) to speed up computation if a crystallographic symmetry
is present in the data. This starts with the work of
Bienenstock and Ewald, 1962 and Ten Eyck, 1973, see also An et al., 1990.
The introduction of Seguel and Burbano, 2003 gives an overview.
This work does not attempt to derive invariant Fourier bases.

In Fourier and PDE analysis.
As we have already explained in some detail, the special case of <ref>
for ${\Pi=[0,1]^n}$ and ${\group=\mathbb{Z}^n}$ yields
the Fourier transform. For this problem, the periodic boundary condition can be
replaced by a Neumann condition, and spectral problems with Neumann conditions
are standard material in textbooks [Brezis, 2011, Lax, 2002].
For shifts that are not axis-parallel, the
periodic boundary condition is known
as a Born-von Karman boundary condition [Ashcroft and Mermin, 1976]. We are not aware of
extensions to crystallographic groups.
An introduction to the PDE techniques used in our proofs
can be found in Brezis, 2011. The conditions imposed there are too restrictive
for our problems, however; a treatment general enough to cover all results we use
is given by McLean, 2000.

In geometry.
<cit.> coined the term orbifold in the
1970s. Commonly cited references include 54, Bonahon and Siebenmann, 1985, Thurston, 1997;
Apanasov, 2000 has a detailed bibliography.
These all focus on general groups, however, for which the theory is much harder than in our case.
The quotient space structure of crystallographic groups was already understood much earlier
by the Göttingen and Moscow schools [Vinberg and Shvartsman, 1993].
A readable introduction to isometry groups and their quotients is given by Bonahon, 2009.
The comprehensive account of Ratcliffe, 2006 is
more demanding, but covers all results needed in our proofs.
Vinberg and Shvartsman, 1993 cover the
geometric aspects of crystallographic groups. <cit.>
explain the geometry of orbifolds heuristically, with many illustrations.

§ SOME OPEN PROBLEMS

Our approach raises a range of further questions well beyond the scope of the present paper,
including in particular those concerning numerical and statistical accuracy. We briefly
discuss some aspects of this problem.

Linear representation.
Suppose we represent a $\group$-invariant continuous function $f$
by evaluating the generalized Fourier basis in <ref>
using the spectral algorithm in <ref>.
The algorithm returns numerical approximations ${\hat{e}_1,\hat{e}_2,\ldots}$ of the
basis functions. We may then expand $f$ as

\begin{equation*}
f\;\approx\;\msum_{i=1}^m c_i\hat{e}_i\;.
\end{equation*}

There are three principal sources of error in this representation:

* The truncation error, since $m$ is finite.

* Any error incurred in computation of the coefficients $c_i$.

* The error incurred by approximating the actual basis functions $\ef_i$ by $\hat{e}_i$.

The truncation error (1) concerns the question how well the vector space ${\text{span}\braces{\ef_1,\ldots,\ef_m}}$
approximates the space $\C_\group$ or $\L_2(\Pi)$. This problem is
studied in approximation theory. Depending on the context, one may choose the first $m$ basis vectors
(a strategy called “linear approximation” in approximation theory), or greedily choose those $m$ basis
vectors that minimize some error measure (“nonlinear approximation”), see DeVore, 1998.
Problem (2) depends on the function $f$, and on how it is represented computationally. If $f$ must
itself be reconstructed from samples, the coefficients are themselves estimators and incur statistical

The error immediately related to our method is (3), and
for the method of <Ref> depends on how well the graph Laplacian used
in <ref> approximates the Laplacian $\Delta$.
This problem has been studied in a number of fields, including machine learning in the context of dimensionality
reduction [Belkin and Niyogi, 2003] and numerical mathematics in the context of
homogenous Helmholtz equations [Harari and Turkel, 1995], and is the subject of a rich literature
[Lafon, 2004, Belkin and Niyogi, 2008, Singer, 2006, Hein et al., 2005, Dunson et al., 2021, García Trillos et al., 2020].
Available results show that, as ${\varepsilon\rightarrow 0}$ in the $\varepsilon$-net, the matrix $L$ converges to $\Delta$,
where the approximation can be measures in different notions of convergence, in particular pointwise and spectral convergence.
The cited results all concern the manifold case. We are not aware of similar results for orbifolds.

For the method of <Ref>, the error depends largely on the choice of basis in $\mathbb{R}^N$ and the accuracy of the numerical integrals, as well as the orbifold map approximation itself (see below).
Error analysis of the Rayleigh-Ritz method has a long history, see, e.g., Weinberger, 1952, Wendroff, 1965, Weinberger, 1974.

Nonlinear representation.
If we define a $\group$-invariant statistical or machine learning model on $\mathbb{R}^n$ by factoring it through an
orbifold, one may ask approximation questions of a more statistical flavor:
Suppose we define a class ${\mathcal{H}=\braces{h_\theta|\theta\in T}}$ of functions ${h_\theta:\mathbb{R}^N\rightarrow\mathbb{R}}$
on the embedding space $\mathbb{R}^N$, with some parameter space $T$.
We then define a class $\mathcal{F}$ of $\group$-invariant functions on $\mathbb{R}^n$ as

\begin{equation*}
\mathcal{F}\;:=\;\braces{f_\theta|\theta\in T}
\quad\text{ where }\quad
\end{equation*}

Depending on the context, we may think of the functions $f_\theta$ e.g., as neural networks or regressors.
The task is then to conduct inference, i.e., to compute a point estimate $\hat{\theta}$ of $\theta$
(say by maximum likelihood estimation or empirical risk minimization), or to compute a posterior
on $T$ in a Bayesian setup.
Since $\mathcal{H}$ and $\mathcal{F}$ share the same parameter space, any such inference task
can be “pushed forward” forward to the embedding space, that is,

\begin{equation*}
\text{inference under }\mathcal{F}\text{ given }x_1,\ldots,x_n
\;=\;
\text{inference under }\mathcal{H}\text{ given }\rho(x_1),\ldots,\rho(x_n)
\end{equation*}

The error can again be separated into components:

* The statistical error associated with fitting $\braces{h_\theta|\theta\in T}$.

* The “forward distortion” introduced by the map ${x_j\mapsto\rho(x_j)}$.

* The “backward distortion” introduced by the map $h_\theta\mapsto h_\theta\circ\rho$.

Problem (1) reduces to the statistical properties of $\mathcal{H}$, and depends on both the
model and the chosen inference method. Problem (2) and (3), however, raise a number of new questions:
The map $\rho$ is, by <ref>, bijective (which means it does not introduce
identifiability problems) and continuous. As the proof of <ref> shows,
it also preserves density properties of certain function spaces, which can be thought of
as a qualitative approximation result. Quantitative results are different matter:
To bound the effect of transformations on
statistical errors typically requires a stronger property than continuity, such as differentiability
or at least a Lipschitz property. In results on manifold learning, the curvature of
$\Omega$ often plays an explicit role. Orbifolds introduce a further challenge,
since smoothness properties fail at the
tips and edges introduced by points with non-trivial stabilizers.
On the other hand, non-differentiabilities of crystallographic orbifolds
have lower-bounded opening angles [Thurston, 1997]—note the tip of the cone in
<ref>, for example, is not a cusp—so it may be possible to mitigate these problems.

§.§ Acknowledgements

The authors would like to thank Elif Ertekin and Eric Toberer for valuable discussions.
RPA is supported in part by NSF grants IIS-2007278 and OAC-2118201. PO is supported by
the Gatsby Charitable Foundation.

[Abdel-Hamid et al., 2014]
O. Abdel-Hamid, A.-R. Mohamed, H. Jiang, L. Deng, G. Penn, and D. Yu.
Convolutional neural networks for speech recognition.
IEEE/ACM Trans. Audio Speech Lang. Process., 22:0
1533–1545, 2014.

[Adams and Fournier, 2003]
R. A. Adams and J. J. F. Fournier.
Sobolev Spaces.
Academic Press, 2003.

[Aliprantis and Border, 2006]
C. D. Aliprantis and K. C. Border.
Infinite-Dimensional Analysis.
Springer, 2006.

[Aliprantis and Burkinshaw, 2006]
C. D. Aliprantis and O. Burkinshaw.
Positive Operators.
Springer, 2006.

[An et al., 1990]
M. An, J. W. Cooley, and R. Tolimieri.
Factorization method for crystallographic Fourier transforms.
Adv. Appl. Math, 11:0 358–371, 1990.

[Apanasov, 2000]
B. Apanasov.
Conformal geometry of discrete groups and manifolds.
de Gruyter, 2000.

[Ashcroft and Mermin, 1976]
N. W. Ashcroft and N. D. Mermin.
Solid State Physics.
Harcourt, 1976.

[Batzner et al., 2022]
S. Batzner, A. Musaelian, L. Sun, M. Geiger, J. P. Mailoa, M. Kornbluth,
N. Molinari, T. E. Smidt, and B. Kozinsky.
E (3)-equivariant graph neural networks for data-efficient and
accurate interatomic potentials.
Nature Communications, 130 (1):0 2453, 2022.

[Becker and Kechris, 1996]
H. Becker and A. S. Kechris.
The descriptive set theory of Polish group actions.
Cambridge Univ. Press, 1996.

[Belkin and Niyogi, 2003]
M. Belkin and P. Niyogi.
Laplacian eigenmaps for dimensionality reduction and data
Neural Comput., 150 (6):0 1373–1396, 2003.

[Belkin and Niyogi, 2008]
M. Belkin and P. Niyogi.
Towards a theoretical foundation for Laplacian-based manifold
J. Comput. System Sci., 740 (8):0 1289–1308,

[Bell, 1954]
D. G. Bell.
Group theory and crystal lattices.
Rev. of Modern Phys., 260 (3):0 311, 1954.

[Bienenstock and Ewald, 1962]
A. Bienenstock and P. P. Ewald.
Symmetry of Fourier space.
Acta Cryst., 15:0 1253–1261, 1962.

[Bloem-Reddy and Teh, 2020]
B. Bloem-Reddy and Y. W. Teh.
Probabilistic symmetries and invariant neural networks.
J. Mach. Learn. Res., 210 (90):0 1–61, 2020.

[Bonahon, 2009]
F. Bonahon.
Low-Dimensional Geometry.
Amer. Math. Soc., 2009.

[Bonahon and Siebenmann, 1985]
F. Bonahon and L. Siebenmann.
Seifert orbifolds and their role as natural crystalline parts of
arbitrary compact irreducible 3-orbifolds.
In London Math. Soc. Lect. Notes, volume 95, pages 18–85.
Cambridge Univ. Press, 1985.

[Bradbury et al., 2018]
J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, C. Leary, D. al Maclaurin,
G. Necula, A. Paszke, J. VanderPlas, S. Wanderman-Milne, and Q. Zhang.
JAX: composable transformations of Python+NumPy programs,
URL <http://github.com/google/jax>.

[Braess, 2007]
D. Braess.
Finite Elements.
Springer, 2007.

[Brezis, 2011]
H. Brezis.
Functional Analysis, Sobolev Spaces, and Partial Differential
Springer, 2011.

[Chung, 1997]
F. Chung.
Spectral graph theory.
Amer. Math. Soc., 1997.

[Cohen and Welling, 2016]
T. Cohen and M. Welling.
Group equivariant convolutional networks.
In Int. Conf. Machine Learning, pages 2990–2999. PMLR, 2016.

[Conway et al., 2016]
J. H. Conway, H. Burgiel, and C. Goodman-Strauss.
The Symmetries of Things.
CRC Press, 2016.

[Cooper et al., 2000]
D. Cooper, C. D. Hodgson, and S. P. Kerckhoff.
Three-dimensional Orbifolds and Cone-Manifolds.
Math. Soc. Japan, 2000.

[DeVore, 1998]
R. A. DeVore.
Nonlinear approximation.
Acta Numer., pages 51–150, 1998.

[Dieleman et al., 2015]
S. Dieleman, K. W. Willett, and J. Dambre.
Rotation-invariant convolutional neural networks for galaxy
morphology prediction.
Monthly Not. Roy. Astronom. Soc., 4500 (2):0
1441–1459, 2015.

[Dresselhaus et al., 2007]
M. S. Dresselhaus, G. Dresselhaus, and A. Jorio.
Group theory: application to the physics of condensed matter.
Springer, 2007.

[Dunson et al., 2021]
D. B. Dunson, H.-T. Wu, and N. Wu.
Spectral convergence of graph Laplacian and heat kernel
reconstruction in $l_\infty$ from random samples.
Appl. Comput. Harmon. Anal., 55:0 282–336, 2021.

[García Trillos et al., 2020]
N. García Trillos, M. Gerlach, M. Hein, and D. Slepčev.
Error estimates for spectral convergence of the graph Laplacian on
random geometric graphs toward the Laplace–Beltrami operator.
Found. Comput. Math., 200 (4):0 827–887,

[Golub, 1973]
G. H. Golub.
Some modified matrix eigenvalue problems.
SIAM review, 150 (2):0 318–334, 1973.

[Gruber, 2007]
P. M. Gruber.
Convex and Discrete Geometry.
Springer, 2007.

[Hahn, 2011]
T. Hahn.
International tables for crystallography, volume A.
Wiley, 2011.

[Harari and Turkel, 1995]
I. Harari and E. Turkel.
Accurate finite difference methods for time-harmonic wave
J. Comput. Phys., 1190 (2):0 252–270, 1995.

[Haussler and Welzl, 1987]
D. Haussler and E. Welzl.
$\varepsilon$-nets and simplex range queries.
Discrete Comput. Geom., 20 (2):0 127–151,

[Hein et al., 2005]
M. Hein, J.-Y. Audibert, and U. v. Luxburg.
From graphs to manifolds–weak and strong pointwise consistency of
graph Laplacians.
In Proc. of COLT, pages 470–485, 2005.

[Hornik et al., 1989]
K. Hornik, M. Stinchcombe, and H. White.
Multilayer feedforward networks are universal approximators.
Neural Networks, 20 (5):0 359–366, 1989.

[Johnston, 1960]
D. Johnston.
Group theory in solid state physics.
Rep. Progr. Phys., 230 (1):0 66, 1960.

[Killingbeck, 1970]
J. Killingbeck.
Group theory and topology in solid state physics.
Rep. Progr. Physics, 330 (2):0 533, 1970.

[Kondor, 2008]
I. R. Kondor.
Group Theoretical Methods in Machine Learning.
PhD thesis, Columbia Univ., 2008.

[Krizhevsky et al., 2017]
A. Krizhevsky, I. Sutskever, and G. E. Hinton.
Imagenet classification with deep convolutional neural networks.
Communications of the ACM, 600 (6):0 84–90,

[Kruskal, 1964]
J. B. Kruskal.
Multidimensional scaling by optimizing goodness of fit to a nonmetric
Psychometrika, 290 (1):0 1–27, 1964.

[Lafon, 2004]
S. Lafon.
Diffusion maps and geometric harmonics.
PhD thesis, Yale University, 2004.

[Lax, 2001]
M. Lax.
Symmetry principles in solid state and molecular physics.
Courier Corporation, 2001.

[Lax, 2002]
P. D. Lax.
Functional Analysis.
Wiley, 2002.

[LeCun et al., 1989]
Y. LeCun, B. Boser, J. Denker, D. Henderson, R. Howard, W. Hubbard, and
L. Jackel.
Handwritten digit recognition with a back-propagation network.
Adv. Neural Inf. Process. Systems, 2, 1989.

[MacKay, 1998]
D. J. C. MacKay.
Introduction to Gaussian processes.
NATO ASI Series F, 168:0 133–166, 1998.

[Marcus and Rosen, 2006]
M. B. Marcus and J. Rosen.
Markov processes, Gaussian processes, and local times.
Cambridge Univ. Press, 2006.

[Maz'ya, 2011]
V. Maz'ya.
Sobolev Spaces.
Springer, 2nd edition, 2011.

[McLean, 2000]
W. McLean.
Strongly Elliptic Systems and Boundary Integral Equations.
Cambridge Univ. Press, 2000.

[Munkres, 2000]
J. R. Munkres.
Prentice Hall, 2nd edition, 2000.

[Pears, 1975]
A. R. Pears.
Dimension theory of general spaces.
Cambridge University Press, 1975.

[Rahimi and Recht, 2007]
A. Rahimi and B. Recht.
Random features for large-scale kernel machines.
Adv. Neural Inf. Process. Systems, 20, 2007.

[Rahme et al., 2021]
J. Rahme, S. Jelassi, J. Bruna, and S. M. Weinberg.
A permutation-equivariant neural network architecture for auction
In Proc. AAAI, volume 35, pages 5664–5672, 2021.

[Ratcliffe, 2006]
J. G. Ratcliffe.
Foundations of Hyperbolic Manifolds.
Springer, 2nd edition, 2006.

[54]
P. Scott.
The geometries of $3$-manifolds.
Bull. London Math. Soc, 5:0 401–487, 83.

[Seguel and Burbano, 2003]
J. Seguel and D. Burbano.
A scalable crystallographic FFT.
In Lecture Notes in Comput. Sci., volume 2840, pages 134–141.
Springer, 2003.

[Singer, 2006]
A. Singer.
From graph to manifold Laplacian: The convergence rate.
Appl. Comput. Harmon. Anal., 210 (1):0
128–134, 2006.

[Slater, 1965]
J. Slater.
Space groups and wave-function symmetry in crystals.
Rev. Modern Phys., 370 (1):0 68, 1965.

[Speiser, 1956]
A. Speiser.
Die Theorie der Gruppen endlicher Ordnung.
Birkhäuser, 1956.

[Steinwart and Christmann, 2008]
I. Steinwart and A. Christmann.
Support Vector Machines.
Springer, 2008.

[Ten Eyck, 1973]
L. F. Ten Eyck.
Crystallographic fast Fourier transforms.
Acta Cryst. A, 29:0 183–191, 1973.

[Tenenbaum et al., 2000]
J. B. Tenenbaum, V. d. Silva, and J. C. Langford.
A global geometric framework for nonlinear dimensionality reduction.
Science, 2900 (5500):0 2319–2323, 2000.

[Thurston, 1997]
W. P. Thurston.
Three-Dimensional Geometry and Topology, Vol. 1.
Princeton Univ. Press, 1997.

[van der Vaart and van Zanten, 2008]
A. van der Vaart and J. H. van Zanten.
Reproducing Hilbert spaces of Gaussian priors.
In Inst. Math. Stat. (IMS) Collect., volume 3, pages 200–222.

[Vaswani et al., 2017]
A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez,
L. Kaiser, and I. Polosukhin.
Attention is all you need.
Adv. Neural Inf. Process. Systems, 30, 2017.

[Vinberg and Shvartsman, 1993]
E. B. Vinberg and O. V. Shvartsman.
Discrete groups of motions of spaces of constant curvature.
In E. B. Vinberg, editor, Geometry II. Spaces of Constant
Curvature. Springer, 1993.

[Weinberger, 1952]
H. Weinberger.
Error estimation in the Weinstein method for eigenvalues.
Proc. Amer. Math. Soc., 30 (4):0 643–646,

[Weinberger, 1974]
H. F. Weinberger.
Variational methods for eigenvalue approximation.
SIAM, 1974.

[Wendroff, 1965]
B. Wendroff.
Bounds for eigenvalues of some differential operators by the
Rayleigh-Ritz method.
Math. Comp., 190 (90):0 218–224, 1965.

[Zaheer et al., 2017]
M. Zaheer, S. Kottur, S. Ravanbakhsh, B. Poczos, R. R. Salakhutdinov, and A. J.
Deep sets.
Adv. Neural Inf. Process. Systems, 30, 2017.

[Ziegler, 1995]
G. M. Ziegler.
Lectures on Polytopes.
Springer, 1995.

Department of Computer Science

Princeton University

<https://www.cs.princeton.edu/ rpa>

Gatsby Computational Neuroscience Unit

University College London

<https://www.gatsby.ucl.ac.uk/ porbanz>

The first three sections of this appendix provide mathematical background on isometries
(<ref>), function spaces and smoothness (<ref>), orbifolds
(<ref>), and spectral theory (<ref>). The proof of the Fourier representation is subdivided into three parts:
We first prove of the flux property, <ref>, in <ref>, and
<ref> on self-adjointness of the Laplacian in <ref>.
Using these results, we then prove the Fourier representation in <ref>.
The proof of the embedding theorem (<ref>) follows in <ref>.
<ref> collects all proofs on kernels and Gaussian processes.

§ BACKGROUND I: ISOMETRIES OF EUCLIDEAN SPACE

Isometries are invertible functions that preserve distance. To define
an isometry between two sets $V$ and $W$, both must be equipped with metrics,
say $d_V$ and $d_W$.
A map ${\phi:\xspace\rightarrow\yspace}$ is then an isometry if it is one-to-one
and satisfies

\begin{equation*}
\qquad\text{ for all }v_1,v_2\in V\;.
\end{equation*}

Since this implies $\phi$ is Lipschitz, isometries are always continuous.
If ${W=V}$, then $\phi$ is necessarily bijective.
An isometry of $\mathbb{R}^n$ is a bijection ${\phi:\mathbb{R}^n\rightarrow\mathbb{R}^n}$ that satisfies

\begin{equation*}
d_n(\phi x,\phi y)\;=\;d_n(x,y)\qquad\text{ for all }x,y\in\mathbb{R}^n\;.
\end{equation*}

Identity (<ref>) shows that every isometry can be uniquely
represented as an orthogonal transformation followed by a shift.
Loosely speaking, an isometry may shift, rotate, or flip $M$, but cannot change its shape or volume.
Recall that a set $\group$ of functions ${\mathbb{R}^n\rightarrow\mathbb{R}^n}$
is a group if
it contains the identity map $\Id$, and if ${\phi,\psi\in\group}$ implies ${\phi\circ\psi\in\group}$
and ${\phi^{-1}\in\group}$.
The set of all isometries of $\mathbb{R}^n$ forms a group, called the
Euclidean group of order $n$.

§.§ More on crystallographic groups

Representation by shifts and orthogonal transformations.
Since every isometry can be decomposed into an orthogonal transformation
and a shift according to (<ref>), every crystallographic group $\group$ has two natural subgroups:
One is the group

\begin{equation*}
\group_o\;:=\;\braces{\phi\in\group\,|\,\phi(x)=Ax\text{ for some }A\in\mathbb{O}_n}\;=\;\group\cap\mathbb{O}_n
\end{equation*}

of purely orthogonal transformations. This is an example of a point group,
since all its elements have a common fixed point (namely the origin).
It is always finite: Fix any $x$ on the unit sphere in $\mathbb{R}^n$. Then ${\phi(x)}$ is also on the sphere for every ${\phi\in\group_o}$, since ${A_\phi}$ is orthogonal. However, discreteness requires there can only be finitely many such points $\phi(x)$ on the sphere.
The other is the group of pure shifts,

\begin{equation*}
\group_t\;:=\;\braces{\phi\in\group\,|\,\phi(x)=x+b\text{ for some }b\in\mathbb{R}^n}\;.
\end{equation*}

One can show there are
linearly independent vectors ${b_1,\ldots,b_n}$ such that

\begin{equation*}
\group_t\;:=\;\braces{x\mapsto x+b\,|\,b=a_1b_1+\ldots+a_nb_n\text{ for }a_1,\ldots,a_n\in\mathbb{Z}}\;.
\end{equation*}

Thus, the generating set for a crystallographic group on $\mathbb{R}^n$ always includes $n$ linearly independent shifts.

§.§ Equivalence to definitions in the literature

Our definition of a crystallographic group in <ref> differs
from those in the literature—we have chosen it for simplicity, but must verify it is equivalent.
There are two standard definitions of crystallographic groups: Perhaps the most common one, used for example by
Thurston, 1997, is as a discrete group of isometries for which ${\mathbb{R}^n/\group}$ is compact in the quotient topology.
Another is as a group of isometries of $\mathbb{R}^n$ such that ${\mathbb{R}^n/\group}$ has finite volume when identified with a subset of Euclidean space
[Vinberg and Shvartsman, 1993]. These are known to be equivalent <cit.>.
Our definition is equivalent to both:

A group $\group$ is crystallographic in the sense of <ref>
if and only if it is a discrete group of isometries of $\mathbb{R}^n$ such
that ${\mathbb{R}^n/\group}$ is compact.

If $\group$ is crystallographic in our sense,
it is discrete (see <ref>), and $\mathbb{R}^n/\group$ is compact
by <ref>, so it satisfies the second definition above.
Conversely, if $\group$ satisfies Thurston's definition,
it tiles $\mathbb{R}^n$ with some set $\Pi$. This set can always be chosen as a convex polytope
<cit.>, so $\group$ is crystallographic in our sense.

We note only en passe that there are tilings that cannot be described
by a group of isometries. That is not at all obvious—the question was one of Hilbert's
problems—but counter-examples of such tilings (with non-convex polytopes) are now
known <cit.>.

§ BACKGROUND II: FUNCTION SPACES

This section briefly reviews concepts from functional analysis that play a role in the proofs.
Helpful references include Aliprantis and Border, 2006, Brezis, 2011 on general functional analysis
and Banach spaces, Brezis, 2011, Adams and Fournier, 2003 on Sobolev spaces,
Steinwart and Christmann, 2008 on reproducing kernel Hilbert spaces, and Aliprantis and Burkinshaw, 2006
on compact operators.

§.§ Spans and their closures

Consider a Banach space $V$ and a subset ${\mathcal{F}\subset V}$.
The span of $\mathcal{F}$ is the set

\begin{equation*}
\text{span}(\mathcal{F})
\braces{\tsum_{i\leq n}c_if_i\,|\,n\in\mathbb{N},c_i\in\mathbb{R},f_i\in\mathcal{F}}
\end{equation*}

of finite linear combinations of elements of $\mathcal{F}$.
Since function spaces are typically infinite-dimensional, we also consider
infinite linear combinations. These are defined with respect to a norm $\|\argdot\|$:

\begin{equation*}
\quad\text{ means }\quad
\|f-\tsum_{i\leq n}c_if_i\|\;\rightarrow\;0\;\text{ as }n\rightarrow\infty\;.
\end{equation*}

In other words, to get from the span to the set of infinite linear combinations,
we take the closure in the relevant norm:

\begin{equation*}
\braces{\tsum_{i\in\mathbb{N}}c_if_i\,|\,c_i\in\mathbb{R},f_i\in\mathcal{F}}
\;=\;
\overline{\text{span}(\mathcal{F})}
\end{equation*}

§.§ Bases

A Hilbert space $\mathbb{H}$ is a Banach space whose norm is induced by an inner product $\sp{\argdot,\argdot}_\mathbb{H}$,
that is,

\begin{equation*}
\|f\|_{\mathbb{H}}\;=\;\sqrt{\vphantom{\sp{\argdot,\argdot}}\smash{\sp{f,f}_\mathbb{H}}}\;.
\end{equation*}

A sequence ${f_1,f_2,\ldots}$ in a Hilbert space is an
orthonormal system if ${\sp{f_i,f_j}=\delta_{ij}}$, where $\delta$ is the Kronecker symbol
(the indicator function of $\braces{i=j}$).
An orthonormal system is complete if its span is dense in $\mathbb{H}$, that is, if

\begin{equation*}
\mathbb{H}\;=\;\overline{\text{span}\braces{f_1,f_2,\ldots}}\;,
\end{equation*}

where the closure is taken in the norm of $\mathbb{H}$. A complete orthonormal system
is also called an orthonormal basis. If ${f_1,f_2,\ldots}$ is an orthonormal
basis, $\mathbb{H}$ can be represented as

\begin{equation*}
\mathbb{H}
\;=\;
\bigl\lbrace
\tsum_{i\in\mathbb{N}}c_if_i
\text{ (convergence in $\|\argdot\|_{\mathbb{H}}$)}\,|\,
c_1,c_2,\ldots\in\mathbb{R}\text{ with }\tsum_ic_i^2<\infty
\bigr\rbrace\;.
\end{equation*}

§.§ $\L_2$ spaces

For any $M$ and a
$\sigma$-finite measure $\nu$ on $M$, the $\L_2$-scalar product and pseudonorm are

\begin{equation*}
\sp{f,g}_{\L_2(\nu)}\;:=\;\mint_{M}f(x)g(x)\nu(dx)
\quad\text{ and }\quad
\|f\|_{\L_2(\nu)}:=\sqrt{\sp{f,f}_{\L_2(\nu)}}\;.
\end{equation*}

To make $\|\argdot\|_{\L_2}$ a norm, one defines the equivalence classes
${[f]:=\braces{g\,|\,\|f-g\|_{\L_2}\!=0}}$ of functions identical outside a null set, and
the vector space

\begin{equation*}
\L_2(\nu)\;:=\;\braces{[f]\,|\,f:M\rightarrow\mathbb{R}\text{ and }\|f\|_{\L_2}<\infty}
\end{equation*}

of such equivalence classes, which is a separable Hilbert space.
Although its elements are not technically functions,
we use the notation ${f\in\L_2}$ rather than ${[f]\in\L_2}$.
We write ${\L_2(\mathbb{R}^n)}$ and ${\L_2(\Pi)}$ respectively if $\nu$ is Euclidean volume
on $\mathbb{R}^n$ or on $\Pi$.
See Aliprantis and Border, 2006 or Brezis, 2011 for background on $\L_2$ spaces.

§.§ Reproducing kernel Hilbert spaces

Consider a set ${M\subseteq\mathbb{R}^n}$.
A symmetric positive definite function ${\kappa:M\times M\rightarrow\mathbb{R}}$
is called a kernel.
A kernel defines a Hilbert space as follows: The formula

\begin{equation*}
\sp{\tsum_{i}a_i\kappa(x_i,\argdot),\,\tsum_{j}b_j\kappa(y_j,\argdot)}_{\mathbb{H}}
\;:=\;
\tsum_{i,j}a_ib_i\kappa(x_i,y_j)\quad\text{ for }a_i,b_j\in\mathbb{R}\text{ and }x_i,y_j\in M
\end{equation*}

defines a scalar product on ${\text{span}\braces{\kappa(x,\argdot)|x\in M}}$.
The closure

\begin{equation*}
\mathbb{H}\;:=\;\overline{\text{span}\braces{\kappa(x,\argdot)|x\in M}}
\quad\text{ with respect to the norm }\quad
\|f\|_{\kappa}\;:=\;\sqrt{\sp{f,f}_{\mathbb{H}}}
\end{equation*}

is a real, separable Hilbert space with inner product ${\sp{\argdot,\argdot}_{\mathbb{H}}}$, called the
reproducing kernel Hilbert space or RKHS of $k$.
Every RKHS satisfies the “reproducing property”

\begin{equation}
\label{eq:reproducing}
f(x)\;=\;\sp{f,\kappa(x,\argdot)}_{\mathbb{H}}\qquad\text{ for all }f\in\mathbb{H}\text{ and all }x\in M\;.
\end{equation}

In particular, ${\kappa(x,y)=\sp{\kappa(x,\argdot),\kappa(y,\argdot)}_{\mathbb{H}}}$.
If ${f_1,f_2,\ldots}$ is an orthonormal basis of $\mathbb{H}$, then

\begin{equation}
\label{RKHS:ONS}
\kappa(x,y)\;=\;\tsum_{i\in\mathbb{N}}f_i(x)f_i(y)\qquad\text{ for all }x,y\in M\;.
\end{equation}

If $\mathbb{H}$ is an RKHS, the map ${f\mapsto f(x)}$ is continuous for each ${x\in M}$.
Conversely, if $\mathbb{H}$ is any Hilbert space of real-valued functions on $M$, and if
the maps are continuous on $\mathbb{H}$ for all ${x\in M}$, there
is a unique kernel satisfying (<ref>) that generates
$\mathbb{H}$ as its RKHS.

§.§ Spaces of continuous functions

For any set $M$, the vector space $\C(M)$ of continuous functions equipped with the
norm ${\|\argdot\|_{\sup}}$ is a Banach space. It is separable if $M$ is compact
[Aliprantis and Border, 2006].
In the proof of the spectral theorem, we must also consider the set

\begin{equation*}
\C_u(M)\;:=\;\braces{f\in\C(M)\,|\,f\text{ uniformly continuous}}\;,
\end{equation*}

and the compactly supported functions

\begin{equation*}
\C_c(M)\;=\;\braces{f\in\C(M)\,|\,f=0\text{ outside a compact set }K\subset M}\;.
\end{equation*}

We recall some basic facts from analysis that are used in the proofs:
[Aliprantis and Border, 2006]

(i) Every continuous function on a compact set is uniformly continuous.
(ii) Every uniformly continuous function $f$ on a set ${M\subseteq\mathbb{R}^n}$ has a unique continuous
extension $\bar{f}$ to the closure $\overline{M}$. Its value at a boundary point ${x\in\partial M}$ is given by
${\bar{f}(x)=\lim_j(x_j)}$ for any sequence of points ${x_i\in M}$ with ${x_i\rightarrow x}$.

§.§ Smoothness spaces

Smoothness spaces quantify the smoothness of functions in terms of a norm.
Two types of such spaces play a role in our results, namely $\C^k$ spaces
and Sobolev spaces. Both define smoothness via derivatives:
We denote partial derivatives as

\begin{equation*}
\partial^{\alpha}f
\;:=\;
\frac{\partial^{|\alpha|}f}
{\partial x_1^{\alpha_1}\cdots\partial x_n^{\alpha_n}}
\qquad\text{ where }\alpha=(\alpha_1,\ldots,\alpha_n)\in\mathbb{N}^n\text{ and }|\alpha|:=\alpha_1+\ldots+\alpha_n\;.
\end{equation*}

If we are taking a derivative with respect to the $i$th coordinate, we use a subscript,

\begin{equation*}
\partial_i f\;:=\;\mfrac{\partial f}{\partial x_i}
\end{equation*}

The set $\C^k$ of $k$ times continuously differentiable functions can then be
represented as

\begin{equation*}
\mathbf{C}^k(M)
\;:=\;
\braces{f\in\mathbf{C}(M)\,|\,\partial^{\alpha}f\in\mathbf{C}(M)\text{
whenever }
|\alpha|\leq k}
\qquad\text{ where }k\in\mathbb{N}\cup\braces{0,\infty}\;.
\end{equation*}

Since that means the norm of $\mathbf{C}$ is applicable to $\partial^{\alpha}f$,
we can define

\begin{equation*}
\|f\|_{\mathbf{C}^k}
\;:=\;
\|f\|_{\sup}
\msum_{|\alpha|\leq r} \|\partial^\alpha f\|_{\sup}\;.
\end{equation*}

It can be shown that this is again a norm, and that it makes $\C^k$ a Banach space
[Brezis, 2011].
$\C^k$ functions are uniformly continuous, and even very smooth functions approximate
elements of $\L_2$ to arbitrary precision:

Let ${M\subseteq\mathbb{R}^n}$ be a set.
(i) If ${f\in\C^k(M)}$ for ${k\geq 1}$, then $f$ and its first ${k-1}$ derivatives are uniformly continuous.
(ii) The set ${\C_c(M)\cap\C^{\infty}(M)}$ is dense in $\L_2(M)$.

The $\C^k$ norms measure smoothness in a worst-case sense. To measure average smoothness instead,
we can replace the sup norm by the $\L_2(M)$-norm: The function

\begin{equation*}
\|f\|_{\H^k}
\;:=\;
\|f\|_{\L_2}
\msum_{|\alpha|\leq k} \|\partial^\alpha f\|_{\L_2}\;,
\end{equation*}

is a norm, called the Sobolev norm of order $k$. It makes the set

\begin{equation*}
\H^k(M):=\braces{f\in\L_2(M)\,|\,\|f\|_{\H^k}<\infty}\;=\;
\braces{f\in\L_2(M)\,|\,\partial^\alpha f\in\L_2(M)}
\end{equation*}

a Banach space, and even a Hilbert space, called the Sobolev space of order $k$.
We will only work with the spaces $\H^1(M)$. A inner product on $\H^1(M)$ is given by

\begin{equation*}
\label{eq:sobolev:sp}
\sp{f,g}_{\H_1}\;:=\;\sp{f,g}_{\L_2}\;+\;\msum_{i\leq n}\sp{\partial_i f,\partial_i g}_{\L_2}\;.
\end{equation*}

The Sobolev norms are stronger than the ${\L_2}$ norm: We have

\begin{equation*}
\|f\|_{\L_2(M)}\;\leq\;c_M\|f\|_{\H_1(M)}\qquad\text{ for all }f\in\L_2(M)\text{ and some }c_M>0\;.
\end{equation*}

Consequently, the approximation property in <ref>(ii) does not
necessarily hold in the Sobolev norm. Whether it does depends on whether the geometry of
the domain $M$ is sufficiently regular:

Let $\Gamma$ be a Lipschitz domain (such as a convex polytope). Then