The Pendulum Data Set

The mathematical pendulum, integrated symplectically and lifted into $\mathbb{R}^4$. It is small, deterministic and needs no download, and it is the smallest system on which a GeometricMachineLearning.SymplecticAutoencoder has something real to do: the data are two-dimensional but they arrive in four dimensions, on a curved submanifold that no linear method recovers exactly.

Everything here is pendulum and angular_to_euclidean — one call to GeometricProblems and GeometricIntegrators, and one coordinate change. The data loader is GeometricMachineLearning's own; this package adds no adapter for it.

The system

The Hamiltonian is the one GeometricProblems.Pendulum defines,

\[H(\theta, p_\theta) = \frac{p_\theta^2}{2m\ell^2} + mg\ell\cos(\theta),\]

with the angle $\theta$, its conjugate momentum $p_\theta$, and the length, mass and gravitational acceleration $\ell$, $m$, $g$.

The pendulum hangs down at θ = π

The potential is $+mg\ell\cos\theta$, not $-mg\ell\cos\theta$, so it is at its minimum at $\theta = \pi$ and at its maximum at $\theta = 0$. The stable equilibrium is therefore $\theta = \pi$ and the inverted one is $\theta = 0$. Trajectories with $H < mg\ell$ librate about $\pi$, trajectories with $H > mg\ell$ rotate, and $H = mg\ell$ is the separatrix between them.

Generating the data

pendulum integrates a grid of initial conditions and hands back the GeometricSolutions.EnsembleSolution:

using GMLDatasets

solution = pendulum()
(length(solution), length(solution.t))

(100, 101)

That is a hundred trajectories of a hundred and one time steps each, spread over $\theta \in [0, 2\pi]$ and $p_\theta \in [-2, 2]$ — the grid GeometricProblems itself uses, which covers both libration and rotation. The keyword arguments are the six grid bounds, the physical parameters, the timespan and timestep, and the integrator. The default integrator is two-stage Gauss collocation, which is symplectic.

The lift into four dimensions

angular_to_euclidean replaces the angle by the position of the bob in the plane and $p_\theta$ by the corresponding Euclidean momentum:

\[q = (\ell\sin\theta,\; \ell\cos\theta), \qquad p = (p_\theta\cos\theta/\ell,\; -p_\theta\sin\theta/\ell).\]

The map is a symplectomorphism onto its image, and its image is the tangent bundle of the circle of radius $\ell$: two dimensions of data inside four, held there by $\|q\| = \ell$ and $q\cdot{}p = 0$.

data = angular_to_euclidean(solution)
size(data)

(4, 101, 100)

The rows are $(q_1, q_2, p_1, p_2)$ — the $q$-components first, then the $p$-components, which is what every symplectic architecture in GeometricMachineLearning assumes. The two constraints hold to rounding:

q, p = data[1:2, :, :], data[3:4, :, :]
(maximum(abs, sum(abs2, q; dims = 1) .- 1), maximum(abs, sum(q .* p; dims = 1)))

(2.220446049250313e-16, 2.220446049250313e-16)

euclidean_to_angular inverts the lift, up to the $2\pi$-periodicity of the angle.

Handing it to a network

There is no pendulum_data_loader here: the array above is already in the shape GeometricMachineLearning.DataLoader reads, so the constructor takes it directly.

using GeometricMachineLearning

dl = DataLoader(data; autoencoder = true, suppress_info = true)
(dl.input_dim, dl.input_time_steps, dl.n_params)

(4, 101, 100)

autoencoder = true is what marks the columns as independent samples rather than a time series, which is what an autoencoder wants. From here on it is an ordinary GeometricMachineLearning training run; scripts/pendulum/train_sae.jl is a complete one.

Train across the separatrix

The standard reduction $(\theta, p_\theta)$ cannot be represented globally by two ordinary real-valued coordinates because $\theta$ is periodic. That does not prevent an SAE from finding a different two-dimensional representation of a bounded part of the cylinder. For example, it can represent the angular direction by an annular geometry in its latent plane. scripts/pendulum/train_sae.jl therefore trains on the default bounded grid, which contains both librating and rotating trajectories and crosses the separatrix. The learned coordinates need not resemble $(\theta, p_\theta)$ to represent those trajectories faithfully.

Checking the integrator

pendulum_energy evaluates $H$ in either set of coordinates — on the canonical solution, on the lifted array, or on a pair of coordinate arrays. That the two agree is the statement that the lift preserved the dynamics:

energy = pendulum_energy(data)
maximum(abs, energy .- pendulum_energy(solution))

8.881784197001252e-16

And because the integrator is symplectic, the energy of each trajectory stays in a bounded band around its initial value rather than drifting away from it over the run:

maximum(abs, energy .- energy[1:1, :])

1.5812078757448944e-6

scripts/pendulum/plot_dataset.jl draws both of these — the phase portrait in each set of coordinates, and the energy error over time.

Other parameters

parameters is the $(\ell, m, g)$ named tuple that GeometricProblems takes, and it travels with the solution: angular_to_euclidean reads $\ell$ off the problem rather than assuming it is one, and pendulum_energy(solution) reads all three.

heavy = pendulum(; parameters = (l = 2.0, m = 3.0, g = 9.81),
                   qsamples = [4], psamples = [3], timespan = (0.0, 2.0))
maximum(abs, sum(abs2, angular_to_euclidean(heavy)[1:2, :, :]; dims = 1))

4.000000000000001

That is $\ell^2 = 4$, not 1: the bob is on the circle of radius two, because angular_to_euclidean took $\ell$ from the problem rather than assuming it away.

For the Euclidean methods of pendulum_energy the parameters are an argument, since a bare array does not carry them:

heavy_parameters = (l = 2.0, m = 3.0, g = 9.81)
maximum(abs, pendulum_energy(angular_to_euclidean(heavy), heavy_parameters) .-
             pendulum_energy(heavy))

2.842170943040401e-14

Zero to rounding — the two coordinate systems give the same energy away from the defaults as well, which is the statement that no factor of $\ell$, $m$ or $g$ was dropped in the lift.