Constraint-based graph network simulator
Authors: Yulia Rubanova
Alvaro Sanchez-Gonzalez
Tobias Pfaff
Peter Battaglia
Venue: ICLR 2022 Submitted (2022)
Keywords: Physical simulations
graph neural network

Abstract: In the rapidly advancing area of learned physical simulators, nearly all methods train a forward model that directly predicts future states from input states. However, many traditional simulation engines use a constraint-based approach instead of direct prediction. Here we present a framework for constraint-based learned simulation, where a scalar constraint function is implemented as a trainable function approximator, and future predictions are computed as the solutions to a constraint satisfaction problem. We implement our method using a graph neural network as the constraint function and gradient descent as the constraint solver. The architecture can be trained by standard backpropagation. We test the model on a variety of challenging physical domains, including simulated ropes, bouncing balls, colliding irregular shapes and splashing fluids. Our model achieves better or comparable performance to top learned simulators. A key advantage of our model is the ability to generalize to more solver iterations at test time to improve the simulation accuracy. We also show how hand-designed constraints can be added at test time to satisfy objectives which were not present in the training data, which is not possible with forward approaches. Our constraint-based framework is applicable to any setting in which forward learned simulators are used, and more generally demonstrates key ways that learned models can leverage popular methods in numerical methods.

Under review as a conference paper at ICLR 2021
CONSTRAINT-BASED GRAPH NETWORK SIMULATOR
Anonymous authors Paper under double-blind review
ABSTRACT
In the rapidly advancing area of learned physical simulators, nearly all methods train a forward model that directly predicts future states from input states. How- ever, many traditional simulation engines use constraint- approach instead of direct prediction. Here we present a framework for constraint-base learned simulation, where scalar constraint function is implemented tion approximator, and future predictions are computed as the solutions to con- straint satisfaction problem. We implement our method using a graph neural net- work as the constraint function and gradient descent as the constraint solver. The architecture can be trained by standard backpropagation. We test the model on a variety of challenging physical domains, including simulated ropes, bouncing balls, colliding irregular shapes and splashing fluids. Our model achieves better comparable performance to top learned simulators. A key advantage of our model is the ability to generalize to more solver iterations at test time to improve the simulation accuracy. We also show how hand-designed constraints can be added at test time to satisfy objectives which were not present in the training data, which is not possible with forward approaches. Our constraint-based framework is applicable to any setting in which forward learned simulators are used, and more generally demonstrates key ways that learned models can leverage popular methods in numerical methods.
trainable fune-
INTRODUCTION
Consider a bowling ball colliding with bowling pin. You might explain this event involving pair of forces being generated, one which causes the pin move, and the other which causes the ball to careen away with different direction and speed. This kind of intuitive cause-and-effect approach is analogous to physical simulators that apply an explicit forward model to calculate a future state directly from the current one, such as when numerically integrating discretized equations of motion.
An alternative, but equally valid, way to explain the collision is in terms of constraint satisfaction: the ball and pin cannot occupy the same location at the same time, and their combined energies and momenta must be conserved, so the post-collision trajectories are the only way the future can unfold without violating these constraints. This constraint-based approach is analogous physical simulators that use an implicit function to model a system of constraints over the current and future states, and which generate a prediction by searching for future state that respects all constraints.
Both families of simulators-those based on explicit, forward functions versus those which define the dynamics implicitly, via constraints-are widely used in physics, engineering, and graphics. In principle they can model the same types of dynamics, however they differ in how their respective predictions are computed and practice strike different trade-off's that determine why one or the other is preferred in different domains. For example, explicit methods are popular for large systems with (mostly) independent local effects whose space and time derivatives are relatively smooth, and their accuracy can often be increased by discretizing space and time more finely. Implicit approaches are often preferred for systems with strong interactions, such rigid and stiff dynamics, and more accurate solutions can often be found by using more sophisticated constraint solvers or by increasing the computational budget (e.g., solver iterations) allocated to searching for solutions. In machine learning (ML). there have been rapid advances recently in methods for learning to simulate complex dynamic processes, however almost all (c.g., Sanchez-Gonzalez et al. (2020); Pfaff et al. (2021)) have focused on explicit forward model approaches, with few exceptions (Yang al., 2020).

---

Under review as a conference paper at ICLR 2021
Here present framework for learning to simulate complex dynamics via constraint satisfaction. Our "Constraint-beee based Graph Network Simulator" (C-GNS) defines a single scalar-valuee constraint function that represents whether future state satisfies the physical constraints, conditioned on the current and previous states. The constraint function is implemented as a Graph Neural Network (GNN) (Bronstein et al., 2017; Battaglia al., 2018), which can model systems with rich composi- tional structure- munntture-uulipll bodies, complex meshes, etc. To predict the next state via the constraint function's implicit representation of the dynamics, a gradient-based solver finds a proposed state which satisfies the constraints. We train it through the solver by backpropagation. We also intro- duce a hybrid approach that proposes and refines the future state using an explicit iterative predictor, rather than solving for learned constraints.
We tested the C-GNS on a variety of challenging physical simulation domains generated by sev- eral different simulation engines: simulated rope, bouncing balls, and bouncing irregular rigid shapes (MuJoCo, Todorov et al. (2012)) and splashing fluids (Flex, Macklin et al. (2014)). We found that the C-GNS's simulated rollouts were more accurate than a state-of-the-ar Graph Net Simula- tor (GNS) (Sanchez-Gonzaler et al., 2020) with comparable number of parameters. At test time, the C-GNS could use additional solver iterations to improve its predictive accuracy, striking desired speed- speed-accuracy trade-offs. It could also satisfy new, dand-desinned constraints jointly alongside its learned constraints. Neither of these capabilities are possible in explicit forward-style approaches.
2 BACKGROUND AND RELATED WORK
Constraint solvers are central to many physics simulators. Most rigid-body and game engines use constraints to model joints, collision and contact (Baraff, 1994). They are used for limiting strain in realistic cloth simulation (Thomaszewski et al., 2009), and are core component in Eulerian incompressible fluid solvers to solve for pressure (Chorin, 1967). Recently, position-based (Müller al., 2007) and projective dynamics methods (Bouaziz al., 2014) have become very popular for interactive simulation. These methods express dynamics purely as constraints, and can simulate a wide range of physical systems from rigids over soft-bodies fluids (Macklin et al., 2014).
Machine learning methods for accelerating scientific simulation of complex systems, such as tur- bulence (Kochkov et al., 2021; Wang et al., 2020) and aerodynamics (Thuerey et al., 2020; Zhang et al., 2018), have grown rapidly in recent years. LNN-based learned simulators, in particular, are a very flexible approach which can model a wide range of systems, from articulated dynam- ics (Sanclez-Gonzalez et al., 2018) to particle-based physics (Mrowca et al., 2018: 8: Li et al., 2019; et al., 2020) and mesh-based continuum systems (Pfaff et al., 2021; De Avila Belbute-Pere PPrree et al., 2020), and generalize well to unseen scenarios. Combining learning algo- rithms with principles from physics and numerical methods, such as auxiliary loss terms and rich inductive biases, can improve sample complexity, computational efficiency, and generalization (Wu et al., 2018; Karniadakis et al., 2021; Chen et al., 2018; Rubanova et al., 2019). Imposing Hamilto- nian (Greydanus et al., 2019; Sanchez- Sanchez-Gonzalez et al., 2019; Chen et al., 2019) and Lagrangian (Lut- ter et al., 2019; Cranmer et al., 2020; Finzi et al., 2020) mechanics in learned simulators offers unique speed/accuracy tradeoffs and can preserve symmetries more effectively.
Recent methods have been proposed for learning constraint functions and solving them in a model's forward pass (Duvenaud et al. (2020)'s "Deep Implicit Layers" tutorial is an excellent hands-on survey). Such models can play games (Amos & Kolter, 2017; Wang et al., 2019), optimize power flow (Donti et al., 2021), support robotic planning (Loula al., 2020), and perform combinato- rial optimization (Bartunov et al., 2020). Solvers such as gradient descent and Newton's 's method are differentiable, and support training by backpropagation, but this can be computationally expen- sive, so approaches such as Deep Equilibrium Models (DEM) (Bai et al., 2019; 2020) use implicit differentiation to compute gradients only at the solution point.
Despite the popularity of constraint-basee based traditional simulators, only single simulator which uses learned constraints has been reported (Yang et al., 2020). Their "Neural Projections" method, based on Goldenthal et al. (2007), iteratively proposes future state with an Euler step, then projects the proposal onto a learned constraint manifold, implemented as a multilayer perceptron (MLP). Crucially, their constraint function only measures how much individual state violates the learned constraints, and thus is not an implicit representation of the dynamics. It is suitable for quasi-statte regimes, but not scenarios such as the elastic collisions in the bowling ball example described above.
2

---

Under review as a conference paper at ICLR 2021
(a)
(b) Forward predictor
(c) Forward iterative predictor
(d) Constraint- based predictor
X
Xx
(d)
Xs1
y(i)
Iterator
Y(i)
Solver
y(
X,
Vy fc
fn
/DJ
fc
X
ONNY
GNN
GAN
Predictor
Updater
6Y
8Y = -V -yy-vo
y(-1) y(i + 8Y
y(+h)l y 6Y
X141
Y(M)
Y(M)
Figure I: Learned simulator schematics. (a) A simulator, 8, maps Xge to a future state X using a PREDICTOR that returns Y. which represents information about the system's temporal evo- lution, and an UPDATER which uses Ý to update X, to X,+1- (b) Forward GNN simulator. The PREDICTOR maps Xgr directly to Y using fD- (c) Iterative GNN simulator. The iterator refines Y by repeatedly applying fDI. (d) Constraint-based Graph Network simulator (C-GNS). The PREDIC- TOR iteratively solves for that satisfies a constraint function, fc. using Vy fc.
MODEL FRAMEWORK
Simulation basics A physical trajectory, measured at discrete time intervals, is sequence of states, (X1...., Xr), where X, represents properties such the positions, velocities, masses, etc, of ele- ments of the system. A physical simulator, is a function that maps current and/or previous state(s), which we term the context, X to predicted future state, Xe+1 a(X (see Figure la)². A simulated physical trajectory, termed a rollout, (Xe, X,+1, Xe42....) can be generated by repeat- edly applying its own predicted state, X s().))
Simulators are often comprised of PREDICTOR mechanism which maps the context Xgr to an update value Ý, that represents information about the system's temporal evolution at the current time. Then Y is used by UPDATER mechanism update the current state the next state: X.+1 U Y), e.g., updating current positions and velocities represented by X, with new velocities and accelcrations represented by Ý, to predict the next state.
Explicit simulators Across science, engineering, and graphics, a popular class of simulators are de- fined explicitly: the state update Y is predicted directly from X using explicit forward function, y = (x) as illustrated in Figure lb. Among the rapidly growing family of learned simulators, the forward function fo is typically implemented using neural network (Sanchez-Gonzalez et al., 2020: Pfaff et al.. 2021).
constraint-based implicit simulators Here we explore learned simulators based on implicit formu- lations of the dynamics. Rather than predicting the desired state directly, as explicit formulations, our implicit simulator uses differentiable constraint function, c = fc(X gt, Y), where e scalar that quantifies how well a proposed state update Ý agrees with X A future prediction is generated by applying solver, such as optimization zero-finding algorithm, to find a Y that satisfies the constraint function, and applying the UPDATER to update X, to X,+1. The fc represent all the physical constraints in the system, including the time dynamics.
Despite that physics Markovian, we use X as input because our framework can also apply to dynamic processes which are non-Matte Markovian. Providing previous states can also often be helpful when there are hidden properties of the system which are only identifiable over sequence of observed states, and when state does represent velocity momentum information.

---

Under review as a conference paper at ICLR 2021
As illustrated in Figure 1d, formulate constraint-solvine procedure via iterative method that starts with an initial proposal, Y On the i-th iteration, the solver uses the gradient of fc w.r.t Y at the current proposal to compute a change to the proposal, oY -- Vy fc( =y This 8Y is then used to revise the proposal to, y = y + oY. This process repeats for N steps, and the final proposal value is treated as the PREDICTOR's s output, Y = Y
Our constraint based model's fc is defined as trainable function approximator which is real-valued and lower bounded at zero, and uses gradient descent to find Y that minimizes it, where À is fixed step size. This induces the semantics that the desired Y = arg miny fc(Xst, (t) Y).
We also explore second constraint-solving procedure, inspired by Yang et al. (2020)'s Neural Projections' use of "fast projection" Goldenthal et al. (2007). Specifically, À
fe(Xg.,
Unlike gradient descent, fast projection is eero-finding algorithm, so in
|| Vy (exfexg..
this case fc is not lower bounded. This induces the semantics that fc(X 8) O.
This general formulation of constraint-based learned simulation can be trained by backpropagating loss gradients through the solver loop". The computational budget of the forward pass can be varied via the number of solver iterations N.
Explicit iterative simulators As a hybrid between forward and constraint-based simulators, we introduced a model which iteratively refines proposed state update, like in the constraint-based approach described above, but using an explicit function to directly output a ōY at each iteration, rather than solving a constraint function (see Figure Ic). See Section 4.3 for details.
4 EXPERIMENTS
4.1 EXPERIMENTAL. TASK DOMAINS
We test our framework on a variety of physical environments, shown Figure 2: ROPF, BOUNCING BALLS and BOUNCING RIGIDS, whose ground truth training and test data were generated by the MuJoCo physics simulator, as well as BOXBATH from Li et al. (2019). These environments demon- strate diverse set of physical constraints: 'hard' constraints (preserving the shape of the rigid object and resolving collisions), and 'soft' constraints on fluid movement, handling gravity and preserving momentum of the rope and bouncing balls. See the Supplementary Materials for details.
ROPE
BOUNCING BALLS
BOUNCING RIGIDS
BOxBATH
Figure 2: Renderings of the physical environments. Videos of the model rollours are available at: http:///wwwww/ostttaaaaaaeeedeid
4.2 MODEL. IMPLEMENTATIONS
Representing the physical system Our experimental domains are physical systems comprised of sets of interacting piint-like elements, c.g., objects, particles, mesh vertices, etc. We represent the state X, (p.)i=1..IIX.,, where [Xel is the number of elements, and y the j-th element's position at time t. There are also other static properties of the physical elements, e.g., masses, material types, etc., which we represent with Z to keep distinct from the dynamic state information represented by X, The input context is X = (2, X; X X Xt-1, X.).
Implicit differentiation at the solution point should he applicable as well, and potentially offer computa- tional henefits mentioned in the Scction 2, though
explore that here.
4

---

Under review as a conference paper at ICLR 2021
In our implementation, Y. represents the predicted changes in position (i.e., the "average velocity" across the time step)", i' Af/+1 p/+1 pl. The UPDATER then computes Xe+1 using p p + A where pf is provided in the input Xgt-
Constructing the input graph Our implementations of the fD- fDI, and fc use GNNs the func- tion approximators, so we need to pack the context, Xgt, and (for the fDI and fc) the proposed state update information, y into an input graph, G. = (Ve, E.). The edges E, represent pos- sible interactions among the elements, such as fully connected edges to represent collisions and rigid attachments in BOUNCING BALLS and BOUNCING RIGIDS, spring constraints in ROPE, and interactions among particles within fixed connectivity radius BOXBATH.
We enforced Translatoo-invaaiinnee by construction, by never providing absolute positions as input the models. Instead, the j-th input node's features are the static properties, and a sequence of the three most recent position changes (i.e. average velocities), vy [a], Ap' Api Api], where, App- = pl. -P---- For fDI and fc. which also take the solver's s current proposed Y we also concatenate the proposed average velocity from the i-th solver iteration, yo:(() pb. as input. For the input edge feature for an edge that connects from node j to k: we also provide the relative displacement vector between the nodes' positions, c/* = ,t 1/.
core We implemented fD. fDI. and fc using Graph Net- Pww---oooos--ooooo architecture. similar works (GN) (Battaglia et al., 2018), arranged in the Sannchez-Gozzlez et al., 2018; 2020; Pfaff et al., to previous work on GN-based learned simulators 2021). The Encoder uses two MLPs encode node and edge features into high-dimensional latent vectors. The Processor applies multiple GNs, with unshared weights, in sequence. with node and edge residual connections each step. We do not use global updates for the GNs. The Decoder uses an MLP to produce an output for each node.
GNN-based Sncode-Process-Der
The fo directly returns Y. The fDI returns a change to the proposed update 5Y for the current iteration. The fc's Decoder returns a scalar for each node to produce constraint value per node {eil = 1... [V] These node-wise wise constraint values are averaged to compute a single scalar c constraint for the entire system, fc(X 8) M |
Solving the constraint For for and fc we initialize (0) Ap/ to the most recent average velocity We used auto-differentiation in JAX to compute the gradient function, Vy fc. and the step size À was specific to the model variant, as described below. During training we used N solver iterations.
4.3 MODEL VARIANTS
The key questions in this work are whether constraint-basee learned simulators can compete with ex- plicit, forward learned simulators, whether implementing the constraint function with GNNs is more effective than with MLPs, and how minima- based constraint functions solved by gradient descent compare to constraints defined as the zeros of a function which are solved by fast projection (Gold- enthal et al., 2007). The following model variants allow us to answer these questions.
Forward GNN This is an explicit, forward GNN-based learned simulator based on the GNS models from Sanchez-Gonzalez et al. (2020); Pfaff et al. (2021). It directly predicts the state update Y from the past time points X
C-GNS Gradient Descent (C-GNS-GD) and C-GNS-Fast Projections (C-GNS-FP) These are our proposed constraint- based GNN models. For the C-GNS-GD, the scalar per-node output oi was squared, to force the overall fc to be non-negative, and gradient descent solver with fixed step size, À = 0.001, was used to minimize it. For C-GNS-FP, the À was based on "fast pro- jection" (Goldenthal et al., 2007: Yang et al., 2020), described in Section 3. Supplementary Figure B.5(c-d) shows ablations.
'For BOXBATH we vary number of modelling choices best match those in Sanchez-Gonzalez al. (2020). The major difference is that we set y to be the average acceleration rather than average velocity. See Supplementary Materials for other differences.
'To ensure analogous information is provided downstream of fb, the update rule also includes the previous average velocity: P4 pd + Apl + +0'
5

---

Under review as a conference paper at ICLR 2021
Rope
Bouncing Balls
Bouncing Rigids
Box Bath
slep position MSE
10"5
10-1
10-6
10"6
10-3
10-7
10-7
90
C-MLP-GD
C-MLP-FP
Forward GNN
10a
10a
T
9 10
C-GNS-GD
10a
10'
C-GNS-FP
Rollout position
10"
8 10
10a
10a
10-1
10"1
10
10-3
10
6x
Figure 3: Test MSE on the node positions predictions across different models. Top row: 1-step position MSE. Bottom row: full 160-step rollout MSE. The bar height represents the median MSEs over random seeds. The black cross marks show the MSE metric for each random seed. The black arrows indicate if MSE metric for random seeds exceeds the upper y limit of the figure.
Iterative GNN We implemented a hybrid between the Forward GNN and C-GNS, shown in Figure 1c. was identical to the C-GNS models, except its fDr directly predicted proposed state updates in fD. rather than being computed via the gradients was done with fc-
ConstraintMLP Gradient Descent (ConstraintMLPPGG) and ConstraintMLPP tions (ConstraintMILP-FF) These GNNs to implement fc, instead concatenated the embeddings of all the input nodes into a single vector and passed them to an MLP implementation of fc. By default, these models cannot han- dle variable-length inputs, so we padded smaller states with zeros up to the maximum state size. The ConstraintMLLP-FPP was the MLP analog to our C-GNS-FP, and was similar to Neural Projec- tions (Yang et al., 2020). The ConstraintMLP-GD used gradient descent, and was the MLP analog to our C-GNS-GD. We omit the results the ConstraintMLP models on BOXBATи (1024 nodes), as MLPs do not generally work well on physical systems with more than few particles (Battaglia et al., 2016; Goabllez et al., 2018).
Projec- MLP-based constraint models, which, rather than using
4.4 TRAINING AND EVALUATION
We trained the models to make next-step predictions, by computing the L2 loss between the pre- dicted X4+1 and the corresponding ground truth X1+1, averaged over nodes. All model weights and biases were trained using standard backpropagation with the Adam optimizer.
At test time, we compute 1-step metrics by evaluating the 1-step enrors along each point of the ground truth trajectory. We also evaluate rollout errors by iteratively applying the learned model starting from initial state, over 160 rollout steps, and computing the error between the predicted and ground truth trajectories.
5 RESULTS
Predictive accuracy Our experimental results show that our C-GNS-GD's performance was gen- erally better than the other model variants. Figure 3 compares the different models on 1-step and rollout position MSE (see Supplementary Table B.1 for numerical results). For each dataset,
Videos of the model rollouts available sites.googll..cmm/..... www.p//www...ooooeoooooootnnnnee

---

Under review as a conference paper at ICLR 2021
10
1 steps
MP steps
MP steps
MP steps
5 MP steps
Forward GNN
Iterative GNN
1.stop MSE
C-GNS-GD
10
Iterations
iterations
Iterations
Iterations
Iterations
steps
Figure 4: Test 1-step MSE error ROPE function of message-passing (MP) steps and solver iterations. The left five subplots shows performance for different numbers of message- passing steps. The blue and red bars represent the Iterative GNN and C-GNS-GD, respectively, where the bar height is Istepp MSE. The green bars show the Forward GNN (it does not use solver iterations). The rightmost subplot shows the Forward GNN with 1 10 MP steps.
used the same number of message-passing steps (MP) for all GN-based models. We used 2 MPs for the ROPE dataset, and MP for all other tasks.
The C-GNS-GD has lower 1-step MSE between the ground truth and predicted positions than other models across all datasets. Qualitatively, we observed that for Forward GNN with a single message- passing step, the box in BOXBATH "melts" over time, the forward model cannot preserve its rigid shape (see Videos). The comparable C-GNS-GD, by contrast, maintains the rigidity more effec- tively. These quantitative results suggest that constraint-basse learned simulators alternative to explicit, forward Iearned simulators. We generally found that the Iterative GNN was fairly competitive with the C-GNS-GD in overall performance and better than the Forward GNN.
competitive
We also found that the C-GNS-FP was generally less stable across seeds, and not as accurate as the C-GNS-GD. The same conclusion holds for ConstraintMLP-FF versus ConstraintMLPP speculate that the fast projection algorithm may make training challenging because the step size À proportional to fc, which may cause poor firoofiniing early training when the fc is not yet informative. Additionally, we find that C-GNS-FP algorithm becomes unstable in the areas with shallow constraint gradients, perhaps because À depends on the inverse of the gradient's norm.
We
We explored how varying the message-passileg steps and solver iterations (N) influenced the rel- ative performance among the models in our ROPE dataset. Figure 4 shows that the C-GNS-GD generally required fewer parameters and message-passing steps achieve comparable 1-step MSE to the other models. Supplementary Figure B.3 shows similar results for the rollout MSE. For most combinations of message- pasgging steps and number of solver iterations, C-GNS-GE GD (green) outper- forms the Iterative GNN (yellow), C-GNS-FP (purple) as well as the Forward GNN (blue) with the same number of MPs (the Forward GNN is not iterative model, so plot it single bar). We hypothesize that the solver iterations in the C-GNS and Iterative GNN may play a similar role to message passing with shared weights.
Interpreting the learned constraints To better understand the learned fc functions in the C-GNS- GD, Figure visualizes the node-wise constraint values function of Y (proposed average veloc- ity) for different nodes in the ROPE dataset while holding the other nodes' proposed update Y fixed.
Node
Node
Node
Node 4
Node 5
25
BRDUO 0
15
10
10
-0
Figure 5: Visualization of the constraint landscape for trained C-GNS-GD. Each subplot cor- responds to a different ROPE node. The heatmap's color shows the constraint value evaluated at different values of Y state. The plot centered around the ground-truth point. The colored points show the five iterations of the constraint solver, from the initial (0) (yellow) final (5) (orange).
7

---

Under review as a conference paper at ICLR 2021
(a) Ground truth sequence
(b) C-GNS
C-GNS with additional spatial constraints
Figure 6: Adding llan-designed constraints. (a) The ground truth sequence of rope states, initial- ized at the cyan-colored state, and simulated over 14 time steps, to the final, purple-colored state. (b) The C-GNS-GID GD's rollout, without added constraints. (c) The C-GNS-GD's rollout, with wall, floor, and disk disk-shaped "forbidden" zones, imposed at test time via hand-designed constraint functions. A video of the trajectories is available at Videos.
We also overlay the sequence of five points that represent the proposed where all nodes were jointly optimized. The figure shows the learned fc has a minimum near the ground truth Y, which the gradient descent steps are able to reach.
(+) steps from the solver
Incorporating novel constraints at test time We next explored a unique advantage the constraint-based model: because the fc measures the degree the physical constraints violated, we can incorporate additional, hesdnde constraints at test time, and use the model to poten- tially satisfy them. For the ROPE dataset, we designed three constraint functions that return positive values which increase quadratically as the rope enters different "forbidden" regions of the space: a vertical wall, horizontal floor, and disk-shaped region. We weighted these constraint terms by coefficient hyperparameter and added each of the hand-designed constraints to the learned fc term of C-GNS-Gi GD and ran the forward evaluation of the model.
As shown in Figure 6, the model was able to simulate the dynamics in a way that the corresponding forbidden region was avoided. In some cases, satisfying the joint constraint resulted in unintuitive behaviors, such as the rope links changing in length to adapt to the obstacle (Videos). However, this is to be expected, as the minimum of the joint constraint may not overlap with the minimum of the learned constraint, which is the one that would otherwise guarantee length preservation. For this example we added further hand-designed constraint that incentivizes maintaining relative distances between nodes. In general this is powerful example of how constraint- models can generalize outside their training data, and solve both for the learned dynamics and arbitrary desired constraints.
Generalizing to larger systems via increased solver iterations In principle, iterative and constraint- simulators should find more accurate solutions by increasing the number solver iterations, N. We investigated whether the C-GNS-GD and Iterative GNN trained on ROPE could generalize from Nitann = on which they were trained, to Nies E [0, 15). We also analyzed whether increased solver iterations could improve generalize performance from training on ropes with 5 10 nodes, to test ropes with 20 nodes.
Figure 7a (top row) shows that for ropes that match the 5- 10 nodes experienced during training. the Iterative GNN (light blue) overfits very heavily to Nest Nurin 5: error increases abruptly for N < and N 6. By contrast, the C-GNS-GD (light red) generalizes much better different Nust- Figure 7a (bottom row) shows that for test ropes with 20 nodes, the Iterative GNN again overfits, while the C-GNS-GD can generalize well longer ropes if Nias is increased.
We also trained the Iterative GNN and C-GNS-GD with additional loss terms that were applied to the (+) on each solver iteration, not only the final one, & = (N), We used an exponential decay factor, a = 0.25, which downweighted this additional loss term more heavily for earlier solver proposals. The dark blue and red curves in Figure 7a show how this additional loss further improves generalization to more solver iterations and larger systems as test time for the Iterative GNN, but especially the C-GNS-GI) Figure 7b visualizes how increasing the solver iterations systematically improves the quality of the long-term rollout accuracy in the ROPE dataset.
Together these results show the C-GNS-GD is effective in making use of additional resources at test time. This opens the exciting possibility of training on small, simple systems, and testing on large, complex systems. See Supplementary Figure B.2 for further details.
8

---

Under review as a conference paper at ICLR 2021
(a)
One step position MSE
Rollout position MSE
(b)
Trairing distribution nodes)
Forward Iterative
State after rollout steps
10-3
C-GNS-GD
Forward Iterative (a-0.25)
C-GNS-GD la=0.25)
Itrraionn
10
T-30
Generalization distribution
10"2
(20 nodes)
10
T-60
10
10
10a
T-90
0
10
15
15
Beertiion (M 5)
Figure 7: Generalization to more solver iterations and larger ROPE systems at test time. (a) Top
row: test accuracy for ropes with the same lengths as those during training (5-10 nodes). Bottom row: test accuracy for larger ropes (20 nodes) than during training. Left column: 1-step MSE. Right column: full 160-step rollout MSE. The x-axes indicate the number of solver iterations at test time (training used 5 iterations). The y-axis represents MSE values. Different line colors represent different models, as indicated by the legend. (b) Example of the rollouts from C-GNS-GD with different number of solver iterations.
DISCUSSION
We presented a general-purpose framework for constraint-based learned simulation, where learned constraint function implicitly represents the dynamics, and future predictions are generated via constraint solver. We implemented our framework using GNNs as the constraint function and gra- dient descent as the constraint solver, and tested it in variety of challenging physical simulation problems. Our results showed that our C-GNS has competitive or better performance compared to previous learned simulators. We demonstrated unique abilities to generalize to novel, hand-designed constraints, and use more solver iterations than experienced during training to improve the accuracy on larger systems.
We can hypothesize about the relationship between explicit, forward learned simulators and implicit, constraint-bast ones in terms of the sharing schemes of these architectures. The C-GNS has stronger inductive bias than the Forward GNN. The transformation of fc in C-GNS effectively ties parameters in the resulting Vy Sc function, and the solver iterations are analogous to how a recurrent neural network's s parameters are shared over iterations. In contrast, the message-passing steps in the Forward GNN used in our work are unshared. principle, the fp of the Forward GNN more expressive because if given enough depth, after training it could learn take parameter values that are equivalent to the shared parameters of C-GNS. Our results shown Figure supports this possibility: the Forward GNN with many more message-passing steps eventually approaches the C-GNS's performance. Moreover, we speculate the C-GNS's inductive biases contribute to its advantages in terms of incorporating novel hand-designed constraints and generalizing more solver iterations and larger systems.
More broadly, the performance, generality and unique advantages of constraint-based learned sim- ulation make it an important new direction in the advancement of machine learning methods for complex simulation problems in science and engineering.
REPRODUCIBILITY STATEMENT
We are committed to open-source the model code after the paper is accepted. Also, we are going to open-source the MuJoCo datasets that generated for this paper. We provide more details on the model implementation as well as the hyperparameters used for each model in the Supplementary Material.

---

Under review as a conference paper at ICLR 2021
REFERENCES
Brandon Amos and Zico Kolter. Optnet: Differentiable optimization as a layer in neural networks. In International Conference on Machine Learning, pp. 136-145. PMLR, 2017.
Jimmy Lei Ba, Jamie Ryan Kiros, and Geoffrey E. Hinton. Layer normalization, 2016.
Shaojie Bai, Zico Kolter, and Vladlen Koltun. Deep equilibrium models. arXiv preprini arXiv:1909.01377, 2019.
Shaojie Bai, Vladlen Koltun, and J Zico Kolter. Multiscale deep equilibrium models, arXiv preprini arXiv:2006.08656, 2020.
David Baraff. Fast contact force computation for nonpenetrating rigid bodies. In Proceedings the 21st annual conference on Computer graphics and interactive techniquex, pp. 23-34, 1994.
Sergey Bartumov, Vinod Nair, Peter Battaglia, and Tim Lillicrap. Continuous latent search for com- binatorial optimization. In Learning Meets Combinatorial Algorithms at NeurIPS2020, 2020.
P. Battaglia, Razvan Pascanu, Matthew Lai, Danilo Jimenez Rezende, and K. Kavukcuoglu. Inter- action networks for learning about objects, relations and physics. ArXiv, abs/1612.00222, 2016.
Peter W Battaglia, Jessica B Hamrick, Victor Bapst, Alvaro Sanzaez- Vinicius Zambaldi, Mateusz Malinowski, Andrea Tacchetti, David Raposo, Adam Santoro, Ryan Faulkner, et al. Relational inductive biases, deep learning, and graph networks. arXiv preprint 012661 2018.
Sofien Bouaziz, Sebastian Martin, Tiantian Liu, Ladislav Kavan, and Mark Pauly. Projective dy- namics: Fusing constraint projections for fast simulation. ACM transactions on graphics (TOG), 33(4):1-11, 2014
Michael M Bronstein, Joan Bruna, Yann LeCun, Arthur Szlam, and Pierre Vandergheynst. Geomet- ric deep learning: going beyond euclidean data. IEEE Signal Processing Magazine, 34(4):18 42, 2017.
Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud. Neural ordinary differ- ential equations. arXiv preprint arXiv:1806.07366 2018
Zhengdao Chen, Jianyu Zhang, Martin Arjovsky, and Léon Bottou. Symplectic recurrent neural networks. arXiv preprint arXiv: 1909. J 13334, 2019.
Alexandre Joel Chorin. The numerical solution of the navier-stokes equations for an incompressible fluid. Bulletin of the American Mathematical Society, 73(6):928-933) 1967.
Miles Cranmer, Sam Greydanus, Stephan Hoyer, Peter Battaglia, David Spergel, and Shirley Ho. Lagrangian neural networks, arXiv preprint arXiv:2003.0463 2020.
Filipe De Avila Belbute-Peres, Thomas Economon, and Zico Kolter. Combining differentiable PDE solvers and graph neural networks for fluid flow prediction. In Hal Daumé III and Aarti Singh (eds.). Proceedings of the 37th International Conference on Machine Learning. volume 119 of Proceedings of Machine Learning Research, pp. 2402-2411. PMLR, 13-18 Jul 2020.
Priya L Donti, David Rolnick, and J Zico Kolter. De3: A learning method for optimization with hard constraints. arXiv preprint 2rXi:22144.12225 2021.
David Duvenaud, Zico Kolter, and Matt Johnson. Deep implicit layers neural odes, deep equilib- -tuttrialll...... org/.
rium models, and beyond, 2020. URL httt::///mplicit-aaays.
Mare Finzi, Ke Alexander Wang, and Andrew Gordon Wilson. Simplifying hamiltonian and la- grangian neural networks via explicit constraints. arXiv preprint arXiv:2010. 135881 2020.
Rony Goldenthal, David Harmon, Raanan Fattal, Michel Bercovier, and Eitan Grinspun. Efficient simulation of inextensible cloth. In ACM SIGGRAPH 2007 papers, pp. 49-es. 2007.
10

---

Under review as a conference paper at ICLR 2021
Samuel Greydanus, Misko Dzamba, and Jason Yosinski. Hamitonnia neural networks. Advances in Neural Information Processing Systems, 32:15379- 15389, 2019.
George Em Karniadakis, loannis Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang. Physics-informed machine learning. Nature Reviews Physics, 3(6):422-440, 2021.
Dmitrii Kochkov, Jamie A Smith, Ayya Alieva, Qing Wang, Michael P Brenner, and Stephan Hoyer. Machine learning-acee accelerated computational fluid dynamics. Proceedings of the National Academy of Sciences, 118(21). 2021.
Yunzhu Li, Jiajun Wu, Russ Tedrake, Joshua B Tenenbaum, and Antonio Torralba. Learning particle dynamics for manipulating rigid bodies, deformable objects, and fluids. In ICLR, 2019.
João Loula, Kelsey Allen, Tom Silver, and Josh Tenenbaum. Learning constraint based planning models from demonstrations. In 2020 TEEF/RSJ International Conference on Intellitgent Robots and Systems (IROS), pp. 5410-5416. IEEE, 2020.
Michael Lutter, Christian Ritter, and Jan Peters. Deep lagrangian networks: Using physics as model prior for deep learning, arXiv preprint 1eXXiv:907.0440 2019.
Miles Macklin, Matthias Müller, Nuttapong Chentanez, and Tae- Yong Kim. Unified particle physics for real-time applications. ACM Transactions on Graphics (TOG), 33(4):1-12, 2014.
Damian Mrowca, Chengxu Zhuang, Elias Wang, Nick Haber, Li Fei-Fei, Joshua B Tenenbaum, and Daniel LK Yamins. Flexible neural representation for physics prediction. arXiv preprint 1eXiv:880600844 2018.
Matthias Müller, Bruno Heidelherger, Mareus Hennix, and John Ratcliff. Position based dynamics. Journal of Visual Communication and Image Representation, 18(2):109-118; 2007.
Tobias Pfaff, Meire Fortunato, Alvaro Gonzalez, and Peter Battaglia. Learning mesh-based simulation with graph networks. In International Conference on Learning Representations, 2021.
URL https:ttp//..c....t. /twp//ppeeeveeeeeeeeeeeeeeeeeeee net/ forum? 2 idd=oNgy1.00_P XP.
Yulia Rubanova, Ricky TQ Chen, and David Duvenaud. Latent odes for irregularly-sampie time series. In Proceedings of the 33rd International Conference on Neural Information Processing Systems, pp. 5320-5330 2019.
Alvaro Sanchez- Sanchez-Gonzales Gonzalez, Nicolas Heess, Jost Tobias Springenberg, Josh Merel, Martin Riedmiller, Raia Hadsell, and Peter Battaglia. Graph networks learnable physics engines for inference and control. In International Conference on Machine Learning, pp. 4470-4479, PMLR, 2018.
Alvaro Sanchez-Gonzale Victor Bapst, Kyle Cranmer, and Peter Battaglia. Hamiltonian graph networks with ode integrators. arXiv preprint arXiv:1909. 12790, 2019.
Alvaro Sanchez- Gonzalez, Jonathan Godwin, Tobias Pfaff, Rex Ying, Jure Leskovee, and Peter Battaglia. Learning simulate complex physics with graph networks. In Hal Daumé III and Aarti Singh (eds.), Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pp. 8468. PMLR, 13-18 Jul 2020. URL html.
https: http:///coceeeeeegeeee mllr p.rr.presss vvvvvv11111 httes//vv19//achheeee wwwanche--gnzaaeeeee
Bernhard Thomaszewski, Simon Pabst, and Wolfgang Strasser. Continuum- based strain limiting. In Computer Graphics Forum, volume 28, pp. 569-576. Wiley Online Library, 2009.
Nils Thuerey, Konstantin Weißennww Lukas Prantl, and Xiangyu Hu. Deep learning methods for aveerae-a navier-stokes simulations of airfoil flows. AIAA Journal, 58(1):25-36, 2020.
Emanuel Todorov, Tom Erez, and Yuval Tassa. Mujoco: A physies engine for model-based control. In 2012 IEFF/RSJ International Conference on Intelligent Robots and Systems, pp. 5026-5033. IEEE, 2012.
Po- Wei Wang. Priya Donti, Bryan Wilder, and Zico Kolter. Sameet: Bridging deep learning and log- ical reasoning using a differentiable satisfiability solver. In International Conference on Machine Learning, pp. 6545-6554. PMLR, 2019.
11

---

Under review as a conference paper at ICLR 2021
Rui Wang, Karthik Kashinath, Mustafa Mustafa, Adrian Albert, and Rose Yu. Towards physics- informed deep learning for turbulent flow prediction. In Proceedings of the 26th ACM SIGKDD International Conference Knowledge Discovery & Data Mining, pp. 1457-1466, 2020.
Jin-Long Wu, Heng Xiao, and Bric Paterson. Physics-informed machine learning approach for augmenting turbulence models: A comprehensive framework. Physical Review Fluids, 3(7): 074602, 2018.
Shuqi Yang. Xingzhe He, and Bo Zhu. Learning physical constraints with neural projec- tions. In H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, and H. Lin (eds.), Ad- vances in Neural Information Processing Systems, volume 33, pp. 5178-5189. Curran Asso- ciates, Inc., 2020. URL https: http://///.pc.eessssss.s.ss wtwprrpp.cc/ppppe/2200000 370057768118003666661118----2223 pdf.
Yao Zhang, Woong Je Sung, and Dimitri N Mavris. Application of convolutional neural network to predict airfoil lift coefficient. In 2018 AIAA/ACCE/AI/S/ACC Structures, Structural Dynamics, and Materials Conference, pp. 1903, 2018.
12

---

Under review as a conference paper at ICLR 2021
SUPPLEMENTARY MATERIAL
A IMPLEMENTATION
A.1 THE DATASETS
We generate the ROPE, BOUNCING BALLS and BOUNCING RIGINS datasets using MuJoCo phys- ical simulator, using a MuloCo timestep of 0.001 and taking every 30th time point. These datasets contain 8000/100/100 train/validation/test trajectories of 160 time points each.
The rope is attached on one end and moves under the force of gravity in 2D space. We randomly sample the number of nodes in [5, 10] interval and the length of the rope links in [0.6, 1.1] interval. The length of the rope links remains constant during the simulation.
ROPE
BOUNCING BALLS The simulation of the balls bouncing inside the box in 2D. The number of bouncing balls is randomly sampled in [5, 10J. The radius of the ball is randomly sampled in [0.11, 0.33. The size of the box is fixed to 5x5 in MuloCo coordinates.
BOUNCING RIGIDS A similar simulation to BOUNCING BALLS, where the nodes are connected into 2D rigid structure that bounces inside the box. We randomly sample the number of nodes between 3 and 6.
BOXBATH (Li et al., 2019) This is 3D particle-based dynamics using the FleX engine of the fluid enclosed inside a box, with a rigid cube floating on the surface of the fluid. Each simulation contains 960 fluid particles and 64 particles representing the cube. The dataset contains 2700/10/100 training/vaidaationnees trajectories with 150 time steps each.
We demonstrate the example of the rollouts for each environment in Supplementary Figure B. and Videos.
A.2 CONSTRUCTING THE INPUT GRAPH
As mentioned, we append history of three" recent velocities (see figure B.5(a-b) for ablations) for each node as node features. We also compute relative positions for each edge by subtracting the most recent positions in the history between the two particles adjacent to the edges, and use this as edge features.? Note do not include the absolute positions node features as the laws of physics are supposed to be invariant of the object's position in the space. We find that including the absolute positions into the nodes features harms the ability generalize larger environment, such longer rope.
Note that we do not provide the 'rest shape' true pairwise distances between the nodes for the ROPE or rigid structures in BOUNCING RIGIDS and BOXBATH. This makes the task harder as the model has to infer the rest shape from the input history, and may gradually deviated from the true shape.
Parameterizing the update Y For ROPE, BOUNCING BALLS, BOUNCING RIGIDS we use the average velocity as the update Y. For BOXBATH, the proposed update represents the normalized acceleration of the particle, initialized to zeros. We use normalized acceleration as the update Y, rather than raw velocity to better match the approach in Sanzaezz et al. (2020); Pfaff al. (2021). This proposed Y is is also concatenated as an extra node feature at each optimization iteration, and optimized to match the target future normalized acceleration. This acceleration is then un-normalized, and then Euler-integrated twice produce the next position. We chose to used normalized acceleration as the update Y, rather than raw velocity to better match the approach in Gonzalez et al. (2020); Pfaff et al. (2021).
We use five Box BATH to match model in Senchez-Gonzalde (2020)
For BOXBATH, we also provide the norm of the relative distances as additional vector to match features Ganzaez- al. (2020)
13

---

Under review as a conference paper at ICLR 2021
Handling walls Similar to Sanchez- Sanchez-Gonzalez al. (2020), include the euclidean distance between the center the node to each of walls as additional node features, clipping at a maximum value to avoid this to become a proxy for absolute position. For BOUNCING BALLS and BOUNCING RIGIDS we clip the distance at 2, and for BOXBATH at 0.08. For the iterative models, we update the distances to the walls after every step of constraint optimization.
Node type We also provide an additional one-hot node feature indicating the node type (e.g. rigid, fluid, fixed). For BOUNCING BALLS, for which the object size varies between the objects, we pro- vide the radius of the object an additional node feature. Other features, like the relative positions are computed between centers of the nodes, and the network has to account for the object size to detect collision.
Graph edges BOUNCING RIGIDS and BOUNCING BALLS we use a fully-connected graph. In ROPE we add edges between nodes that adjacent within the rope. In BOXBATH add edges between particles that are within a radius of 0.08 from within each other, and then recompute these edges at every step of rollout according the updated positions (as in Sanchez-Gonzalez et al. (2020)).
Normalization For the large scale datasets BOXBATH found it important to normalize inputs and targets to zero-mean unit-variance (as in Sanchez-Gonzalez et al. (2020)). In the other datasets, BOUNCING BALLS, BOUNCING RIGIDS, and ROPE, the scale of the features was al- ready close to zero-mean unit-variance, except for the input/target velocities in BOUNCING BALLS, BOUNCING RIGIDS which we applied a scaling factor 100.
Noise
To stabilize rollouts in BOXBATH, we added noise to the input sequences in the same man- ner and with the same magnitude as Sanchez-Gonzalez al. (2020).
Fixed particles Some of the datasets contain the fixed nodes that do change the their position in the simulation, like the first node in the ROPE. We prevent the update for those nodes by using stop-ge -radderai
A.3 MODEL IMPLEMENTATION
Computing the constraint gradients To compute the gradients of the constraint scalars for the batch of graphs, we use the vector-jacobian product (vjp) function. VJP does explicitly con- struct a jacobian, and its asymptotic computational cost is the same as the forward evaluation of the constraint function.
Constraint function We use the mean aggregation for the per-node outputs to obtain the scalar constaint value for the entire graph. For gradient descent, we additionally take a square of per- node outputs before aggregating them. We use a fixed learning rate of 0.001 for Gradient Descent constraint solver.
ConstraintMIIP For the ConstraintML.P Yang et al. (2020). We concatenate the features for each node into a single vector and run an ML. to produce scalar constraint output. The node features include the absolute positions, velocities and distances to the walls for each node for the optimized state and for the past three time points. Note that this model does not have access to the position differences between the nodes (this information is in the edge features for the graph network). Therefore we add the absolute positions to the node features, which simplifies the task for the model. Additionally, we adapt the ConstraintMLP to handle the scenes with the variable number of the nodes. To do we pad the missing nodes with zeros vectors up to the maximum number of nodes in the dataset.
similar setup in the original paper by
A.4 HYPERPARAMETERS
Rope We used message-paseee passing step for the graph- seachh-nttoor-bbaseeee models. The MLPs for node and edge processing consist have 3 hidden layers with 256 hidden units, and hidden node and edge latent size of 64. We use "softplus" activation and LayerNorm Ba et al. (2016).
14

---

Under review as a conference paper at ICLR 2021
For the ConstraintMLPP we performed the hyperparameter search and chose the best performing configuration. We use the MLP with 10 hidden layers and 2048 hidden units with 'softplus' activa- tion with no LayerNorm for the constraint function. We include the absolute positions into the node features, otherwise the ConstraintMI.P model does not contain any information about the relative positions of the nodes (that we add into the edges for the graph model).
Bouncing Balls We used message-passiinn step for the yraph-network-basse models. The MLPs for node and edge processing consist have 3 hidden layers with 256 hidden units, and hidden node and edge latent size of 64. We use 'softplus' activation and LayerNorm after every MLP, except the final decoder.
For ConstraintMILP used the best-performing configuration from hyperparamter search. We use the MLP with 10 hidden layers and 2048 hidden units with 'softplus' activation with no Layer- Norm.
Bouncing Rigids We used message-passing step for the graph-network-basee models. The MLPs for node and edge processing consist have 3 hidden layers with 256 hidden units, and hidden node and edge latent size of 64. We 'tanh' activation and LayerNorm after every MLP, except the final decoder.
For ConstraintMLPP we use the MLP with 10 hidden layers and 2048 hidden units with 'softplus' activation. On this dataset, we found that adding a LayerNorm to ConstraintMLP-FP model helps the stability of the model.
Box Bath We used I message-passing step for the graph-network-based models. All other hy- perparameters are as in Lunchez-Gonzale et al. (2020). The MLPs for node and edge processing consist have hidden layers with 128 hidden units, and hidden node and edge latent size of 128. We use "softplus" activation and LayerNorm after every MIP, except the final decoder.
Training We train the models for million steps on ROPE, BOUNCING BALLS and BOUNCING RIGIDS. We used the initial learning rate 0.0001, with the decay factor of 0.7 and learning schedule of (1e5, 2e5, 4e5, 8e5). We use batch size of 64. We trained for 2.5 million steps for the experiments studying the number of solver iterations. On Box BATH we trained for 2.5M steps with batch size of 2. With learning rate starting at 0.001 and decaying continuously at rate of 0.1 every I 1M steps, as in Gonzaez et (2020).
Loss on multiple iterations For the variable solver iterations results we experimented training models with a loss imposed not only the last iteration output (NV but at intermediate outputs (+) the iterative models. We used exponentially decaying relative weights (from last to first), such that the relative weight w; for the loss term for the P output was w; = a N-i The goal of a is to encourage the model to make progress towards the minimum at each iteration, without penalizing the model heavily for not reaching the solution < N iterations. More details on the choice of are provided Supp. Fig. B.2.
15

---

Under review as a conference paper at ICLR 2021
Rope (GT)
Rope
(Rollout)
Bouncing Balls (GT)
Bouncing Balls (Rollout)
Bouncing Rigids (GT)
Bouncing Rigids (Rollout)
BoxBath
(GT)
BoxBath (Rollout)
Figure B.I: Examples of the rollouts for our simulation environments,
16

---

Under review as a conference paper at ICLR 2021
B SUPPLEMENTARY PLOTS AND TABLES
sapp
Bullos postiion
arbiisary of(we)
12-1
Cid cD 12w0
duu di (o.4.
105
10P
C-SAS CD towd
12"4
107
107
Cenven Genersti weseeaaan bastioniio
10'
101
300-
12.
150
ham
Figure B.2: Generalization to more solver iterations at test time N as function of a. Imposing loss only at the last iteration (a=1), causes the model to
steps
steps
MP
Forward
Iterative GNN
Rollout MSE
10P
C-GNS-GD
Forward GNN
10-1
3
3 iterations
3 iterations
1
5 iterations
12344.5 7
8 910
iterations
iterations
Figure B.3: Comparison of C-GNS to the baselines with different number of message-passing layers and number of constraint solver iterations (Test MSE on the full rollout). Last facet: 1-step rollout error of the state-of-the-an Forward GNN with to 10 message message-passing steps for reference.
MP stops
2 MP stops
MP stops
4 MP steps
5 MP steps
10
MSE
10-4
10-5
2
5
2
5
2
5
2
5
2
5
iter ations
iterations
iterations
iterations
# iterations
10a
MP steps
2 MP steps
MP steps
steps
5 steps
100
Roilout
10-1
10-2
2
5
2
5
2
5
2
5
2
5
iterations
iterations
iterations
Reerttions
iterations
Figure B.4: Test 1-step MSE and Full Rollout MSE of the C-GNS with Neural Projection with different number of message-passing layers and number of constraint solver iterations on the Rope dataset. Top row: position MSE on the full rollout.
17

---

Under review as a conference paper at ICLR 2021
Table B.1: Median performance of the models on different datasets. *ConstraintMII. results shown for the 'BoxBath' dataset as it is infeasible to run MLP constraint on 1000s of nodes. The standard deviation from the median shown over 5 random seeds. We do not show the results for the models where the median value, or standard deviation is more than 1000 times larger than the best model in each column. Results for all seeds are shown in figure B.6 Note that the tables use different scales to demonstrate the errors on 1-step enror, 10-step rollouts and full rollouts.
Model
One-step position MSE
Rope (le-4)
Bouncing
Balls
Bouncing
Rigids
Box Bath (le-7)
(1e-6)
(le-7)
ConstraintML.PPOO
0.740 0.021
7.716 1.610
8.017 = 1.029
ConstrainMMLP
32.332 + 90,408
31.336 + 0.000
25.221 0.000
Forward GNN
0,406 = 0.005
0.126 = 0.009
1.152 = 0.051
1.440 = 0.038
lterative GNN
0.135 0.001
0.168 0.019
1.172 = 0.298
1.042 0.103
C-GNS-FP
0.698 0.297
5.896 5.930
0.954 1.885
C-GNS-GD
0.110 0.005
0.103 1 0.032
0.884 + 0.163
0.998 1 0.038
Model
Rollout position MSE (10 steps)
Rope (Ie-3)
Bouncing (1c-5)
Balls
Bouncing
Rigids
Box Bath (le-5)
(Ic-4)
MontrrinLMLP-GDD
1.463 0.231
0.998 + 0.436
ConstraintMLP-FF
692.940 + 0.000
4.460 = 0.000
Forward GNN
27.585 = 7.335
2.853 0.277
1.217 = 0.346
0.386 0.061
Iterative GNN
0.251 0.145
2.260 = 0.429
0.238 + 0.026
0.174 0.020
C-GNS-FP
1.458 2.323
661.72 + 372.49
0.200 = 0.321
C-GNS-GD
0.226 + 0.058
0.613 + 0.333
0.248 + 0.038
0.288 0.055
Model
Rollout position MSE
Rope (le-1)
Bouncing Balls
Bouncing
Rigids
Box Bath (le-2)
(le-1)
ConstraintMLP-CG
20.619 + 94.212
2.161 0.217
4.241 = 1.970
ConstrainMLLPII
6.664 1 0.000
13.128 0.000
Forward GNN
18.305 + 2.340
0.389 = 0.022
1.174 = 0.486
0.756 - 0.089
Iterative GNN
0.546 + 0.026
0.445 = 0.044
0.306 + 0.253
0.609 + 0.015
C-GNS-FP
2.804 3.432
0.689 0.099
C-GNS-GD
0.602 0.088
0.308 1 0.142
0.374 + 0.171
0.654 1 0.002
18

---

Under review as a conference paper at ICLR 2021
Rope
Bouncing Balls
Bouncing Rigids
Box Bath
10-1
4 10-7
10-2
98 BR+a
31e-07
Onc-step posttin MSE
10
10-3
10-3
10:
10-7
10-4
10-7
10-7
10
10"5
2072
10a
80-38'1
109
10-1
6 10
Rollout pocition
20
10-1
10-2
10-0
1004
(10
10-2
10
10-1
3 10-6
10-3
10
XBXX
xx 10
10
10
10"
204
1011
9 10-3
Rollout positon MSE
C-MLP.CC GD
10a
10a
107
8 10
C-MLP-FP
10a
Forward GNN
Iterative
10!
109
10a
x xx-00-
C-GNS-GD
109
10-1
20
10-1
6 10-2
Figure B.6: Test MSE on the node positions predictions across different models. First row: 1- step position MSE. Second row: position MSE eтror for the first 10 steps of a rollout. Third row: position MSE error full 160-step rollout MSE. The bar height represents the median MSEs over random seeds. The x marks show the MSF metric for each random seed. The black arrows indicate the MSE metric for random seed exceeding the upper y limit of the figure.
(a) # past time point
(b) # past time point the history: constraint MSE
(c) Squaring per-node outputs: position MSE
(d) Squaring per-node
in the history: position
outputs: constraint
MSE
MSE
Figure B.5: Ablations of the modelling choices.
19


Reviews:
[{"review_id": "mB78WNl5yMA", "reviewer": "Reviewer_TMn2", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "1. I find the notations X and Y with and without hats very confusing from the beginning. It would be great if the paper could clarify the semantic meaning of these notations.\n* Fig. 1 (a): its caption says the UPDATER uses $\\hat{Y}$ to update $X_t$, but the figure shows an UPDATER using $Y$ to update $\\hat{X}\\_t$. The notations in Fig. 1 (a) itself is extremely confusing: it starts with $X_1$ (without a hat), proceeds to $\\hat{X}\\_t$ (with a hat), but then uses $X\\_{\\leq t}$ (without a hat), and finally predicts $\\hat{X}_{t+1}$ (with a hat). How did you determine when to use a hat and when not?\n* When defining the rollout, $X$ with and without hats are mixed. From the rollout’s definition, it starts with $X_t$ but is followed by $\\hat{X}\\_{t+1}$, $\\hat{X}\\_{t+2}$, …. . Does this mean you manually choose t and separate the whole rollout into two parts: before $t$ every $X$ is without a hat and after $t$ every $X$ is with a hat (generated by the simulator)?\n\n2. I do not see a strong reason for introducing the explicit iterative simulator for comparison. It looks like $f\\_{DI}$ can be rewritten as $f_D$ plus a few linear operators: $f\\_{DI}(X_{\\leq t}) = (f_D(X\\_{\\leq t}) - X_t) / N$. Therefore, I would expect $f_D$ and $f\\_{DI}$ to have very similar capabilities. What insight would we expect to get from comparing them?\n\n3. It looks like the proposed constraint-based predictor runs a fixed number of gradient-based iterations and a fixed step size without guaranteeing the constraint is satisfied eventually. Therefore, I feel it is a bit too much to claim that the proposed approach is doing constraint-based simulation and to claim that the proposed network is a constraint solver. I think it would be necessary to either tone down this claim in the title, abstract, introduction, etc., or state explicitly at the beginning of the paper that the proposed approach does not guarantee the constraints are always satisfied.\n\n4. The network design ensures translation-invariance by taking as inputs position offsets instead of absolute positions. Does the network also ensure rotation-invariance?\n\n5. “These node-wise constraint values are averaged to compute a single scalar c constraint for the entire system”. This seems questionable to me if the goal is to solve constraints $f_C(X\\_{\\leq t},  \\hat{Y}) = 0$, in which case I would expect that using an average of $|c|$ or $c^2$ makes more sense.\n\n6. Looking at the experimental task domains, I feel there are two types of constraints involved in the simulator: equality constraints (e.g., end points of the consecutive rope segments must share the same location) and inequality constraints (e.g., bouncing balls must have nonnegative distances to the boundary). Does this paper deal with both equality and inequality constraints?\n\n7. Both N and lambda seem to be crucial hyperparameters that need to be chosen for different environments individually. If this approach is applied to a new environment, how would you determine a proper N and lambda?\n\n8. “In principle, iterative and constraint-based simulators should find more accurate solutions by increasing the number of solver iterations, N.” I am not sure I fully agree with this claim because the iterative solver is not equipped with a line search algorithm that adaptively changes the step size.\n\n9. The large rollout MSE (orders of magnitude larger than #timesteps x One-step MSE) in Fig. 3 bottom seems to imply the learned simulator is not a good replacement of a numeric simulator because it accumulates errors from all time steps. I understand this may be a common issue that many other neural-network-based simulators also suffer from (all baselines in Fig. 3 have accumulated substantial errors and produced large rollout MSE), but I am still wondering whether the authors could give people a strong reason why it is useful to develop such a neural-network-based simulator if it is not accurate.\n\n10. Similarly, I wonder if this paper could provide more discussions on the generalizability of the learned network on the static properties of the environment (the Z vector in the main paper), e.g., density, material types, time step size, etc. My understanding is that the network needs to be retrained if Z is updated, or the training set needs to be augmented to see various Z values. This does not seem to be very ideal for a simulator. Again, I understand this may be a common problem for many neural-network-based simulation papers, so I won’t hold it against this paper too much.\n", "summary_of_the_paper": "This paper presents a neural-network simulator that learns to solve constraints inspired by classic physics-based simulation. The proposed simulator uses a graph neural network to encode a constraint solver and shows results in a number of simulation environments.\n\nThe main contribution of the paper is its idea of using graph networks to encode constraint-based simulation.\n", "summary_of_the_review": "I recommend rejection based on the current status of the paper. My concerns are listed above.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "d9lyggIg0Eq", "reviewer": "Reviewer_kLRf", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "[Strength]\n\nThis paper tackles an important question of how we can incorporate prior knowledge in the form of explicit physical constraints in the learning-based simulators to enable better generalization.\n\nThe experiments on four environments show that the proposed method can deliver better prediction results than unconstrained baselines.\n\n\n[Weakness]\n\nWhile I like the direction this paper is going, I have concerns regarding the expressiveness of the experiments and the missing details of the method.\n\nPrior methods on learning-based physics simulators have shown results on a set of much larger-scale and more complex environments involving fluids on rough terrain, fabrics with novel geometries, etc. [1, 2]. The experiment environments used in this paper may be a bit too simple compared to what's out there in the literature, making it hard to know how the method works in larger and more complicated scenarios.\n\nContinuing from my previous point, [1] showed generalization to environments with drastically different geometry than seen during training, and [2] showed that the model could scale up to significantly larger and more complex cloth than seen during training. Adding constraints based on our understanding of physics is supposed to improve the model's generalization ability. As a result, I don't think the current experiments in the paper are enough to demonstrate the benefit of the constrained optimization process. The authors should consider including concrete experimental evidence on how the incorporated constraints may lead to even better/larger-scale generalization than what's already shown in the literature.\n\nThe authors should also consider including more details on how they construct the constraint function, f_C, e.g., for the fluid, rigid object, boundary conditions, etc. Without further details, it is hard for me to imagine how they are defined and implemented.\n\nRelated to my previous point, do we need any assumptions on f_C other than being differentiable? For example, for discontinuous events like contacts, I imagine we can differentiate through the LCP constraints, but how useful are the gradients?\n\nHow much more computing resources and time are needed to apply gradient descent on the constraint function? Multi-step message passing and solver iterations do not come for free. It may significantly increase forward prediction time for each time step. Therefore, it is essential to provide the time spent on each forward pass for different design choices and discuss the trade-off between the performance gain and the decrease in computational efficiency.\n\nAccording to Figure 3, it is a bit hard to know which of the following three methods works better: (i) Iterative GNN, (ii) C-GNS-GD, and (iii) C-GNS-FP. When giving a new scenario, is there a way to know which one we should use, or should we try all of them and choose the one that works the best?\n\nHow did you choose the weight on the constraint term when incorporating novel constraints at test time? From the video, the rope is jittering. What might be the reason, and is it possible to resolve it? [1] shows generalization results of fluid simulation on unseen terrains much different from what the model was trained on. [3] also showed examples of generalization to unseen obstacle configurations. How would you compare the way your paper and [1, 3] incorporate new boundary constraints during testing?\n\n\n[1] Benjamin Ummenhofer, Lukas Prantl, Nils Thuerey, Vladlen Koltun, \"Lagrangian Fluid Simulation with Continuous Convolutions\"\n[2] Tobias Pfaff, Meire Fortunato, Alvaro Sanchez-Gonzalez, Peter W. Battaglia, \"Learning Mesh-Based Simulation with Graph Networks\"\n[3] Alvaro Sanchez-Gonzalez, Jonathan Godwin, Tobias Pfaff, Rex Ying, Jure Leskovec, Peter W. Battaglia, \"Learning to Simulate Complex Physics with Graph Networks\"\n\n\n===================\n\n[Post Rebuttal]\n\nI thank the authors for the detailed feedback, which addressed most of my concerns. I hope the authors can incorporate the response into the manuscript to improve its clarity. I'm happy to raise my score from 5 to 6.", "summary_of_the_paper": "This paper aims to add explicit/human-defined constraints to learning-based simulation frameworks, where a learned constraint function implicitly regularizes the dynamics, and future predictions are generated via a constraint solver. The authors built the framework on top of graph neural networks (GNNs) to capture the compositionality of the underlying system and enforce the constraint using an implicit constraint function optimized via gradient descent.\n\nThe authors tested the proposed method in four physical simulation environments, including rope, bouncing balls, bouncing rigids, and BoxBath. Experiment results show that the proposed C-GNS has a competitive or better performance compared to prior learned simulators. In addition, they have also demonstrated that C-GNS can generalize to unseen, hand-designed constraints by applying more solver iterations than experienced during training to improve the accuracy on larger systems.", "summary_of_the_review": "While I like the direction this paper is going, I have concerns regarding the expressiveness of the experiments compared to previous work on learning-based physics simulation, e.g., extreme extrapolation generalization on systems much larger than what the model was trained on.\n\nThere are also a lot of missing implementation details, making it hard to know (1) how the authors define and solve the constraints, e.g., for fluids, rigid objects, boundary conditions, etc., (2) how hard the constrained optimization problem really is, and (3) how does it affect the final performance.\n\nAs a result, I currently lean towards the rejection side.\n", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "cRGVUfas0Jd", "reviewer": "Reviewer_AXAA", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Strengths:\n\n1.\tThis paper is overall well written. The authors have clearly demonstrated the pipeline of the proposed strategy.\n\n2.\tAlthough the novelty of combining GNN with constraint projection is weak (see the weaknesses below), it is valuable to check if this method can outperform those typical forward approaches (such as the work by Sanchez-Gonzalez et al., 2020). The experimental evaluations generally serve this purpose. \n\n\nWeaknesses:\n\n1.\tThe biggest concern is that the novelty is weak. At its core, this paper applies the pipeline in (Sanchez-Gonzalez et al., 2020), including the strategy by first computing the predictor and then updater, and the usage of GNN for interaction modeling. The main difference is that for the predictor, it replaces the traditional forward prediction with the iterative gradient-based solver of the constraint approximator, which is interesting. Yet, this idea has been proposed by Yang et al. (2020) in terms of the iterative projections along the gradient direction of the constraint network. Although the authors have further augmented the input of the constraint function with Y to take the dynamics into account, this modification seems minor and a straightforward enhancement.  \n\n2.\tRegarding the constraint satisfaction. The authors first predict the changes in position (Y) and then update the next state X_{t+1} via an Euler integrator. Even the prediction of Y is derived via the constraint solver, the constraint will be broken after the following updater from X_{t} to X_{t+1}, which is problematic to maintain the hard constraints (such as the case in Figure 6 (c) where the movement within walls and floors is forbidden). How does the proposed method tackle this issue? In the work by Yang et al. (2020), the authors use an opposite order by first updating X_{t} and then projecting the positions onto the constraint manifold, which is able to meet any kind of constraints.\n\n3.\tFor the comments above, there are several important baselines that are not tested in the experiments: 1) using the method by Yang et al. (2020) but with the GNN projector; 2) first updating X_{t} and then predicting Y with other setting unchanged in the current framework; 3) augmenting the explicit simulators (both the iterative and non-iterative versions) with a regulation loss to enforce certain hard-crafted constraints such as the cases in Figure 6 (c). \n\nOther comments:\n\n1.\tThis paper is almost well organized, but there are still some confusions in the current version. \n1.1 In introduction, the authors mention that both families of simulators (explicit forward vs implicit constraint-based). Is this statement discussed in previous papers? Is there any citation?\n1.2 In related work, the authors state that the neural projection by Yang et al. (2020) is the first work that uses learned constraints. However, they also introduce that “Recent methods have been proposed for learning constraint functions and solving them in a model’s forward pass” such as “Deep Implicit Layers” and “Deep Equilibrium Models”. These two statements seem self-contradictory.\n1.3 In section 4.3, the authors claim that in Figure 3 the C-GNS-GD’s performance was generally better than the other model variants, which is not true. In terms of the rollout MSE, Iterative GNN outperforms C-GNS-GD in three out of four cases. The authors are suggested to provide more explanations here.\n\n2.\tIt is good that the proposed method is translation-invariance. Yet, besides this symmetry, there are other cases, such as rotation invariance/equivariance. This is important for improving the generalization ability of the simulator, given the fact that if we rotate the input states under a certain angle, the output changes in the same way. Have the authors taken this symmetry into account? \n", "summary_of_the_paper": "This paper presents to simulate physics via a constraint-based approach instead of direct prediction. In particular, the authors first employ GNN taking as input the history positions and dynamics, to capture the interaction between different particles within the system, whose output is considered as the constraint satisfaction scalar. Then, the gradient with respect to the constraint function is applied as the update of the dynamics over a certain number of solver iterations. The experiments are conducted on a variety of challenging physical domains, including simulated ropes, bouncing balls, colliding irregular shapes and splashing fluids.", "summary_of_the_review": "Overall, given the limited novelty and insufficient experimental evaluations, I initially suggest weak rejection.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "CmsQ9GVZuxV", "reviewer": "Reviewer_WhyX", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "This manuscript extends an idea from [1] that the learning of physical simulation can be viewed as a constrained optimization problem. There are a few interesting points:\n\n1. The authors extend it using graph neural networks. This modification generalizes this method to the indefinite number of states in the physical system.\n2. The authors claim that this method trade-off between inference time and accuracy dynamically so that in the test time they are possible to increase the number of steps for better convergence.\n\nDespite these respectful contributions, I still have a few questions and reservations:\n\n1. I notice different activation functions are applied to different environments. Are they intentional? Is there a guideline for choosing?\n2. LayerNorm will be affected by the magnitude of the static information. For example, if you simulate under the continuous mechanism, Young's Modulus can be so huge that kills other entries in the vector. I wonder how should the user deal with it.\n3. How does this method generalize to the static information and time step?\n4. Second-order optimization usually provides super-linear convergence when it is close to the optimum, which is a similar case referring to figure 5. I will find a comparison between the proposed first-order method and Newton's or quasi-Newton's method helpful.\n5. The authors claim that adding an $\\alpha$ improves the generalizability to unseen parameters including articulation length and optimization steps. Is there a reason why it is not a part of the standard model given it has an advantage?\n6. 3 and 5 are separately used in different environments. Some discussion on the choice will be helpful (I understand 5 is from previous work). Is it because of the time integration method used in Mujoco?\n7. My major reservation is the technical novelty compared to its preceding work [1]. I find the extend to graph neural networks exciting, but in the meantime, incremental. I wonder if the authors can give me more explanation on the delta between them, maybe supported with a more impressive application.\n\n[1] Yang, Shuqi, Xingzhe He, and Bo Zhu. \"Learning Physical Constraints with Neural Projections.\" Advances in Neural Information Processing Systems 33 (2020): 5178-5189.", "summary_of_the_paper": "This manuscript proposes to learn the numerical solutions to Lagrangian physical simulation with a constraint-based inference method on graph neural networks. This method involves an iterative update during inference and thus enables test-time dynamical correction. The contributions are: (1) this manuscript builds a scalar predictor to indicate how well the constraint is agreed. (2) this manuscript proposes using the graph neural networks as the backbone to deal with a variable length of the physical domain. They also examine the effectiveness with a bunch of experiments including the state prediction experiment on four different environments and multiple ablation studies towards different hyper-parameters.", "summary_of_the_review": "I appreciate the idea of using graph neural networks with a constraint-based learning paradigm. However, there are also a few technical questions remaining unsolved. I personally find this work is not satisfying me with solid novelty, but I am open to change if the authors address my concerns.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}]

Author rebuttal:

Meta-review: The paper uses graph-based neural networks to ensure constraint-based simulation.  Even though the approach is a good one, it is only incremental w.r.t. the work published by Yang et al at NeurIPS in 2020; then, the experimental section is not convincing enough.

While the authors indicate their dissatisfaction with one of the reviewers' assessment, the overall reviews of the paper are not very positive.
Decision: Reject