{'fraction_non_alphanumeric': 0.029300518134715025, 'fraction_numerical': 0.0075647668393782385, 'mean_word_length': 4.782055122828041, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
KeyCLD: Learning Constrained Lagrangian Dynamics in Keypoint Coordinates from Images

We present KeyCLD, a framework to learn Lagrangian dynamics from images. Learned keypoints represent semantic landmarks in images and can directly represent state dynamics. Interpreting this state as Cartesian coordinates coupled with explicit holonomic constraints, allows expressing the dynamics with a constrained Lagrangian. Our method explicitly models kinetic and potential energy, thus allowing energy based control. We are the first to demonstrate learning of Lagrangian dynamics from images on the dm_control pendulum, cartpole and acrobot environments. This is a step forward towards learning Lagrangian dynamics from real-world images, since previous work in literature was only applied to minimalistic images with monochromatic shapes on empty backgrounds. Please refer to our project page for code and additional results: https://rdaems.github.io/keycld/


Introduction and Related Work
Learning dynamical models from data is a crucial aspect while striving towards intelligent agents interacting with the physical world. Understanding the dynamics and being able to predict future states is paramount for controlling autonomous systems or robots interacting with their environment. For many dynamical systems, the equations of motion can be derived from scalar functions such as the Lagrangian or Hamiltonian. This strong physics prior enables more data-efficient learning and holds energy conserving properties. Graydanus et al. [1] introduced Hamiltonian neural networks. By using Hamiltonian mechanics as inductive bias, the models respect exact energy conservation laws. Lutter et al. [2,3] pioneered the use of Lagrangian mechanics as physics priors for learning dynamical models from data. Cranmer et al. [4] expanded this idea to a more general setting. By modelling the Lagrangian itself with a neural network instead of explicitly modelling mechanical kinetic energy, they can model physical systems beyond classical mechanics. Zhong et al. [5] included external input forces and energy dissipation, and introduced energy-based control by leveraging the learned energy models. Finzi et al. [6] introduced learning of Lagrangian or Hamiltonian dynamics in Cartesian coordinates, with explicit constraints. This enables more data efficient models, at the cost of providing extra knowledge about the system in the form of a constraint function.
It is often not possible to observe the full state of a system directly. Cameras provide a rich information source, containing the full state when properly positioned. However, the difficulty lies in interpreting the images and extracting the underlying state. As was recently argued by Lutter and Peters [7], learning Lagrangian or Hamiltonian dynamics from realistic renderings remains an open challenge. The majority of related work [1,8,9,10,11] use a variational autoencoder (VAE) framework to represent the state in a latent space embedding. The dynamics model is expressed in this latent space. Zhong and Leonard [12] use interpretable coordinates, however their state estimator module needs full knowledge of the kinematic chain, and the images are segmented per object. All of these methods use minimalistic images, empty backgrounds and monochromatic shapes with very distinct colors, which are hard to extend to the general setting. Our method works with more realistic rendering, including 3D models, lighting, reflections, complex backgrounds and shadows. We hope this paves the way towards learning Lagrangian dynamics from real-world images.
Instead of using VAE inspired latent embeddings, our method leverages fully convolutional keypoint estimator models to observe the state from images. Because the model is fully convolutional, it is also translation equivariant. This can lead to a higher data efficiency. Objects can be represented with one or more keypoints, fully capturing the position and orientation. Zhou et al. [13] used keypoints for object detection, with great succes. Keypoint detectors are commonly used for human pose estimation [14]. More closely related to this work, keypoints can be learned for control and robotic manipulation [15,16]. Minderer et al. learn unsupervised keypoints from videos to represent objects and dynamics [17]. Jacques et al. [18] leverage keypoints for system identification and dynamic modelling. Jakab et al. [19] learn a keypoint representation unsupervised by using it as an information bottleneck for reconstructing images. The keypoints represent semantic landmarks in the images and generalise well to unseen data. It is the main inspiration for the use of keypoints in our work.

Contributions (1)
We introduce KeyCLD, a framework to learn constrained Lagrangian dynamics from images. We are the first to use learned keypoint representations from images to learn Lagrangian dynamics. We show that keypoint representations derived from images can directly be used as the state for constrained Lagrangian dynamics, expressed in Cartesian coordinates. (2) We show how to control constrained Lagrangian dynamics in Cartesian coordinates with energy shaping, where the state is estimated from images. (3) We are the first to demonstrate succesful empirical results on the pendulum, cartpole and acrobot system from dm_control. Images are rendered with 3D shapes with lighting effects, shadows, reflections and background. This outperforms previous work that only shows results on minimalistic images with monochromatic shapes on empty backgrounds and paves the way towards learning Lagrangian dynamics from real-world images.

Constrained Lagrangian Dynamics
For a dynamical system with m degrees of freedom, a set of independent generalized coordinates q ∈ R m represents all possible kinematic configurations of the system. The time derivativesq ∈ R m are the velocities of the system.
If the system is fully deterministic, its dynamics are described by the equations of motion, a set of second order ordinary differential equations (ODE): where u(t) are the external forces acting on the system. From a known initial value (q,q), we can integrate f through time to predict future states of the system. It is possible to model f with a neural network, and train the parameters with backpropagation through an ODE solver [20].
However, by expressing the dynamics with a Lagrangian we introduce a strong physics prior [3]: T is the kinetic energy and V is the potential energy of the system. For any mechanical system the kinetic energy can be described as where M(q) ∈ R m×m is the positive semi-definite mass matrix. Ensuring that M(q) is positive semi-definite can be done by expressing M(q) = L(q)L(q) , where L(q) is a lower triangular matrix. It is now possible to describe the dynamics with two neural networks, one for the mass matrix and one for the potential energy. Since both are only in function of q and notq, and expressing the mass matrix and potential energy is more straightforward than expressing the equations of motion, it is generally much more simple to learn dynamics with this framework. In other words, adding more physics priors in the form of Lagrangian mechanics, makes learning the dynamics more robust and data-efficient [2,3,4,7].
The Euler-Lagrange equations (4) allow transforming the Lagrangian into the equations of motion by solving forq. W is the external work done on the system, e.g. forces applied for control.
Adding external forces u ∈ R l to the model is possible with the introduction of the input matrix g ∈ R m×l . If the external forces and torques are aligned with the degrees of freedom q, g can be a diagonal matrix or even an identity matrix. More generally, if no prior knowledge is present about the relationship between u and the generalized coordinates q, g(q) : R m → R m×l is a function of q and can be modelled with a third neural network [5]. If the system is fully actuated l = m, if it is underactuated l < m.
Finzi et al. [6] showed that expressing Lagrangian mechanics in Cartesian coordinates x ∈ R k instead of generalised coordinates q ∈ R m has several advantages. The mass matrix no longer changes in function of the state, and is thus static. This means that a neural network is no longer required to model the mass matrix, simply the values in the matrix itself are optimized. The potential energy V (x) and input matrix g(x) are in function of x. Expressing the potential energy in Cartesian coordinates can often be simpler than in generalised coordinates. E.g. for gravity, this is simply a linear function.
The Euler-Lagrange equations require that the system is expressed in generalised coordinates, meaning that all possible values of q correspond to possible states of the system. Since we are now expressing the system in Cartesian coordinates, this requirement no longer holds. We need additionally a set of n holonomic constraint functions Φ(x) : R k → R n , where n is the number of constraints so that the degrees of freedom are correct: m = k − n (see Appendix A. 6 for examples of constraint functions). Deriving the equations of motion including the holonomic constraints yields (see Appendix A.1 for the full derivation and details): with D being the Jacobian operator. Since time derivatives of functions modelled with neural networks are no longer present, equation (6) can be implemented in an autograd library which handles the calculation of gradients and Jacobians automatically. See Appendix A.2 for details and the implementation of equation (6) in JAX [21].
Note that in equation (6) only the Jacobian of Φ(x) is present. This means that there is no need to learn explicit constants in Φ(x), such as lengths or distances between points. Rather that constant dinstances and lengths through time are enforced by DΦ(x)ẋ = 0. We use this property to our advantage since this simplifies the learning process.
The constraint function Φ(x) adds extra prior information to our model. Alternatively, we could use a mapping function x = F(q). This leads directly to an expression of the Lagrangian in Cartesian coordinates usingẋ = DF(q)q: from which the equations of motion can be derived using the Euler-Lagrange equations, similar to equation (6). In terms of explicit knowledge about the system, the mapping x = F(q) is equivalent to the kinematic chain as required for the method in Zhong and Leonard [12]. Using the constraint function is however more general. Some systems, such as systems with closed loop kinematics, can not be expressed in generalized coordinates q, and thus have no mapping function [22]. We therefore argue that adopting the constraint function Φ(x) is more general and requires less explicit knowledge injected in the model.

Relationship between Lagrangian and Hamiltonian
Both Lagrangian and Hamiltonian mechanics ultimately express the dynamics in terms of kinetic and potential energy. The Hamiltonian expresses the total energy of the system H(q, p) = T (q, p) + V (q) [1,8]. It is expressed in (q, p), the position and the generalized momenta, instead of generalized velocities. Using the Legendre transformation it is possible to transform L into H or back. We focus in our work on Lagrangian mechanics because it is more general [4] and observing the momenta p is impossible from images. See also Botev et al. [11] for a short discussion on the differences.
All images z i in the sequence are processed by the keypoint estimator model, returning each a set of heatmaps s i representing the spatial probabilities of keypoint positions. s i consists of m heatmaps s i k , one for every keypoint x i k , k ∈ {1, . . . , m}. The keypoint estimator model is a fully convolutional neural network, maintaining a spatial represention from input to output. This contrasts with a model ending in fully connected layers regressing to the coordinates directly, where the spatial representation is lost [12]. Because a fully convolutional model is equivariant to translation, it can better generalize to unseen states that are translations of seen states.
To distill keypoint coordinates from the heatmaps, we define a Cartesian coordinate system in the image (see for example Fig. 1). Based on this definition, every pixel p corresponds to a point x p in the Cartesian space. The choice of the Cartesian coordinate system is arbitrary but is equal to the space of the dynamicsẍ(ẋ, x, t, u) and the constraint function Φ(x) (see Section 2). We use spatial softmax over all pixels p ∈ P to distill the coordinates of keypoint x k from its probability heatmap: Spatial softmax is differentiable, and the loss will backpropagate through the whole heatmap since x k depends on all the pixels. Cartesian coordinates x k of the different keypoints are concatenated in vector x which serves as the state representation of the system. This compelling connection between image keypoints and Cartesian coordinates forms the basis of this work. The keypoint estimator model serves directly as state estimator to learn Lagrangian dynamics from images.
Similar to Jakab et al. [19], x acts as an information bottleneck, through which only the Cartesian coordinates of the keypoints flow to reconstruct the image with the renderer model. First, all x k are transformed back to spatial representations s k using Gaussian blobs, parameterized by a fixed σ.
The renderer model can more easily interpret the state in this spatial representation, as it lies closer to its semantic meaning of keypoints as semantic landmarks in the reconstructed image. Additionaly to the input s k , the renderer model learns a constant feature tensor (inspired by HoloGAN [23]), which it can use to learn static background information. The keypoint estimator model and renderer model have similar architectures, utilizing down-and upsampling and skip connections as in [24,25] (See Appendix A.4 for details). Finally, a reconstrucion loss is formulated over the reconstructed images z i and original images z i .
Dynamics loss function The sequence {x i }, corresponding to the sequence of given images {z i }, and the constant input u is used to calculate the dynamics loss. A fundamental aspect in learning dynamics from images is that velocities can not be directly observed. A single image only captures the position of a system, and contains no information about its velocities 1 . Other work uses sequences of images as input to a model [8], or a specific velocity estimator model trained to estimate velocities from a sequence of positions [26]. Zhong and Leonard [12] demonstrate that for estimating velocities, finite differencing simply works better.
We experimented with the approach of Zhong and Leonard, but experienced difficulty training the models. Predicting trajectories with an ODE solver during training was cumbersome, because of the error in velocity estimation. Since the velocities are approximated with finite differencing, we found a collocation technique for learning the dynamics more suiteable. Additionaly, this approach is significantly faster than making predictions using an ODE solver. Roesch et al. [27] employ a similar method for first degree neural ODE's.
Concretely, the dynamics loss function is expressed as: whereẍ is the dynamics model as expressed in equation (6) and h the timestep. Equations (11) and (12) estimate time derivatives with central method differentials. Formulating a finite difference over a longer timespan could be benificial for noisy data but was not necessary in our experiments. The total loss is the weighted sum of L r and L d , with a weighing hyperparameter λ: L = L r + λL d .
To conclude, the keypoint estimator model, renderer model and dynamics models (mass matrix, potential energy and input matrix) are jointly trained end-to-end on sequences of images {z i } and constant inputs u with stochastic gradient descent.
Rigid bodies as rigid sets of point masses By interpreting a set of keypoints as a set of point masses, we can represent any rigid body and its corresponding kinetic and potential energy. Additional constraints are added for the pairwise distances between keypoints representing a single rigid body [6]. For 3D systems, at least four keypoints are required to represent any rigid body [28]. We focus in our work on 2D systems in a plane parallel to the camera plane. 2D rigid bodies can be expressed with a set of 2 point masses, which can further be reduced depending on the constraints and connections between bodies (see Appendix A.3 for more detail and proof). In our framework, the keypoint model is free to choose the relative placement of keypoints on the different moving parts of the dynamic system, enabling the choice of distinct landmarks that also express the state accurately, e.g. the endpoint of a beam.
The interpretation of rigid bodies as sets of point masses allows expressing the kinetic energy as the sum of the kinetic energies of the point masses. Corresponding with equation (3), the mass matrix for a 2D system is defined as a diagonal matrix with masses m k for every keypoint x k . To enforce positive values, the masses are parameterized by their square root and squared.

Control
A major argument in favor of expressing dynamics in terms of a mass-matrix and potential energy is the straightforward control design via passivity based control and energy shaping [29]. Recent works of Zhong et al. [5,12] use energy shaping in generalized coordinates. In Cartesian coordinates, energy shaping can still be used. This is easily seen from the fact that for the holonomic constraints Φ(x) ≡ 0, we have the derivative DΦ(x)ẋ = 0, which means that the constraint forces in equation (6) are perpendicular to the path and hence do no work nor influence the energy [30].
Energy shaping control makes sure that the controlled system behaves according to a potential energy Changing the behavior of the kinetic energy is also possible [31], but was deemed not necessary for this work. The passivity-based controller is (see Appendix A.5 for full derivation and details): Many model-based reinforcement learning algorithms require the learning of a full neural network as controller. Whilst in this work, due to knowledge of the potential energy, we only need to tune two parameters k p and k d .  Evaluating how well a dynamics model, learned from images, captures the true underlying system is an open research problem. In an effort to measure progress in this field, Botev et al. [11] recently published a benchmark with 17 datasets. None of these cover a control aspect or external forces, and the images are minimalistic monochromatic shapes with empty backgrounds. This is similar to other reported results in literature [1,9,8,10]. All these works furthermore train a latent space embedding, which does not allow explicit interpretation of the state as spatial coordinates, and do not model external inputs or control actions.
Zhong and Leonard [12] use interpretable state representation and model external forces. Their coordinate-aware encoder explicitly interprets the state from the different bodies by using the known kinematic chain, and it requires seperate images per rigid body. We use the more general constraint Lagrangian approach which does not need the kinematic chain (see Section 2) and more importantly, does not need seperate images for each body.
We adapted the pendulum, cartpole and acrobot environments from dm_control [32,33] with minor changes for our experiments. See Appendix A.6 for details about the environments, their constraint functions and the data generation procedure. The exact same model architectures, hyperparameters and control parameters were used across the environments. This further demonstrates the generality and robustness of our method. See Appendix A.7 for more details.
Future frame prediction Fig. 4 shows long-term predictions of images with trained KeyCLD models. See Appendix A.8 for more extensive qualitative results and the project page for movies.
Ablation study Various components of KeyCLD were ablated to investigate their effects. All models were trained according to the training scheme in Section 3. Hyperparameters and other training details were identical for all ablations.

KeyCLD
The full framework as described in Sections 2 and 3.

KeyCLD-NC
The constraint function is omitted.

NODE2
A second order neural ODE modelling the acceleration is used instead of the Lagrangian. The keypoint estimator and renderer model are identital to KeyCLD. See A.4 for details about the model architecture.

FC
The dynamical model is identical to KeyCLD, but the keypoint estimator and renderer model represent coordinates directly, respectively ending and starting with fully connected layers. See A.4 for details about the model architectures.
Since KeyCLD is trained directly on image observations, quantitative metrics can only be expressed in the image domain. The mean square error (MSE) in the image domain is not a good metric of long term prediction accuracy [17,12]. A model that trivially learns to predict a static image, which is the average of the dataset, learns no dynamics at all. Yet, this model could report a lower MSE than a model that did learn the dynamics but started drifting from the groundtruth after some time.  Therefore, we adopt the valid prediction time (VPT) score from Botev et al. [11]. VPT measures how long the predicted images stay close to the groundtruth images of a sequence: where z i are the groundtruth images, z i are the predicted images and is the error threshold. is determined seperately for the different environments because it depends on the relative size in pixels of moving parts. It is equal to the MSE of the averaged image of the respective validation dataset. Thus it is the lower bound for a model that would simply predict a static image. Fig. 3 shows the MSE over time for all systems and ablations. Error bands were calculated by sampling 20 sequences in the validation set.
The VPT metric in Table 1 and MSE results in Fig. 3 indicate that KeyCLD managed to model cartpole and acrobot very accurately. Interestingly, results for pendulum are worse. Though it is a more simple dynamical system, the pendulum has a thinner pole than the cartpole and the acrobot. Interpretable energy models and control Since KeyCLD models the potential energy explicitly, we can use it for energy shaping control. See Fig. 5 for results of succesfull swing-up of the pendulum, cartpole and acrobot system. The same control parameters k p = 5.0 and k d = 2.0 are used for all systems, demonstrating the generality of the control method.

Conclusion and Future Work
We introduced the use of keypoints to learn Lagrangian dynamics from images. Learned keypoint representations derived from images can be directly used as positional state vector for learning constrained Lagrangian dynamics. We are the first to demonstrate succesful empirical results on dm_control environments, with more realistically rendered images than other results in literature.
Our work focusses on 2D systems, where the plane of the system is parallel with the camera plane. Elevation to 3D, e.g. setups with multiple cameras, is an interesting future direction. Secondly, our work focusses on energy-conserving systems. Modelling energy dissipation is necessary for realworld applications. Several recent papers have proposed methods to incorporate energy dissipation in the Lagrangian dynamics models [5,34]. However, Gruver et al. [35] argue that modelling the acceleration directly with a second order differential equation and expressing the system in Cartesian coordinates, is a better approach. Further research into both approaches would clarify the benefit of Lagrangian and Hamiltonian priors on real-world applications.

Broader impact
A tenacious divide exists between control engineering researchers and computer science researchers working on control. Where the first would use known equations of motion for a specific class of systems and investigate system identification, the latter would strive for the most general method with no prior knowledge. We believe this is a spectrum worth exploring, and as such use strong physics priors as Lagrangian mechanics, but still model e.g. the potential energy with an arbitrary neural network. The broad field of model-based reinforcement learning could benefit from decades of theory and practice in 'classic' control theory and system identification. We hope this paper could help bridge both worlds.
Using images as input is, in a broad sense, very powerful. Since camera sensors are consistently becoming cheaper and more powerful due to advancements in technology and scaling opportunities, we can leverage these rich information sources for a deeper understanding of the world our intelligent agents are acting in. Image sensors can replace and enhance multiple other sensor modalities, at a lower cost. This work demonstrates the ability to efficiently model and control dynamical systems that are captured by cameras, with no supervision and minimal prior knowledge.

A.1 Derivation of constrained Euler-Lagrange equations
The Lagrangian of a mechanical system described in Cartesian coordinates x ∈ R k is: with M a static mass matrix, not depending on x, and V (x) the potential energy. If the system has m degrees of freedom, additionally n holonomic constraints are necessary such that m = k − n. These are described by a constraint function Φ(x) : R k → R n . Including the input matrix g(x) ∈ R k×l and external inputs u(t) ∈ R l , the constrained Euler-Lagrange equations can be expressed with a vector λ(t) ∈ R n containing Lagrange multipliers for the constraints [6,30]: Because the mass matrix is static 2 , this is simplified to: The Lagrange multipliers λ(t) are solved by substitutingẍ from equation (20) in equation (21): We use the chain rule a second time to get rid of the time derivative of DΦ(x): Substituting λ(t) in (20) we arrive at equation (6).

A.2 Implementation of constrained Euler-Lagrange equations in JAX
It could seem a daunting task to implement the derivation of the constrained Euler-Lagrange equations (6) in an autograd library. Therefore, we provide an implementation in JAX [21]. For more context, please see the full code base on the project page (keycld/models.py).  Proof. To find conditions such that the kinetic energy expressed in two point masses should be equal to the rigid body representation, we start by expressing general 3D-movement: Where the vector xc are the coordinates of the center of mass and the vector x i/c is the position of the point mass relative to the center of mass. Since this relative position x i/c has fixed length, only a rotation is possible and hence the equation of the velocity is: where ω is the rotational velocity of the body. Substituting this in the kinetic energy of the point masses, we get: Where we calculated the square and used the circular shift property of the triple product on the last term.
For movement in the 2D-plane (i.e. ω = ezωz and xi = exxi,x + eyxi,y), this becomes: Matching the kinetic energy of the 2 point masses (equation (27)) with that of the rigid body representation (left hand side of Figure 6), we get following conditions: Since the last equation is a vector equation, this gives us four equations in six unknowns (m1,m2,x1,x,x1,y,x2,x,x2,y), which leaves us the freedom to choose two.
It follows from the third condition of (28) that points x1, x2 and xc should be collinear. To conclude, we can freely choose the positions of the point masses (as long as xc is on the line between them), and will be able to model the rigid body as a set of two point masses. In practice, KeyCLD will freely choose the keypoint positions to be able to model the dynamics. Depending on the constraints in the system, it is possible to further reduce the number of necessary keypoints. See Appendix A. 6 for examples.

A.4 Neural network model architectures
Neural network model architectures below are given in Flax code [36]. For more context, please see the full code base on the project page (keycld/models.py).

A.5 Energy shaping control
Energy shaping control makes sure that the controlled system behaves according to a potential energy V d (x) instead of V (x): where ypassive can be any passive output, the easiest choice being ypassive = k d g ẋ, where k d is a tuneable control parameter. The proposed potential energy V d should be such that: Where g ⊥ is the left-annihilator of g, meaning that g ⊥ g = 0. For fully actuated systems, the first condition of equation (30) is always met and the easiest choice is: where kp is a tuneable control parameter. The desired equilibrium position x * is obtained by putting a picture of the desired position of the system into the keypoint estimator model. Finally, the passivity-based controller that is used is: A.6 Details about the dm_control environments and data generation We adapted the pendulum, cartpole and acrobot environments from dm_control [32] implemented in MuJoCo [33]. Both are released under the Apache-2.0 license. Following changes were made to the environments to adapt them to our use-case 3 : Pendulum The camera was repositioned so that it is in a parallel plane to the system. Friction was removed.
Torque limits of the motor are increased.
Cartpole The camera was moved further away from the system to enable a wider view, the two rails are made longer and the floor lowered so that they are not cut-off with the wider view. All friction is removed. The pole is made twice as thick, the color of the cart is changed. Torque limits are increased and actuation is added to the cart to make the system fully actuated.
Acrobot The camera and system are moved a little bit upwards. The two poles are made twice as thick, and one is changed in color. Torque limits are increased and actuation is added to the upper part to make the system fully actuated.
Data generation For every environment, 500 runs of 50 timesteps are generated with a 10% validation split. Initial state for every sequence is at random position with zero velocity. The control inputs u are constant throughout a sequence, and uniform randomly chosen between the force and torque limits of the input. We set u = 0 for 20% of the sequences. We found this helps the model to learn the dynamics better, discouraging confusion of the energy models with external actions. Figure 7: From left to right the pendulum, cartpole and acrobot dm_control environments. The respective constraint functions are given below each schematic.
The constraint function for each of the environments are given in Fig. 7. As explained in Appendix A.3, every rigid body needs to be represented by two keypoints. But due to the constraints it is possible to omit certain keypoints, because they do not move or coincide with other keypoints. As experimentally validated, we can thus model all three systems with a lower number of keypoints, where the number of keypoints equals the number of bodies.
Pendulum One keypoint is used to model the pendulum. The second keypoint of this rigid body can be omitted because it can be assumed to be at the origin. Due to the constraint function, this point will provide no kinetic energy since it will not move. Since the other keypoints position and mass is freely chosen, any pendulum can be modelled. The constraint function expresses that the distance l1 from the origin to x1 is fixed. The value of l1 in the implementation is irrelevant because it vanishes when taking the Jacobian.
Cartpole Two keypoints are used to model the cartpole. The constraint function expresses that x1 does not move in the vertical direction and the distance l1 between x1 and x2 is constant. Again, the values of l1 and l2 in the implementation are irrelevant.
Acrobot Two keypoints are used to model the acrobot. The constraint function expresses that lengths l1 and l2 are constant through time. Again, the values are irrelevant in the implementation.

A.7 Training hyperparameters and details
All models were trained on one NVIDIA GTX RTX 2080 Ti GPU. We used the Adam optimizer [37], implemented in Optax [38]. Training one model takes approximately 75 minutes. Hyperparameters for all reported experiments are listed in Table 2. We used the exact same hyperparameters for all the environments and did not tune them individually.

A.8 Additional experimental results
Here we present additional qualitative results. Please refer to the project page for movies.
Explicit potential energy models Since the potential energy V is explicitly modelled, we can plot values throughout sequences of the state space. A sequence of images is processed by the learned keypoint estimator model, and the states are then used to calculate the potential energy with the learned potential energy model. Absolute values of the potential energy are irrelevant, since the potential is relative, but we gain insights by moving parts of the system separately. See Figure 8 for results for the pendulum, Figures 9 and 10 for the cartpole and Figures 11 and 12 for the acrobot.  The pendulum makes a full rotation. As expected, the potential energy follows a smooth sinusoidal path throughout this sequence. The maximum value is reached when the pendulum is upright, and the minimum value is reached when the pendulum is down.
Future frame predictions Additional future frame predictions are presented here. We generate predictions of the 48 next frames, given the first 2 frames of the ground truth sequence. Please compare these qualitative results for the pendulum, cartpole and acrobot environment (respectively Figures 13, 14 and 15). Every fourth frame of the sequence is shown. See also the project page for movies of these sequences.

A.9 Integrating constrained Lagrangian dynamics
The equations of motion (6) are integrated through time with an ODE solver. Due to numerical errors, the constraint function can start to drift over longer time predictions [6]. Integration can be improved by using numerical integrators specifically designed for constrained ODE's [22,39]. In our work, training is done using finite differences and there is no need for an ODE solver. To generate predictions of future frames, we use an ODE solver which we adapted with a constraint solver. Within the differential equation, the state is optimized so thatΦ(x) = 0. See the full code base on the project page for more details. The pole is fixed and the cart moves from left to right. As expected, the change in potential energy in this sequence is very low (compare to Fig. 9 with the same axis). A horizontal movement has no impact on the gravity potential.  Figure 11: Potential energy of the trained KeyCLD model of the acrobot environment. The first link makes a full rotation, the second link is fixed relative to the first link. As expected, the potential energy follows a smooth sinusoidal path throughout this sequence.  Figure 12: Potential energy of the trained KeyCLD model of the acrobot environment. The first link is fixed and the second link makes a full rotation. Again, the potential energy follows a smooth sinusoidal path throughout this sequence. Please compare with Fig. 11, where both links are moving.
Here the potential energy changes less, because the first link is not moving.  Figure 13: Predictions of future frames for the pendulum environments, for the different models in the ablation study. Top row is the ground truth sequence. Every predicted sequence is based on the first two frames of the ground truth sequence (a), since at least two frames are necessary to estimate the velocities. Every fourth frame of 50 frames in the sequences is shown. (b) KeyCLD is capable of making accurate long-term predictions, including reflections and shadow. Interestingly, the pendulum renders poor reconstructions in the bottom half. Since the constraint is not explicitly learned because only the constraint Jacobian is used in the loss, the model encodes a slightly moving keypoint position relative to the tip of the pendulum. Thus when integrating long predictions, the keypoint position reaches states that are never seen by the render model on the other half of the rotation. (c) Removing the constraint function still produces good short-term predictions but is detrimental for long-term predictions. The keypoint predictions start drifting and reach positions which the renderer model never encounters during training. (d) Modelling the dynamics with a second order neural ODE leads to predictions that are only accurate on the short-term. (e) Modelling the keypoints with fully connected layers instead of spatial heatmaps, yields static predictions. The model learns no dynamics at all.  Figure 14: Predictions of future frames for the cartpole environments, for the different models in the ablation study. Top row is the ground truth sequence. Every predicted sequence is based on the first two frames of the ground truth sequence (a), since at least two frames are necessary to estimate the velocities. Every fourth frame of 50 frames in the sequences is shown. (b) KeyCLD is capable of making accurate long-term predictions, including reflections and shadow. (c) Removing the constraint function still produces good short-term predictions but is detrimental for long-term predictions. (d) Modelling the dynamics with a second order neural ODE leads to predictions that are only accurate on the short-term. (e) Modelling the keypoints with fully connected layers instead of spatial heatmaps, yields static predictions. The model learns no dynamics at all.  Figure 15: Predictions of future frames for the acrobot environments, for the different models in the ablation study. Top row is the ground truth sequence. Every predicted sequence is based on the first two frames of the ground truth sequence (a), since at least two frames are necessary to estimate the velocities. Every fourth frame of 50 frames in the sequences is shown. (b) KeyCLD is capable of making accurate long-term predictions, including reflections and shadow. (c) Removing the constraint function still produces good short-term predictions but is detrimental for long-term predictions. (d) Modelling the dynamics with a second order neural ODE leads to predictions that are only accurate on the short-term for the first link, but the model fails to accurately represent the second link. (e) Modelling the keypoints with fully connected layers instead of spatial heatmaps, yields static predictions. The model learns no dynamics at all.