Predicting PDEs Fast and Efficiently with Equivariant Extreme Learning Machines
Extreme Learning Machines for PDEs

Hans Harder Sebastian Peitz
Paderborn University, Paderborn, Germany

{hans.harder<EMAIL_ADDRESS>
We utilize extreme learning machines for the prediction of partial differential equations (PDEs). Our method splits the state space into multiple windows that are predicted individually using a single model. Despite requiring only few data points (in some cases, our method can learn from a single full-state snapshot), it still achieves high accuracy and can predict the flow of PDEs over long time horizons. Moreover, we show how additional symmetries can be exploited to increase sample efficiency and to enforce equivariance.

§ INTRODUCTION

The efficient modeling, prediction, and control of complex systems governed by partial differential equations (PDEs) via data-driven and machine-learning-based techniques has received massive attention in recent years. PDEs are an extremely important class of equations that govern the physics of many applications such as heat transfer, continuum mechanics, fluid dynamics, or quantum mechanics. Classical approaches such as finite differences or finite elements transform the problem into a (potentially very high-dimensional) set of coupled ordinary differential equations. The cost of these solution methods increases drastically with the complexity of the dynamics, and it is the goal of data-driven methods to automatically identify underlying low-dimensional structures and thereby mitigate the curse of dimensionality [9].

In the past decades, we have seen a plethora of various surrogate-modeling techniques for PDEs. Popular examples are the Proper Orthogonal Decomposition [33, 12, 26]—where we project the known equations onto a low-dimensional subspace spanned by basis functions learned from data—as well as many non-intrusive methods such as the Dynamic Mode Decomposition [32], i.e., the Koopman operator framework [31, 19, 17, 24, 25].
Quite unsurprisingly, deep learning plays a central role as well, see, e.g., [41, 21, 43, 39].

More recently, the notion of physics-informed machine learning [11] has become increasingly popular. There—instead of treating the system as a black box—one encodes system knowledge into the training process, for instance in the form of the governing equations [7, 14, 38, 16, 29]. This significantly reduces the amount of required training data, while at the same improving robustness.
Alternatively, one may also enforce known symmetries—also referred to as geometric deep learning [1] in the machine learning community—to increase the efficiency, see, for instance, [3, 6] for translational and rotational symmetries in image classification. In reinforcement learning control of PDEs, the exploitation of translational symmetries gives rise to very efficient multi-agent concepts [42, 8, 35, 27].

When it comes to symmetry exploitation in data-driven modeling of PDEs, only translational equivariance has been exploited until now [23, 28, 22, 36]. These works are based on reservoir computing, i.e., recurrent architectures, and it is well-known that these can be challenging to train and require a substantial initialization phase to adjust the latent variables before they can make useful predictions.
The goal of this paper is to increase the sample efficiency while avoiding recurrence. To this end, we are going to make use of Extreme Learning Machines (ELMs) [37, 4, 10, 5] to learn the flow map of chaotic PDEs. ELMs are two-layer networks that can be trained using simple linear regression.
We demonstrate that when exploiting translational equivariance—similarly to [23]—ELMs can be trained from few time steps and make accurate predictions over long time horizons (<Ref>). In some instances, a single time step suffices to learn a precise approximation of the PDE (<Ref>). Including additional symmetries (i.e., rotations) further increases the efficiency (also <Ref>).

All simulations have been performed using the spectral PDE solver  [18].
Our  code for the experiments is freely available under <https://github.com/graps1/ELM-PDEs>.

§ PRELIMINARIES

§.§ Partial Differential Equations

Consider a function space (e.g., Sobolev space) $\Fc$ containing state functions of the form $v : \R^n \rightarrow \R$ for $v \in \Fc$, where $n$ is the number of spatial dimensions.
We are interested in PDEs where the temporal derivative is given by

\begin{align*}
\partial_t u(\cdot, t) = f(u(\cdot, t)),
\end{align*}

where $u(\cdot, t) \in \Fc$ is the state at time $t \geq 0$ and $f: \Fc \rightarrow \Fc$ is such that $f(v)(x)$ for $v \in \Fc, x \in \R^n$ depends only on the coordinate $x$ and the (spatial) derivatives of $v$ evaluated at $x$.
We assume that there is a unique flow

\Phi : \Fc \times [0,\infty) \rightarrow \Fc, (v,t) \mapsto \Phi_t(v)

associated with the PDE so that $f(\Phi_t(v)) = \partial_t \Phi_t(v)$.
In other words, we must have $u(\cdot, t) = \Phi_t(u(\cdot, 0))$ for any $u(0) \in \Fc$ and $t \geq 0$.
Since our method serves as an approximation to the flow and not the temporal derivative, the time step $\Delta t$ will be fixed depending on the context, allowing us to drop the $\Delta t$ from $\Phi_{\Delta t}(v) =: \Phi(v)$.

In this regard, we will study three variants of the Kuramoto–Sivashinsky equation under periodic boundary conditions. The first one is used for the case of two spatial dimensions, defined by

\begin{equation}
\partial_t u(x,t) =
-\Laplace u(x,t) - \Laplace^2 u(x,t)
- \frac 1 2 | \nabla u(x,t) |^2,
\label{eq:KS_Nd}
\end{equation}

where $\Laplace$, $\Laplace^2$ and $\nabla$ are the Laplace, biharmonic and gradient operators. Another variant is given for one spatial dimension by

\begin{equation}
\partial_t u(x,t) =
-\partial_{xx}^2 u(x,t) - \partial_{xxxx}^4 u(x,t)
- \partial_x u(x,t) u(x,t),
\label{eq:KS_1d}
\end{equation}

which is obtained by taking (<ref>) for the case of one spatial dimension, differentiating with respect to $x$ and then substituting $\partial_x u$ by $u$. This form is popular in fluid dynamics, and, due to its simplicity but still chaotic behavior, often used as a test bed for learning algorithms [23, 28, 36]. Finally, we consider a variant of (<ref>) studied by [23] that introduces an explicit dependence on the spatial coordinate:

\begin{equation}
\partial_t u(x,t) =
-\partial_{xx}^2 u(x,t) - \partial_{xxxx}^4 u(x,t)
- \partial_x u(x,t) u(x,t) + \mu \cos\big(2 \pi x / \lambda\big),
\label{eq:KS_1d_spatial}
\end{equation}

where $\mu$ and $\lambda$ are real-valued parameters.

There are multiple ways to handle the spatial dimensions of a PDE, but one common approach is to discretize the space into a grid in which every node stores an estimate of the PDE's state at that coordinate.
For a one-dimensional system like (<ref>), the state function is approximated by a vector, in the case of two dimensions by a matrix, etc.
To obtain spatial derivatives, one approach it to use finite-difference coefficients (or stencils) that linearly aggregate the values of neighboring nodes.
In our setting,
this means that the value of the temporal derivative at a node's position only depends on a local neighborhood and spatial coordinate.
(When it comes to learning these PDEs, this is commonly exploited using CNN-based architectures [29, 39, 15].)

To give an example, (<ref>) can be defined on $x \in (0,L)$ with periodic boundary conditions, so we may assume that its domain is all of $\R$. Discretizing $v \in \Fc$ near a center $\overline x \in \R$ gives a $k$-dimensional slice

\begin{align*}
\sl(v; \overline x, k) = (v(\overline x_{-(k-1)/2}),v(\overline x_{-(k-1)/2+1}),\dots, v(\overline x_{+(k-1)/2}))^T \in \R^k,
\end{align*}

where $\overline x_i = \overline x + \Delta x i$, with $\Delta x > 0$ being a step size.
Finite-difference coefficients [13, 34] are $k$-dimensional vectors[We assume here that there is one universal $k$ for all finite-difference vectors. This is justified since we use only finitely many.] used to estimate spatial derivatives by taking inner products with $\sl(v; \overline x, k)$.
That is, for $i \geq 1$, there is $c_i \in \R^{k}$ so that

\partial_x v(\overline x) \approx c_1^T  \sl(v; \overline x, k),\quad
\partial^2_{xx} v(\overline x) \approx c_2^T  \sl(v; \overline x, k), \dots

Hence, for the system (<ref>),

\begin{align*}
f(v)(\overline x) \approx -c_2^T \sl(v; \overline x, k) -c_4^T \sl(v; \overline x, k)
- c_1^T \sl(v; \overline x, k) v(\overline x) + \mu \cos(2 \pi \overline x / \lambda ).
\end{align*}

Notice that the right-hand side
describes a function that depends only on $v(\overline x)$'s immediate neighbors $v(\overline x - \Delta x (k-1)/2), \dots, v(\overline x + \Delta x (k-1)/2)$
and $\overline x$.

The simplest way to achieve discretization in time is to use the explicit Euler method, defined by the approximation $\Phi_{\Delta t}(v) \approx v + \Delta t f(v)$.
Combining this with the spatial discretization, notice that $\Phi_{\Delta t}(v)(\overline x)$ is also approximated by a function of $v(\overline x)$'s neighbors and $\overline x$.

This justifies our approach, which is similar to that of [23]; we learn a mapping that approximates $\Phi_{\Delta t}(v)(\overline x)$ from $\overline x$ and $\sl(v; \overline x, k)$.
For more than one spatial dimension, this approach can be generalized straightforwardly: In two dimensions, the discretization is a matrix, and the value at a certain location is approximated by a function accepting a two-dimensional window as its input.

§.§ Extreme Learning Machines

Extreme learning machines (ELMs) [37, 4, 10] are two-layered neural networks. The first layer's weights are chosen randomly, and only the second layer is trained with its bias set to zero. In contrast to the usual gradient-based methods, this is done by directly solving a least-squares problem. Given a pair of random variables $(\xf,\yf)$ modeling the input-output relationship with values in $\R^N$ and $\R^n$ respectively, extreme learning machines solve the problem

\begin{align*}
\min_{\Theta \in \R^{n \times H}} \E[|| \Theta \phi(\xf)  - \yf ||^2_2],
\end{align*}

where $H$ is the hidden dimension and $\phi : \R^N \rightarrow \R^{H}$ composes the first layer with the activation function. Therefore, only the linear output layer, parameterized by $\Theta$, is learned. Setting the derivative of the loss with respect to $\Theta$ zero, we may express the optimal solution as

\begin{equation}
\Theta^* = \E[\yf \phi(\xf)^T] \E[\phi(\xf) \phi(\xf)^T]^{-1} = CD^{-1}.
\label{eq:optimal_Theta}
\end{equation}

The $n\times H$ and $H \times H$-matrices $C$ and $D$ are estimated by averaging over known input-output pairs. Technically, it is not necessary to frame $C$ and $D$ as expectations, but we found it useful since it allows us to flexibly stop the training process when the estimates stop changing.

The mapping $\phi$ can be seen as a random $H$-dimensional dictionary of observables. That is, it takes the data and maps it into a potentially higher dimension, making it easier to learn a linear mapping from embedded input to output.
This training procedure is shared by, or at least similar to, many other approaches such as sparse-identification of nonlinear dynamics (SINDy) [2], radial-basis function networks [20] or extended dynamic mode decomposition [40]. Due to the random initialization, an advantage of ELMs is that the hyperparameter tuning is omitted up to choosing the activation function and an initialization strategy for the first layer. Its advantage over gradient-based methods is that the optimal solution is recovered instantaneously and with high precision.

We will show how one can construct ELMs that approximate solutions to PDEs. The simplest case is the one-dimensional Kuramoto–Sivashinsky system, which comes in a spatially inhomogeneous and homogeneous variant. In the former case, i.e., in the face of an explicit dependence on the spatial coordinate, we show how one can use positional encodings to effectively train an ELM. More complex systems such as (<ref>) for two spatial dimensions are shown to be learnable; moreover, we demonstrate the exploitability of multiple inherent symmetries.

Let us also stress that what we propose is not the same as convolutional extreme learning machines [30]. These models initialize the convolution kernels of CNNs in a certain way (e.g., randomly or with Gabor filters) and only train the last layer in an ELM-style fashion.

§ PREDICTING (IN)HOMOGENEOUS PDES

[thin, gray] (0,0) – (0,7*ß);
[thin, gray] (ß,0) – (ß,7*ß);
ıin 0,...,7
[thin, gray] (0,ı*ß) – (ß,ı*ß);

[fill=cyan] (0, 2*ß) rectangle ++(ß,2*ß) ;
[very thick, black] (0,ß) – (0,5*ß);
[very thick, black] (ß,ß) – (ß,5*ß);
ıin 1,...,5
[very thick, black] (0,ı*ß) – (ß,ı*ß);

[thin, gray] (,0) – (,7*ß);
[thin, gray] (+ß,0) – (+ß,7*ß);
ıin 0,...,7
[thin, gray] (,ı*ß) – (+ß,ı*ß);

[fill=cyan] (, 2*ß) rectangle ++(ß,2*ß) ;
[very thick, black] (,2*ß) – (,4*ß);
[very thick, black] (+ß,2*ß) – (+ß,4*ß);
ıin 2,...,4
[very thick, black] (,ı*ß) – (+ß,ı*ß);

[thick] (*ß, 1.5*ß) rectangle
(*ß+*ß,4.5*ß) node[midway] ELM;

[thick, black] (ß,ß) – (*ß, 1.5*ß);
[thick, black] (ß,5*ß) – (*ß, 4.5*ß);
[thick, black] (*ß+*ß,1.5*ß) – (,2*ß);
[thick, black] (*ß+*ß,4.5*ß) – (,4*ß);

ıin 0,...,6
node at (-1.5*ß,ı*ß+0.5*ß)
$\overline x_{
\ifthenelse{\i<3}
{ {\the\numexpr(-2*(\i-3)-1)}/2}

ıin 2,...,3
node at (+2.5*ß,ı*ß+0.5*ß)
$\overline x_{
\ifthenelse{\i<3}
{ {\the\numexpr(-2*(\i-3)-1)}/2}

(*ß+*0.5*ß, +7*ß) circle (.4) node $p(\overline x)$;
[->] (*ß+*0.5*ß, +7*ß-0.4) – (*ß+*0.5*ß, 4.5*ß);

node at (0.5*ß,-ß)  $v$ ;
node at (+0.5*ß,-ß)  $\Phi(v)$ ;

[decorate, decoration = calligraphic brace]
(-3*ß,1*ß) –  (-3*ß,5*ß) node[midway, left] $\sl(v; \overline x, 2\ell+s)\, $ ;

[decorate, decoration = calligraphic brace,mirror]
(+ 4*ß,2*ß) –  (+4*ß,4*ß) node[midway, right] $\, \sl(\Phi(v); \overline x, s)$ ;

Sliding-window approach with an extent of $\ell = 1$ and a stride of $s = 2$. Because of the even window sizes, the centering coordinate $\overline x$ is exactly between two grid nodes. These nodes can be obtained by adding resp. subtracting $\Delta x/2$ from $\overline x$.

The Kuramoto–Sivashinsky equation has been a frequent test bed for learning algorithms. Most approaches are based on recurrent neural networks or reservoir computing [23, 28, 36].
Falling into the latter category, the method by [23] trains one (in case of spatial homogeneity) or multiple (in case of spatial inhomogeneity) reservoirs that handle different parts of the space. We employ a similar approach but always train a single ELM that is applied everywhere. Our implementation iterates over windows centering around different spatial coordinates $\overline x \in \R$ and predicts the PDE's state at the succeeding and smaller window. We use a stride $s$ that denotes the size of the output window, and an extent $\ell$ that denotes the number of nodes that the input window additionally takes into account left and right.
The ELM $E : \R^{2\ell + s} \rightarrow \R^s$ then computes the approximation

\begin{align*}
E(\sl(v; \overline x, 2\ell + s)) \approx \sl(\Phi(v); \overline x, s),
\end{align*}

see <Ref> for a graphic illustration of this approach.
If the PDE is spatially inhomogeneous, we found it best to work with an additional binary-like encoding of the position:

\begin{align*}
p(x) = (\cos(\pi x/L),\cos(2\pi x/L),\cos(4\pi x/L),\dots,\cos(2^k \pi x/L)),
\end{align*}

with the idea that the model could “branch” on the components of $p(x)$ whenever it needs to know its current location.

Performance of the ELM on the homogeneous Kuramoto–Sivashinsky system compared to the true trajectory.

The homogeneous variant.
The homogeneous Kuramoto–Sivashinsky system (<ref>) is simulated on a state-space of size $L = 200$ and a grid with $512$ nodes. After simulating the PDE in steps of $0.05\sec$ until it settles onto its attractor, we continued the simulation for another $200\sec$. We then trained an ELM without positional encodings, hidden dimension $H = 150$, softplus nonlinearity, extent of $\ell=7$ and stride of $s=4$ on the first $20$ snapshots (i.e., the first second). (That is, the ELM's input dimension is $2\ell+s = 18$ and its output dimension is $4$.) We iteratively and randomly pick one of these snapshots, select a random window of input-dimensional size, compute its embedding and the target window, and improve the estimates of $A$ and $G$ from (<ref>) using running averages. Additionally, we found that adding zero-mean Gaussian noise with standard deviation $10^{-4}$ to the input window substantially improved the model's stability.
The prediction error on a separate test set is not substantially reduced after 200 iterations, i.e., samples. Finally, we used the model to predict the PDE's trajectory–which means that we apply the same ELM to $512/4=128$ windows–and compared it to the true trajectory. The results are depicted in <Ref>.

Performance of the ELM with and without positional encodings on the inhomogeneous Kuramoto–Sivashinsky system compared to the true trajectory.

The inhomogeneous variant.
The parameters of the simulation are as in the previous section; only the spatially inhomogeneous term of (<ref>) is introduced with $\lambda = 50$ and $\mu = 0.05$. This additional term causes the typical worm-like patterns of the Kuramoto–Sivashinsky system to concentrate around and merge with four equidistant points located approximately at $0$, $L/4$, $L/2$ and $4L/3$. We train an ELM with $p(x)=(\cos(\pi x/L), \cos(2\pi x/L), \cos(4\pi x/L))$ and without positional encodings on the first $200$ snapshots (i.e., the first $10$ seconds). The remaining parameters stay the same when compared to the previous section. The results are depicted in <Ref>.

The ELM-based approach performs well on both variants. Due to the small number of parameters and the way the model is trained, the algorithm converges quickly. The prediction quality is high and comparable to state-of-the-art algorithms based on recurrent networks [23, 28, 36], but without requiring any initialization phase. Moreover, tuning of the hidden layer is much less challenging than selecting a proper reservoir. <Ref> shows that including the positional encoding allows the model to recover typical patterns.

§ SAMPLE-EFFICIENT LEARNING IN MULTIPLE SPATIAL DIMENSIONS

Trajectory of the two-dimensional Kuramoto–Sivashinsky system compared to the ELM's prediction.

We consider (<ref>) in two spatial dimensions.
To predict $\Phi(v)$ from $v$, our algorithm uses windows of side length $2\ell + s$ to predict windows of side length $s$. If $\overline x \in \R^2$ is the window's center, the coordinates of the input window take the form $\overline x_{i,j} := \overline x + \Delta x \big(\substack{i \\ j}\big)$, where both $i$ and $j$ range from $-(\ell+(s-1)/2)$ to $\ell + (s-1)/2$ in unit steps. The output window contains the coordinates where $i$ and $j$ range from $-(s-1)/2$ to $+(s-1)/2$ instead. (Notice that $i,j$ might not be integer.) More formally, we introduce the square of length $k$ centering around $\overline x$ for a $v \in \Fc$ as

\begin{align*}
\sq(v; \overline x, k) = \begin{bmatrix}
v(\overline x_{-(k-1)/2,-(k-1)/2}) & \dots & v(\overline x_{+(k-1)/2,-(k-1)/2}\big) ) \\
\vdots & \ddots & \vdots  \\
v(\overline x_{-(k-1)/2,+(k-1)/2}) & \dots & v(\overline x_{+(k-1)/2,+(k-1)/2}\big))
\end{bmatrix} \in \R^{k \times k}.
\end{align*}

The ELM $E : \R^{(2\ell+s) \times (2\ell+s)} \rightarrow \R^{s \times s}$ is trained to give the approximation

\begin{equation}
E ( \sq(v; \overline x, 2\ell+s) )
\approx \sq(\Phi(v); \overline x, s).
\label{eq:ELM-approximation}
\end{equation}

The system is simulated on a grid of size $256 \times 256$ with periodic domain defined on $(0,60\pi) \times (0,60\pi)$ in steps of $\Delta t = 0.01\sec$ until it settled onto its attractor. We then recorded the trajectory for another step, which gave rise to a single input-output pair that is used to train the ELM. The model has an extent of $\ell = 2$, stride of $s = 4$, $400$ hidden units with softplus activation function, and is trained by uniformly selecting random windows and updating the matrices in (<ref>) using running averages. For evaluation purposes, we ran a second simulation that recorded the trajectory for another $30\sec$.

The ELM's prediction together with the true trajectory and their deviation is shown     in <Ref>. The model performs well and achieves high prediction accuracy. Notice that there are $750$ frames between $0\sec$ and $7.5\sec$, so it can make accurate predictions even if frequently bootstrapped on itself.

§.§ Exploiting symmetries and enforcing equivariance

[->] (0,5) – (-2.75,4) node[midway, left, shift=(-0.3,0.1)] $\Kc_{A_1}$;
[->] (0,5) – (-1.75,4);
[->] (0,5) – (-1,4);
[->] (0,5) – (-.3,4);
[->] (0,5) – (+.3,4);
[->] (0,5) – (+1,4);
[->] (0,5) – (+1.75,4);
[->] (0,5) – (+2.75,4) node[midway, right, shift=(+0.2,0.1)] $\Kc_{A_8}$;
[inner sep=0pt] at (0,5)

ıin 0,...,7
[->] (-100pt + 12.5pt + ı*200pt/8, 3.5) – ++(0, -0.75-0.5) ;

node at (-100pt, 2.75) ELM;
node at (-70pt, 2.75) …;
node at (+100pt, 2.75) ELM;

ıin 0,...,7
[->] (-100pt + 12.5pt + ı*200pt/8, 2) – ++(0, -0.95);

node at (-100pt, 1.5) $\Kc_{A_1^{-1}}$;
node at (-70pt, 1.5) …;
node at (+100pt, 1.5) $\Kc_{A_8^{-1}}$;

[<-] (-16pt,-0.3) – (-2.75,0.5);
[<-] (-12pt,-0.2) – (-1.75,0.5);
[<-] (-8pt,-0.1) – (-1,0.5);
[<-] (-4pt,-0.05) – (-.3,0.5);
[<-] (4pt,-0.05) – (+.3,0.5);
[<-] (8pt,-0.1) – (+1,0.5);
[<-] (12pt,-0.2) – (+1.75,0.5);
[<-] (16pt,-0.3) – (+2.75,0.5);

node at (+75pt, 0) $\frac 1 8 \times \sum$;

[inner sep=0pt] at (0,3.5)

[inner sep=0pt] at (0,2)

[inner sep=0pt] at (0,.8)

[inner sep=0pt] at (0,-.5)

Multiple symmetries can be used to improve a prediction.

The Kuramoto–Sivashinsky system (<ref>) for two spatial dimensions admits many symmetries that can be exploited.
For $x \in \R^2$, we introduce transformations of the form

\begin{align*}
R(x) = \overline x + A(x-\overline x), \quad R^{-1}(x) = \overline x + A^{-1}(x - \overline x),
\end{align*}

where $A \in \R^{2\times 2}$ is an orthogonal matrix and $\overline x \in \R^2$ the center of a window. Orthogonality
essentially restricts $A$ to rotations and reflections, so $R$ rotates or reflects points across the center coordinate.
Additionally, we take $A$ to be one of the matrices representing the symmetry group of the square.
This particular group consists of eight elements: The identity element, rotations of 90, 180 and 270 degrees, and finally reflections along the horizontal, vertical and two diagonal axes; we denote them as elements of $G = \{ A_1,\dots, A_8 \}$.
This restriction on $A$ is convenient to work with since it simply permutes the neighboring coordinates around $\overline x$: One can implement a transformation $\Kc_{A}$ taking one square of neighbors to the other using axis-flips and transpositions:

\begin{equation}
\sq(v \circ R; \overline x, k) = \Kc_{A} \sq(v; \overline x, k).
\label{eq:transforming-the-square}
\end{equation}

For example, if $A$ is the $90$-degree rotation matrix, we have $R(\overline x_{i,j}) = \overline x_{-j,i}$.
So $\sq(v; \overline x, k)$ can be transformed into $\sq(v \circ R; \overline x, k)$ by transposing and then inverting the first axis (cf. [3]).

When it comes to the Kuramoto–Sivashinsky system, notice that the operators $\Laplace$, $\Laplace^2$ and $| \nabla (\cdot)|$ are invariant with respect to $R$ (appendix, <Ref>). Consequently, we have $f(v) \circ R = f(v \circ R)$ for $v \in \Fc$.
This argument extends to the PDE's solution (appendix, <Ref>):

\begin{equation}
\Phi(v \circ R) = \Phi(v) \circ R.
\label{eq:invariance-solution}
\end{equation}

Intuitively, this establishes that transforming a state (rotating, reflecting, translating, etc.) and evolving it forward in time is the same as first evolving it and then applying the transformation.

Putting (<ref>), (<ref>) and (<ref>) together, we see that

\begin{align*}
\sq(\Phi(v); \overline x, s)
&= \sq(\Phi(v \circ R) \circ R^{-1}; \overline x, s) \\
&= \Kc_{A^{-1}} \sq(\Phi(v \circ R); \overline x, s) \\
&\approx \Kc_{A^{-1}} E( \sq(v \circ R; \overline x, 2\ell + s)) \\
&= \Kc_{A^{-1}} E(  \Kc_{A} \sq(v; \overline x, 2\ell + s)),
\end{align*}

where $A \in G$.
We can get some practical use out of this in two ways: First, we can artificially add new samples by mapping known input-output pairs $(\sq_{\text{in}}, \sq_{\text{out}})$ to artificial samples $(\Kc_A \sq_{\text{in}}, \Kc_A \sq_{\text{out}})$. This decreases the overall number of “real” samples needed. Second, by using multiple symmetries of one sample to generate a single output and therefore improving the prediction quality. That is, we predict $\sq_{\text{out}}$ by averaging over $\Kc_{A^{-1}} E(  \Kc_{A} \sq_{\text{in}} )$, where $A$ ranges over some elements in $G$. If $A$ ranges over all eight, this enforces full equivariance. This averaging procedure is shown in <Ref>.

Symmetries can be used to increase data efficiency (left) and to improve the prediction quality (right). (One sample = one input window and one output window.)

<Ref> shows that both approaches work well. Adding artificially altered training samples significantly reduces the overall number of necessary samples, and using multiple symmetries for prediction improves the model's accuracy. Notice that both steps can be executed separately: It is possible to train a model on multiple symmetries but to make predictions using less, and vice versa.

§ CONCLUSION

Our proposed approach appears to be highly promising for the task of modeling PDEs. This is due to several reasons.
First, it has low data requirements because a single ELM is used to predict all parts of the state space using only local information. In the case of spatial inhomogeneities, we have shown how to keep the sliding-window approach by including additional positional information.
Second, it is not recurrent like previous state-of-the-art approaches, which significantly simplifies the setup including training process and prediction.
Third, we have shown that it is possible to exploit system-inherent symmetries to decrease the number of necessary training samples and to improve the prediction quality.

The authors acknowledge financial support by the project “SAIL: SustAInable Life-cycle of Intelligent Socio-Technical Systems” (Grant ID NW21-059D), which is funded by the program “Netzwerke 2021” of the Ministry of Culture and Science of the State of Northrhine Westphalia, Germany.

§ APPENDIX

Define $R(x) = Ax + d$ for $x \in \R^n$, where $A \in \R^{n\times n}$ and $d \in \R^n$. Set $J = \{1,\dots,n\}$. For $f : \R^n \rightarrow \R$ and $\ifat \in J^m$, denote $\partial_{\ifat}^m f := \partial_{\ifat_1,\dots,\ifat_m}^m f$.

If $f : \R^n \rightarrow \R$ is $m \geq 1$-times differentiable and $\ifat \in J^m$, we have

\begin{equation}
\partial^m_{\ifat} (f \circ R)(x)
= \sum_{\kfat \in J^m} (\partial^m_{\kfat} f \circ R)(x) \prod_{s=1}^m A_{\kfat_s,\ifat_s}.
\label{eq:thr-1}
\end{equation}

This is done by induction on $m$. For $m = 1$, notice that

\begin{align*}
\partial_i (f \circ R)(x)
&= \textstyle \lim_{\delta \rightarrow^+ 0} \textstyle\frac 1 \delta (f(R(x+\delta {\bm 1}_i))-f(R(x)) ) \\
&= \textstyle \lim_{\delta \rightarrow^+ 0} \textstyle\frac 1 \delta (f((Ax+d)+\delta A {\bm 1}_i)-f(Ax+d) ) \\
&= \textstyle \nabla_{A{\bm 1}_i} f(Ax+d)
=  (A{\bm 1}_i)^T \nabla f(Ax+d)
=  \sum_{j=1}^n A_{ji} (\partial_j f \circ R)(x),
\end{align*}

where $\nabla_{A \bm{1}_i}$ denotes the directional derivative with respect to $A \bm{1}_i$.
For $m \geq 1$, we have inductively for $\ifat \in J^{m+1}$, $\ifat' = (\ifat_1,\dots,\ifat_m)$, $i = \ifat_{m+1}$:

\begin{align*}
\partial_{\ifat}^{m+1} (f \circ R)(x)
&= \textstyle \partial_{i} \partial^{m}_{\ifat'} (f\circ R)(x) \\
&= \textstyle \partial_{i} \sum_{\kfat \in J^m}^n (\partial^m_{\kfat} f \circ R)(x) \prod_{s=1}^m A_{\kfat_s,\ifat'_s}  \\
&= \textstyle \sum_{\kfat \in J^m} \big( \sum_{k=1}^n A_{k, i} (\partial^{m+1}_{\kfat,k} f\circ R)(x)\big)\prod_{s=1}^m A_{\kfat_s,\ifat'_s} \\
&= \textstyle \sum_{\kfat \in J^{m+1}}  (\partial^{m+1}_{\kfat} f \circ R) (x) \prod_{s=1}^{m+1} A_{\kfat_s,\ifat_s}.
\end{align*}

If $f : \R^n \rightarrow \R$ is four times differentiable and $A A^T = I$, we have

\begin{align*}
| \nabla (f \circ R) | = | \nabla f \circ R |, \quad
\Laplace (f \circ R) = \Laplace f \circ R, \quad
\Laplace^2 (f \circ R) = \Laplace^2 f \circ R.
\end{align*}

Notice that $A A^T = I$ implies $\sum_{i=1}^n  A_{k,i} A_{j,i} = 1$ if $k = j$ and $0$ otherwise. Denote this fact by $(*)$.

\begin{align*}
|\nabla (f\circ R)(x) |^2
&= \textstyle \sum_{i=1}^n (\partial_i (f\circ R)(x))^2 \\
&= \textstyle \sum_{i=1}^n (\sum_{j=1}^n A_{ji} (\partial_j f \circ R)(x))^2
& \eqref{eq:thr-1} \\
&= \textstyle \sum_{i=1}^n \sum_{j_1,j_2=1}^n A_{j_1, i} A_{j_2,i} (\partial_{j_1} f\circ R)(x)(\partial_{j_2} f\circ R)(x) \\
&= \textstyle \sum_{j_1,j_2=1}^n (\partial_{j_1} f \circ R)(x)(\partial_{j_2} f \circ R)(x) \sum_{i=1}^n A_{j_1, i} A_{j_2,i} \\
&= \textstyle \sum_{j=1}^n  (\partial_{j} f \circ R)(x)^2 = | (\nabla f \circ R)(x)|^2. & (*) \\
\Laplace (f \circ R)(x)
&= \textstyle \sum_{i=1}^n \partial_{i,i}^2 (f\circ R)(x) \\
&= \textstyle \sum_{i=1}^n \sum_{j_1,j_2=1}^n A_{j_1,i} A_{j_2,i} (\partial^2_{j_1,j_2} f \circ R)(x)
& \eqref{eq:thr-1} \\
&= \textstyle \sum_{j_1,j_2=1}^n (\partial^2_{j_1,j_2} f \circ R)(x) \sum_{i=1}^n A_{j_1,i} A_{j_2,i} \\
&= \textstyle \sum_{j=1}^n (\partial^2_{j,j} f \circ R)(x) = (\Laplace f \circ R)(x).
& (*) \\
\Laplace^2 (f \circ R)(x)
&= \textstyle \sum_{i=1}^n \sum_{j=1}^n \partial_{i,i,j,j}^4 (f \circ R)(x) \\
&= \textstyle \sum_{i=1}^n \sum_{j=1}^n \sum_{\kfat \in J^4} A_{\kfat_1,i} A_{\kfat_2,i} A_{\kfat_3,j} A_{\kfat_4,j}
(\partial^4_{\kfat} f \circ R)(x)
& \eqref{eq:thr-1} \\
&= \textstyle \sum_{\kfat \in J^4} (\partial^4_{\kfat} f \circ R)(x)
\sum_{i=1}^n A_{\kfat_1,i} A_{\kfat_2,i}  \sum_{j=1}^n A_{\kfat_3,j} A_{\kfat_4,j} \\
&= \textstyle \sum_{k_1=1}^n \sum_{k_3=1}^n (\partial^4_{k_1,k_1,k_3,k_3} f \circ R)(x) = (\Laplace^2 f \circ R)(x).
& (*)
\end{align*}

Define the explicit Euler method for $v \in \Fc$, $t \geq 0$, $\delta > 0$ by:

\begin{align*}
\Phi^{\delta}_t(v) := \begin{cases}
v + t f(v) & t \leq \delta  \\
\Phi^{\delta}_{t-\delta}(v) + \delta f( \Phi^\delta_{t-\delta}(v) ) & t > \delta,
\end{cases}
\end{align*}

Suppose that it converges for every $v \in \Fc$ and $t \geq 0$: $\lim_{\delta \rightarrow^+ 0} \Phi^\delta_t(v) = \Phi_t(v)$.
Suppose moreover that $f(v \circ R) = f(v) \circ R$ for $v \in \Fc$.

\begin{equation}
\Phi_t(v \circ R) = \Phi_t(v) \circ R.
\end{equation}

We show $\Phi^\delta_t(v \circ R) = \Phi^\delta_t(v) \circ R$ for $\delta > 0$, $t \geq 0$, by induction on the up-rounded number of times $\delta$ fits into $t$; $M = \lceil t / \delta \rceil$. For $M \leq 1$, we have $1 \geq M \geq t/\delta$ and hence $\delta \geq t$. Using the invariance of $f$ and linearity of the composition operator,

\begin{align*}
\Phi^\delta_t(v \circ R) = (v \circ R) + t f(v\circ R)
= (v + t f(v)) \circ R
= \Phi^\delta_t(v) \circ R.
\end{align*}

For the induction hypothesis, suppose there is a fixed $M \geq 1$ such that $\Phi^\delta_t(v \circ R) = \Phi^\delta_t(v) \circ R$ for all $\delta > 0, t \geq 0$ with $\lceil t / \delta \rceil \leq M$. For the induction step, suppose $\lceil t / \delta \rceil = M + 1$. Then

2 \leq M+1 = \lceil t / \delta \rceil < t/\delta + 1,

and therefore $\delta < t$. Thus,

\begin{align*}
\Phi^\delta_t(v \circ R)
&= \Phi^\delta_{t-\delta}(v \circ R) + \delta f(\Phi^\delta_{t-\delta}(v \circ R)).
\end{align*}

By induction hypothesis,

\begin{align*}
\Phi^\delta_{t-\delta}(v \circ R) = \Phi^\delta_{t-\delta}(v)  \circ R
\end{align*}

since $\lceil (t-\delta) / \delta \rceil = \lceil t/\delta - 1 \rceil = \lceil t / \delta \rceil - 1 = M$. Using invariance of $f$ and linearity of the composition operator, we obtain:

\begin{align*}
\Phi^\delta_t(v \circ R)
= (\Phi^\delta_{t-\delta}(v) + \delta f(\Phi^\delta_{t-\delta}(v))) \circ R
= \Phi^\delta_t(v) \circ R.
\end{align*}

This concludes the induction.
Taking $\delta$ to zero, we finally have

\begin{align*}
\Phi_t(v\circ R)
&= \textstyle \lim_{\delta \rightarrow^+ 0} \Phi^\delta_t(v \circ R)  \\
&= \textstyle \lim_{\delta \rightarrow^+ 0} (\Phi^\delta_t(v) \circ R) \\
&= \textstyle (\lim_{\delta \rightarrow^+ 0} \Phi^\delta_t(v)) \circ R
= \textstyle \Phi_t(v) \circ R.
\end{align*}

[1]
M. M. Bronstein, J. Bruna, T. S. Cohen, and Veličković. P.
Geometric deep learning: Grids, groups, graphs, geodesics, and
arXiv:2104.13478, 2021.

[2]
Steven L. Brunton, Joshua L. Proctor, and J. Nathan Kutz.
Discovering governing equations from data by sparse identification of
nonlinear dynamical systems.
Proceedings of the National Academy of Sciences,
113(15):3932–3937, April 2016.

[3]
T. S. Cohen and M. Welling.
Group equivariant convolutional networks.
In International Conference on International Conference on
Machine Learning, volume 33, pages 2990–2999, 2016.

[4]
Shifei Ding, Xinzheng Xu, and Ru Nie.
Extreme learning machine and its applications.
Neural Computing and Applications, 25(3-4):549–556, September

[5]
V. Dwivedi and B. Srinivasan.
Physics Informed Extreme Learning Machine (PIELM)–A rapid method
for the numerical solution of partial differential equations.
Neurocomputing, 391:96–118, 2020.

[6]
C. Esteves, C. Allen-Blanchette, A. Makadia, and K. Daniilidis.
Learning SO(3) Equivariant Representations with Spherical CNNs.
In Proceedings of the European Conference on Computer Vision
(ECCV), 2018.

[7]
N. Geneva and N. Zabaras.
Modeling the dynamics of PDE systems with physics-constrained
deep auto-regressive networks.
Journal of Computational Physics, 403:109056, 2020.

[8]
L. Guastoni, J. Rabault, P. Schlatter, H. Azizpour, and R. Vinuesa.
Deep reinforcement learning for turbulent drag reduction in channel
The European Physical Journal E, 46(4), 2023.

[9]
J. Han, A. Jentzen, and W. E.
Solving high-dimensional partial differential equations using deep
Proceedings of the National Academy of Sciences,
115(34):8505–8510, 2018.

[10]
Guang-Bin Huang, Qin-Yu Zhu, and Chee-Kheong Siew.
Extreme learning machine: Theory and applications.
Neurocomputing, 70(1-3):489–501, 2006.

[11]
G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, and
L. Yang.
Physics-informed machine learning.
Nature Reviews Physics, 3(6):422–440, 2021.

[12]
K. Kunisch and S. Volkwein.
Control of the burgers equation by a reduced-order approach using
proper orthogonal decomposition.
Journal of Optimization Theory and Applications,
102(2):345–371, 1999.

[13]
Randall J. LeVeque.
Finite Difference Methods for Ordinary and Partial
Differential Equations: Steady-State and Time-Dependent Problems.
Society for Industrial and Applied Mathematics, January 2007.

[14]
Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and
A. Anandkumar.
Fourier neural operator for parametric partial differential
arXiv:2010.08895, 2020.

[15]
Zichao Long, Yiping Lu, Xianzhong Ma, and Bin Dong.
Pde-net: Learning pdes from data.
In International Conference on Machine Learning, pages
3208–3216. PMLR, 2018.

[16]
L. Lu, P. Jin, G. Pang, Z. Zhang, and G. E. Karniadakis.
Learning nonlinear operators via DeepONet based on the universal
approximation theorem of operators.
Nature Machine Intelligence, 3(3):218–229, 2021.

[17]
A. Mauroy.
Koopman operator framework for spectral analysis and identification
of infinite-dimensional systems.
Mathematics, 9(2495), 2021.

[18]
Mikael Mortensen.
Shenfun: High performance spectral Galerkin computing
Journal of Open Source Software, 3(31):1071, November 2018.

[19]
H. Nakao and I. Mezić.
Spectral analysis of the Koopman operator for partial
differential equations.
Chaos: An Interdisciplinary Journal of Nonlinear Science,
30(11), 2020.

[20]
Mark JL Orr et al.
Introduction to radial basis function networks, 1996.

[21]
S. Pandey, P. Teutsch, P. Mäder, and J. Schumacher.
Direct data-driven forecast of local turbulent heat flux in
Rayleigh–Bénard convection.
Physics of Fluids, 34(4):045106, 2022.

[22]
Sandeep Pandey and Jörg Schumacher.
Reservoir computing model of two-dimensional turbulent convection.
Physical Review Fluids, 5(11):113506, November 2020.

[23]
Jaideep Pathak, Brian Hunt, Michelle Girvan, Zhixin Lu, and Edward Ott.
Model-Free Prediction of Large Spatiotemporally Chaotic
Systems from Data: A Reservoir Computing Approach.
Physical Review Letters, 120(2):024102, January 2018.

[24]
S. Peitz, H. Harder, F. Nüske, F. Philipp, M. Schaller, and K. Worthmann.
Partial observations, coarse graining and equivariance in Koopman
operator theory for large-scale dynamical systems.
arXiv:2307.15325, 2023.

[25]
S. Peitz and S. Klus.
Koopman operator-based model reduction for switched-system control of
Automatica, 106:184–191, 2019.

[26]
S. Peitz, S. Ober-Blöbaum, and M. Dellnitz.
Multiobjective optimal control methods for the navier-stokes
equations using reduced order modeling.
Acta Applicandae Mathematicae, 161(1):171–199, 2019.

[27]
S. Peitz, J. Stenner, V. Chidananda, O. Wallscheid, S. L. Brunton, and
K. Taira.
Distributed control of partial differential equations using
convolutional reinforcement learning.
Physica D: Nonlinear Phenomena, 461:134096, 2024.

[28]
Mushegh Rafayelyan, Jonathan Dong, Yongqi Tan, Florent Krzakala, and Sylvain
Large-Scale Optical Reservoir Computing for Spatiotemporal
Chaotic Systems Prediction.
Physical Review X, 10(4):041037, November 2020.

[29]
Pu Ren, Chengping Rao, Yang Liu, Jian-Xun Wang, and Hao Sun.
PhyCRNet: Physics-informed convolutional-recurrent network
for solving spatiotemporal PDEs.
Computer Methods in Applied Mechanics and Engineering,
389:114399, 2022.

[30]
Iago Richard Rodrigues, Sebastião Rogério Da Silva Neto, Judith Kelner,
Djamel Sadok, and Patricia Takako Endo.
Convolutional Extreme Learning Machines: A Systematic Review.
Informatics, 8(2):33, May 2021.

[31]
C. W. Rowley, I. Mezić, S. Bagheri, P. Schlatter, and D. S. Henningson.
Spectral analysis of nonlinear flows.
Journal of Fluid Mechanics, 641:115–127, 2009.

[32]
P. J. Schmid.
Dynamic mode decomposition of numerical and experimental data.
Journal of Fluid Mechanics, 656:5–28, 2010.

[33]
L. Sirovich.
Turbulence and the dynamics of coherent structures Part I:
Coherent structures.
Quarterly of Applied Mathematics, XLV(3):561–571, 1987.

[34]
John C. Strikwerda.
Finite Difference Schemes and Partial Differential
Equations, Second Edition.
Society for Industrial and Applied Mathematics, January 2004.

[35]
C. Vignon, J. Rabault, J. Vasanth, F. Alcántara-Ávila, M. Mortensen,
and R. Vinuesa.
Effective control of two-dimensional Rayleigh–Bénard
convection: Invariant multi-agent reinforcement learning is all you need.
Physics of Fluids, 35(6):065146, 2023.

[36]
P.R. Vlachas, J. Pathak, B.R. Hunt, T.P. Sapsis, M. Girvan, E. Ott, and
P. Koumoutsakos.
Backpropagation algorithms and Reservoir Computing in Recurrent
Neural Networks for the forecasting of complex spatiotemporal dynamics.
Neural Networks, 126:191–217, June 2020.

[37]
Jian Wang, Siyuan Lu, Shui-Hua Wang, and Yu-Dong Zhang.
A review on extreme learning machine.
Multimedia Tools and Applications, 81(29):41611–41660,
December 2022.

[38]
Rui Wang, Karthik Kashinath, Mustafa Mustafa, Adrian Albert, and Rose Yu.
Towards Physics-informed Deep Learning for Turbulent Flow
In Proceedings of the 26th ACM SIGKDD International
Conference on Knowledge Discovery & Data Mining, pages
1457–1466, Virtual Event CA USA, August 2020. ACM.

[39]
S. Werner and S. Peitz.
Learning a model is paramount for sample efficiency in reinforcement
learning control of PDEs.
European Control Conference (preprint: arXiv:2302.07160), 2024.

[40]
Matthew O Williams, Ioannis G Kevrekidis, and Clarence W Rowley.
A data–driven approximation of the koopman operator: Extending
dynamic mode decomposition.
Journal of Nonlinear Science, 25:1307–1346, 2015.

[41]
N. Winovich, K. Ramani, and G. Lin.
ConvPDE-UQ: Convolutional neural networks with quantified
uncertainty for heterogeneous elliptic partial differential equations on
varied domains.
Journal of Computational Physics, 394:263–279, 2019.

[42]
K. Zeng and M. D. Graham.
Symmetry reduction for deep reinforcement learning active control of
chaotic spatiotemporal dynamics.
Physical Review E: Statistical Physics, Plasmas, Fluids, and
Related Interdisciplinary Topics, 104(1):04210, 2021.

[43]
K. Zeng, A. J. Linot, and M. D. Graham.
Data-driven control of spatiotemporal chaos with reduced-order neural
ODE-based models and reinforcement learning.
Proceedings of the Royal Society A: Mathematical, Physical and
Engineering Sciences, 478(2267):20220297, November 2022.