CROM: Continuous Reduced-Order Modeling of PDEs Using Implicit Neural Representations
Authors: Peter Yichen Chen
Jinxu Xiang
Dong Heon Cho
Yue Chang
G A Pershing
Henrique Teles Maia
Maurizio M Chiaramonte
Kevin Thomas Carlberg
Eitan Grinspun
Venue: ICLR 2023 notable top 25% (2023)
Keywords: PDE
implicit neural representation
neural field
latent space traversal
reduced-order modeling
numerical methods

Abstract: The long runtime of high-fidelity partial differential equation (PDE) solvers makes them unsuitable for time-critical applications. We propose to accelerate PDE solvers using reduced-order modeling (ROM). Whereas prior ROM approaches reduce the dimensionality of discretized vector fields, our continuous reduced-order modeling (CROM) approach builds a low-dimensional embedding of the continuous vector fields themselves, not their discretization. We represent this reduced manifold using continuously differentiable neural fields, which may train on any and all available numerical solutions of the continuous system, even when they are obtained using diverse methods or discretizations. We validate our approach on an extensive range of PDEs with training data from voxel grids, meshes, and point clouds. Compared to prior discretization-dependent ROM methods, such as linear subspace proper orthogonal decomposition (POD) and nonlinear manifold neural-network-based autoencoders, CROM features higher accuracy, lower memory consumption, dynamically adaptive resolutions, and applicability to any discretization. For equal latent space dimension, CROM exhibits 79$\times$ and 49$\times$ better accuracy, and 39$\times$ and 132$\times$ smaller memory footprint, than POD and autoencoder methods, respectively. Experiments demonstrate 109$\times$ and 89$\times$ wall-clock speedups over unreduced models on CPUs and GPUs, respectively. Videos and codes are available on the project page: https://crom-pde.github.io

Published as a conference paper at ICLR 2023
CROM: CONTINUOUS REDUCED-ORDER MODELING OF PDES USING IMPLICIT NEURAL REPRESENTATIONS
Peter Yichen Chen
Jinxu Xiang'
Dong Heon Choo
Yue Chang
G A Pershing'
Henrique Teles Maia!
Maurizio M. Chiaramonte
Kevin Carlberg
Eitan Grinspun*1
"Columbia University
Meta Reality Labs Research
MIT CSAIL
UUniversity of Toronto
ABSTRACT
The long runtime of high-fideligg partial differential equation (PDE) solvers makes them unsuitable for time-critical applications. We propose to accelerate PDE solvers using reduced-order modeling (ROM). Whereas prior ROM approaches reduce the dimensionality of discretized vector fields, our continuous reduced-order modeling (CROM) approach builds a nwwdimensioonll embedding of the contin- nous vector fields themselves, not their discretization. We represent this reduced manifold using continuously differentiable neural fields, which may train any and all available numerical solutions the continuous system, even when they are obtained using diverse methods or discretizations. We validate approach extensive range of PDEs with training data from voxel grids, meshes, and point clouds. Compared prior discretizatioon-eeennde subspace proper orthogonal decomposition (POD) and nonlinear manifold neural- network-basee autoencoders, CROM features higher accuracy. lower memory con- sumption, dynamically adaptive resolutions, and applicability to any discretization. For equal latent space dimension, CROM exhibits 79x and 49 better accuracy, and 39x x and 132x smaller memory footprint, than POD and autoencoder methods, respectively. Experiments demonstrate 109x and 89 x wall-clock speedups over unreduced models on CPUs and GPUs, respectively. Videos and codes are available on the project page: https: crom ceerommmpee.... github. io.
ROM methods, such linear
INTRODUCTION
Many scientific and engineering models are posed as partial differential equations (PDEs) of the form
F V.f.vf.V'f
f(x. t) Ω Rd
(1)
subject to initial and boundary conditions. Here f is spatiotemporal dependent, multidimensional continuous vector field, such temperature, velocity, or displacement; V and (-) are the spatial and temporal gradients; 52 C Th' and T C IR are the spatial and temporal domains, respectively.
We may solve for f by discretizing in space, f(z, t) x f2, t) = Ee a' (t)) transforming the continuous spatial representation to a (P dd-dimensional vector whose coefficients a a'(t) T+ dd and the corresponding basis functions N'(x)) + R (c.g., polynomial basis, fourier basis) approximate the continuous solution. For instance, if N is the linear finite element basis, the coefficients a' (t) = f("', t) are field values at spatial samples x (Hughes) 2012).
After introducing temporal samples (t,, -T we temporally evolve the solution by solving for P unknowns (a (Ln+1)} given the previous state (a'(1,.)) Unfortunately, when P is large, processing and memory costs of these full-order solves become intractable. To alleviate this computational burden, prior model reduction techniques (Berkooz et al. 1993 Willeox & Peraire [2002] Benner et al.
2015) construct manifold caatfffldd-parmmmtriiiataii function gr R. Ira with < Pd, such that every low-dimensional latent space vector q(t) I² maps to discrete field gr(q) + (a ,...,, ?
For instance, for linear finite elements BBariič & James 2005), gp( (f .).... f(aP, 1))

---

Published as a conference paper at ICLR 2023
(a) Prior approach
Discretized vector field
(x),,1)
g+
gp
gr
((x.11)
(ta)
((ty)
((c) RR, <d
latent space trajectory
x Q 0
g
f(x, E IR
9
9
Continuous vector field
(b) Our approach
Figure I: Model reduction solves PDEs via temporal evolution of the fow-dimensional latent space vector (t).. (a) Prior work assumes that the low-dimensional representation gr is built for the already-discretizz
vector field; (b) our approach constructs the monifold-prrameterizatt function g directly for the continuous vector field itself. In this case, the vector field f represents the twisting material governed by the elastodynamics equation.
as depicted in Figure Da. ROM saves computation because it requires evolving only < Pd latent space variables
Since existing ROM approaches apply already-discreetzee fields, model training and PDE solving are tied to the dimension and discretization type of the training data, causing key limitations:
Discretization dependence. If alter the training simulation resolution (P) or the discretization types (c.g., meshes to point clouds), we must also alter the architecture and numbers of parameters.
Memory scaling. Memory footprint grows with discretization resolution P.
Fixed discretization. We cannot dynamically adapt spatial resolution P. discretization type, or basis function N during latent-space-Pii PDE solves, e.g., dynamic remeshing (Peraire al. 1987).
Altogether these problems arise because the architecture of gr(g) is tied to the discretization (a (a'....,a")) a...."") "²)"
Introducing a dlicretization-ineeeeede train a manifold- itself, not its discretization (see Figure (p). Note that the domain and co-domain of g are continuous domains: they do not depend on the choice of discretization(s) used at any stage of the process, i.e., during preparation of training data, nor during latent-space-PDE solving. In this sense, the manifold- parameterization architecture discretization independent. In our implementation, g is embodied as an implicit neural representation (Park et al. 2019 Chen & Zhang 2019 Mescheder et al. (2019), also known as a neural field, yielding a smooth and amalyticaly-ddiferee This representation's 's memory footprint depends the complexity of fields produced by the PDE, not the discretization resolution.
architecture In an alternative point of departure, we manifold-parammttriiaatte parameterization function g(x, q) f(x, to approximate the continuous field
manifol- paramcterization.
After training, we evolve the latent variables, as governed by the PDE. for previously-unexplored parameters. Unlike approaches that discard the PDE after training, we evaluate the original PDE at a small number of domain points every time integration step. We validate our approach classic PDEs with discretized data from voxels, meshes, and point clouds. In comparison to the full-order model, our approach reduces the number of spatial degrees of freedom, memory, and computational cost. In comparison to prior linear and nonlinear our method exhibits higher accuracy and consumes less memory. To highlight another benefit of being discretization discretization-annostiii, we demonstrate an elasticity simulation that readily adapts mesh resolution.
model reduction methods,
2 RELATED WORK
Reduced-Order Modeling for PDEs. Early works identifying low-dimenssonall latent space focused on linear methods (Berkooz et al. 1993. Holmes et al. 2012), e.g., proper orthogonal decomposition (POD) or principal component analysis (PCA). Recent nonlinear manifolds (Fulton) ct al. 2019 Lee & Carlberg 2020). often constructed via autoencoder neural networks, have been
'The latent space vector is also known the feature, subspace, state vector; or the generalized coordinates.
2

---

Published as a conference paper at ICLR 2023
shown to significantly outperform their linear counterparts slowly decaying Kolmogorov n-width problems (Peherstorfer, (2022. Most of these prior works exclusively focus on building a latent space for the diready-dicretized vector fields. Chen et al. (2023); Pan et al. (2023) and our work attempt to construct the latent space for the continuous vectors themselves. However, Chen et al. (2023)'s 's treatment specializes in the material point method (discretization) and elasticity (PDE); our general treatment is both discretization and PDE agnostic. Likewise, Pan et al. (2023) train a iiscretization-gnosiic latent space for PDE data; we further solve the PDES in the reduced space via rapid latent space traversal. Specifically, we evolve the latent space in PDE-constraint manner while Yin et al. 2023) take data-driven approach. Additional references to ROM are listed in Appendix AA]
Implicit neural representations use (fully- couley-co neural networks to represent arbitrary vector fields. While the Euclidean space spatial coordinates always form part of the input to the network, it is also common to have a latent space vector to complement the rest of the input. Different latent space vectors correspond to different states of the continuous vector field (see Figure Tp), e-g. different geometries (Park et al., 2019; Chen & Zhang 2019 Mescheder et al. 2019) or different radiance fields (Mildenhall et al. 2020). A key contribution of our work is nonlinearly traversing the latent space of neural representations under explicit PDE constraint.
Machine learning (ML) for PDEs. Physics-informed neural networks (PINNs) (Raissi et al. 2019) 2020b) demonstrate that PDEs can be accurately solved via neural representations. Notably, PINN enables prediction and discovery from incomplete models and incomplete data (Karniadakis al. 2021). However, the degrees of freedom involved in their approaches still Pd, and the underlying gradient-descent-based solver is often computationally more expensive than traditional solvers (see Table by Zehnder al. (202I)). By contrast, our goal is building computationally more efficient solution that solves for only r degrees of freedom (r < Pd). In fact, our approach can be viewed as an extension of PINN for model reduction. Setting the latent space vector (t) in our formulation (see Figure[Ip) to the time variable recovers the exact formulation of PINN. addition to PINN, PINN.(Sanchez- Sanchez-Gonzaler et al. (2020) show graph neural network (GNN) architectures are also capable of learning PDEs. Yet, like PINNs, GNNs also do not offer dimension reduction.
Sitzmann et al.
METHOD: OVERVIEW AND MANIFOLD CONSTRUCTION
Overview Our goal is to efficiently obtain the solution of Equation (I). We begin by constructing a low-dimensional embedding and the corresponding manifold- parameterization function (see below), after which we solve PDEs by time-integrating the dynamics of the manifold's latent space vector (see Section 4). As we demonstrate in examples from various scientific disciplines (see Section 5). this general method is applicable regardless of the discretization of the training data (e.g.. voxel grids, meshes, point clouds) or the discretization deemed most useful for evaluating the gradients (e.g.. physical forces) when solving PDEs on the constructed manifold.
Manifold Construction As depicted in Figure Ib, we seek manifold-parameteriine function g(x,q), q).
g(x, q(t: p)) se f(x, (4) Vo EN, VEET Vpe
(2)
that well approximates the continuous field f(x, t: μ1) throughout the spatiotemporal domain f x T, and for a workable range of problem parameters μ E D. For ease of exposition, we omit the dependencies of q and f on the problem parameters μ.. Here D is an arbitrary parameter space (e.g., material properties, external forces, user settings). For scenarios that do not feature trivial parameterizations, e.g., external force via crowd-sourcing (Barbië & James 2005), D can also be implicitly defined.
What distinguishes our approach from typical model reduction is that takes the position æ E S as an input. Thus, unlike prior approaches that infer only discrete coefficients, our approach infers field values at arbitrary domain positions æ.
We parameterize g with a neural network Io, whose weights 6 satisfy the minimization problem
T
Σ
l190, (a', qttt))) - f(0", tn)
(3)
1=1 n= ME Derain
3

---

Published as a conference paper at ICLR 2023
08g
4 (t.)
g0g
g0g
Encoder network ea.
Low-dimensional manifold gR,
Figure 2: Constructing the ranifolddparameetrrizaaioo function neural network trained via supervised learning. We pass each snapshot (time step) from the training dataset into an encoder to obtain a latent space vector q. We then concatenate q with the spatial coordinates and pass that into the
function with the goal of reconstructing f for each individual spatial sample. The same q shared among all spatial samples this time step.
where Dtrain D is the training set, and (t,,. is the latent space vector shared among all spatial samples. This objective aims to reproduce all the field values present in the training data, generated limit the discretization strategy of
via full-order PDE solutions. Notably, approach imposes the PDE solver. For instance, this framework is applicable to training data from both finite difference methods and finite element methods as well both voxel grids and meshes.
There are commonly two approaches to define the latent space vector q: the auto-decoder approach (Park al. 2019) that trains q along with 9e, and the encoder approach (Chen & Zhang 2019 Mescheder al.) 2019) that trains separate network to output q. While both approaches work for our application, we adopt the latter.
The encoder network Co, with weights 9, takes input vector constructed by concatenating all the discrete degrees of freedom from the training data and outputs a latent space vector (see Figure 2
eo, (f (t)) = q(t), where 7(1) (f(x), .),.. f(x' .t,.................. f(x 1))"
We emphasize that this discretization-depeeee smoothly varying latent space and for determining the initial latent space vector from the discretized training data. Otherwise, the implicit neural representation Is, remains discretivation-agnoiii architecture.
encoder (Xie al. 2021 is tool for training the
Adding the encoder, Equation now becomes
min 0g.0,
Σ Derann
11,,, (x' ea.
ee. (T())) f(x", f
Figure illustrates the training pipeline. Please refer to Appendix D for network and training details and Appendix F for hyperparameter selection.
4 METHOD: LATENT SPACE DYNAMICS
After the manifold is constructed, we compute the low-dimensional latent space dynamics (q,, H> In++) in three steps (see Figure 3 and appendix Q (1) network inference, (2) PDE time-stepping, and (3) network inversion. The neural network strictly serves as a kinematic spatial representation in the maps from and to the latent space (steps 1 and 3, respectively). The time integration (step 2) itself uses the original PDE, not a neural network approximation thereof. The ROM community has demonstrated that these three steps
yield strong long-time stability even stiff and chaotic (PaaeesssSdddddPedddddeeaa
dynamical systems (Carlberg et al., 20133 2017). Appendices O and dedemnntttatee smooth latent space trajectories and empirical stability analysis using this approach.
4

---

Published as a conference paper at ICLR 2023
An+1
1. Network inference
3. Network inversion
9
96
2. PDE time-stepping with boundary conditions
F Vy M
Figure 3: Latent space dynamics: temporally evolve from one latent space vector to another, governed by the PDE. The entire pipeline only involves degrees of freedom from a small spatial subset M, where M << P.
Commonly shared among all three steps are "integration samples", a finite set of spatial domain points M := {y' 12 S j < |M]} chosen at the user's discretion (see Section 4.4). These samples need not coincide with the previously mentioned full-order finite element discretization samples (x')
4.1 STEP NETWORK INFERENCE
We first aim to gather all the full-space spatiotemporal information (vy E M) necessary for PDE time integration. The function value f itself can be evaluated via inferencing of the neural network (y, 4n) 9o, (y. q..). The spatial and temporal gradients computed either by differentiating
the network, (f 4n) = y9e, and fly, (,) = 09%, mation. Higher-order gradients may be generalized in similar manner. Further details on gradient computation are listed in Appendix
ag (n. respectively, or by numerical approxi-
4.2 STEP 2: PDE TIME-STEPPING
We now evolve time from t, to tn+1. Unlike end-to-end learning-basse latent space dynamics methods (Lusch al. 2018), we evaluate time derivatives using the exact PDE (T. not a learned surrogate. At each integration point y, we evaluate the temporal derivative by solving the PDE (1) for f (n+1(y):
F(f.. V m. att f.t.....) = 0.
(5)
We evolve the configuration to time fn+1 ta + At using the chosen explicit time integration method If subject to given boundary conditions, e.g., Runge-Kutta methods (Dormand & Prince 1980):
fn+1 = II f J f..1....) Vy EMM
While this work focuses on explicit time integration, we can also extend the framework for implicit time integration (Carlberg et al. 2017).
4.3 STEP 3: NETWORK INVERSION
We project back onto the reduced manifold by finding the corresponding input J that best matches the evolved configuration f f+1 in a least-squares sense: (Quarteroni et al. 2014
min 2.+1ER ERR
l1go, (y. 9n+1) f(y. 111.
DEM
The objective is similar to the training loss found in Equation 3 but with two dimensions significantly reduced: the dimension of the unknown In+1, and the summation bound [MV. Consequently, instead of using stochastic gradient descent type of method, such as the auto-decoder scheme by Park et al. (2019), we achieve rapid inversion using the Gauss-Newton algorithm (Nocedal & Wright (006) with conditionally quadratic convergence. Further details are listed in Appendix
4.4 SPATIAL SAMPLE REDUCTION
necessary condition for the well well-posednesss the least squares formulation from Equation 7 is M|. Since the manifold- parmmeiidppaaaetterizaaiooo function construction guarantees that r Pd,
5

---

Published as a conference paper at ICLR 2023
ggoundd mith B 501 sunnclee IT's emer groady (oard) snic
Scand Energy Exoped Energy
Time Stape
Figure 4: Thermodynamics. Left: integration samples (M). Greedily selecting the samples (blue) allows us to use significantly fewer degrees of freedom than the full-order simulation (purple) while getting a much higher accuracy than naive uniform sampling (red). The depicted field is temperature governed by the heat equation after 100 time steps. Right: as dictated by the PDE, the reduced system conservative: total energy (stored thermal energy plus cumulative flux boundary) conserved.
choose ş S |Ml < P. To obtain the next-time step fly fn+1), Vy E M necessary for the least squares solves, we only require PDE updates (Section 4.2) and spatiotemporal data (Section 4.1 at these [M] samples. As such, the entire latent space dynamics framework (Figure 3) requires only [M] samples, compared to the full-order solver's s P samples. Hyper-reduction approaches like this have captured wide range of real-world scenarios, including massive elasticity deformations (Fulton et al. (019)) and large turbulent flows (Grimberg et al.) 202T). Unlike our dusereiizationniddeeeeeee approach, prior methods only support hyper-reduction samples that coincide with the full-order discretization.
naive selection of the integration samples can lead to inaccurate latent space dynamics, even if [M] 2 / refer to Figure 4 for the failure case of uniform sampling. As noted in the hyper-reduction literature, stochastic sampling can eliminate such errors (Carlberg) [201T). To better control hyper- reduction eтror, we draw inspirations from the cubature approach by An et al. (208) and propose a greedy algorithm that augments the sample set to meet target residual (see Appendix B). This strategy iteratively adds spatial samples to the sampling set until the target residual is met. In this case, the user decides the target accuracy, not the sample set size [M]. Alternatively, the user can terminate the algorithm early once it reaches chosen sample set size [Ml. Results of our sampling approach are shown in Figures and 20
EXPERIMENTS
We analyze the proposed framework on classic PDEs, with training data produced using a variety of discretizations (voxel grids, meshes, and point clouds). Unless otherwise noted, for each PDE, we delincate a testing set where Dteat C. D with Dirain n Preat = We construct the manifold (Section 3 with data from Derain and then validate the latent space dynamics (Section 4 Drest. We compare our approach with prior discretization-deee deeenddnt ROM methods, including POD (Berkooz et al. 1993 Holmes et al. 2012)) and neural -ewwwrk-wwoebee based autoencoder approaches (Fulton et al. 2019 Lee & Carlberg] 2020] (Shen et al.] (02T). under the same latent space dimensions. Additional implementation and reproducibility details are listed from Appendix to Appendix L Experiment statistics are summarized in Table [ The temporal evolutions of the PDEs are best illustrated via the supplementary video.
v(x) 84y = 0 Temperature u is governed by a one-dimensional heat equation. describes the apatiallyvarriin diffusion speed. Figure 4 displays our approach's 's ability to use very few integration samples and to capture conservation of energy.
Thermodynamics, के
((x)) 22 uu== 0. We model image blurring with the 2D diffusion equation (Perona & Malik 1990). Figure 5 shows that under the same latent space dimension (r 3), our method is more accurate than POD (Berkooz et al., 1993] Holmes et al.. (0122), both visually and quantitatively. With its architecture independent of pixel count, CROM uses an order of magnitude less memory than POD (see Figure 55))
Image Processing,
Transport dominated systems. Next, examine two transport-dominated slowly decaying Kol- mogorov n-widths problems, the Advection Equation and Burgers' Equation, where classic model reduction techniques often struggle (Peherstorfer 2022):
Advection Equation, du + (a V)u = 0. Here u is the advected quantity and a is the advection öt velocity. Figure 66 depicts CROM's favorable trajectory tracking relative both POD and the convolutional autoencoder (CAE) of Lee & Carlberg (2020), for equal latent space dimension (also Figure 21.

---

Published as a conference paper at ICLR 2023
(8)
ADCURACT: CURS VS. POO Srfes JPOO
Intical state
2248
(b) Accuracy comparison (*)
state
Memory comperison
Ground truth
POD
Figure 5: Image processing: comparison with POD. Ground truth solution uses / 65, 536 pixels. (a) Visually, our approach better captures the sharp initial state and the smoothed finial state than POD. (b) Quantitatively, CROM obtains a higher PSNR than POD (f the higher the better). (c) Our approach also uses order-of -fagnniuue less memory than POD (1 the lower the better). Both POD and our approach use the same latent space dimension (r 3).
I00
c.l1
(8) Advaciion equelion
(b) Burger's equation
Figure Breaking the Kolmogorov barrier: CROM outperforms POD and convolutional autoencoders
(Lee & Carlberg (2200) in tracking both the Advection and Burgers' trajectories. For both cases, use intrinsic colution-manifol dimension (Lee & Carlberg 2020). the lower bound of r, as latent space dimension To isolate the source of error, no hyper-reduction method is applied.
accurately captures the nonlinear dynamics than both POD and CAE (also see Figure 22p. Additionally, CROM
Burgers' Equation, Ow 5v 0. Figure 60 shows that CROM
uses 12x less memory than CAE (see Figure 27 These accuracy and memory advantages
consistent with the implicit neural representation literature (Chen & Zhang 2019).
Whereas CAE applicable only voxel grid data, CROM is applicable to data from any discretiza- tion, as highlighted by the tetrahedral mesh and point cloud data examples (see Figure
Vju = - Vp + v + fext: V u = 0. CROM captures the challenging Karman vortex street (see Figure governed by the incompressible Navier-Stokes equations. It even generalizes well beyond the training temporal range.
Su + (u Fluid Dynamics, øt
Steady state
Bifurcation
Vortex shedding regime
Ground Truth
Ours
Training temporal range
Generalize beyond training temporal range
Figure 7: The Karman vortex street. Our latent space dynamics solver accurately captures the flow both before and after the bifurcation: the steady state and the periodic vortex shedding regime. The visualized quantity is the the curl the velocity field. Our approach even generalize well beyond the training temporal range, thanks to the kinematies-apprxxiattor
nature of CROM.
7

---

Published as a conference paper at ICLR 2023
stp
OVER PRIOR ROM 48X
Wall-clock Speedup
REMESHI
35X
CPU thread
CPU
200
lociigg
high-res
Adaattive
visualization
Inaccurate
accurate slow
Reccurate
Prior approach
Our appreach
Figure 8: Adaptive discretization in solid mechanics. A falling deformable body impacts two static objects. (a and b) Prior approaches (Barbič & James 2005 Fulton et al. 2019; Shen et al. 2021 only support reduced-order dynamics with fixed discretizations, forcing a tough choice between accuracy and computational cost; (c) approach allows dynamic mesh adaptivity to balance accuracy and cost as problem difficulty varies over time. (d) Our resolution-indepedeent encoding of the deformation map allows to visualize the field using any method (e.g.. a high resolution mesh) without regard to the computational mesh. CROM's speedup over prior discretization discretization-deeeneeet ROM methods observed for 1-32 threads, and emphasized the case limited computational resources. To ensure representativeness of prior ROM architectures, we report timing strictly for the PDE-time-stepping stage (Step 2), which is shared among all prior ROM architectures.
CROM'S SPEEDUP
OVER FULL-ORDER MODEL
Wall-clock Speedup
a1-109X fanter than truth 45%
26X then
B8X
5.7..
approach
41X
NVIDIA V100
2.30GHZ
mesh and Integration camples (d) Intial condtion
(c) Natre downerselution
(1 Ours different elastc modul
(g) Point dloud
Figure 9: Solid mechanics (a) The ground truth is generated via the full-order PDE solver. (b) Our approach is 41-109x faster than the ground truth while capturing detailed shearing and volume- preserving hehaviors (1.46% error). (c) Prior approach (Barbič 2012) of the same latent space dimension consumes 26x more memory while suffering from volume-gain artifacts (5. 76% error). (d) These simulations adopt a tetrahedral discretization (d, pink mesh, /' = 66, 608). Instead of using the expensive high-resolution mesh, our approach computes dynamics using very few integration samples (d, colorful spheres, |M| = 40). (e) Naive down-resolution of the ground truth simulation yields similar runtime but leads to significantly worse quality. (f) After training, our model can capture wide range of material properties. (g) The same manifold-parameee parameterization function architecture (different network weights) can also be used for model-reducing point-cloud based simulation (reproduced from the work by (Chen et al., 2023)). The speedup plot on the right demonstrates the effectiveness of our approach on diverse computing platforms. Disclaimer: the authors do not support animal cruelty.
Solid Mechanics, oo0 V P(VØ) + B We solve the second-order elastodynamics equation for the deformation map & of soft bodies, where Po is the initial density, P is the first Piola-Kirchhoff
stress, and B the body force density (see Figures 8 to 10 and 25
Adaptive resolution. Prior ROM methods (see Figure Th) fix the discretization at onset, precluding the dynamic adaptivity approaches that benefit problems with evolving complexity, When modeling falling block (see Figure 8) prior ROM methods must begin with a high-resolution mesh, despite that the high-resolution required only during contact (Figure & and b). Our wwwpetiiaaaiionnndeeeeeeeee representation allows the time integrator to freely adapt the spatial discretization throughout the dynamic evolution. For instance, we can employ a coarse mesh (P = vertices) to economize computation during the rigid falling phase (Figure (b)) and high-resolution mesh (P 065) to capture details during contact (Figure 8b), yielding 25x-48 x8x wall-clock speedups over fixed discretization ROM methods. To ensure generalizability of the comparison, the timing is done strictly for the PDE time-stepping stage (Step 2). Prior ROM architectures (i.e, both POD and autoencoder approaches) share the same PDE time-stepping routine.
8

---

Published as a conference paper at ICLR 2023
Mannry (umnceed Paramelend
10047
Recirsmation Shot Rasenni
1004
10044
0.0001
10044
10040
0.00001
Space
Lonert Space D monsion
--OROM (Duss)
--POD
--Fuuton
2019] -(Shin
(Duns)
sondl
|Sthen 20e1)
Figure 10: Our approach uses significantly less memory and simultancously yields higher accuracy than POD and the autoencoder-based ROM methods Fulton et al. (2019) and SSen et al. (2021) for the solid mechanics experiment depicted in Figure For example, with 3, our approach uses 39x x and x32x less memory, while simultaneously offering x9x and 49x more accuracy, compared POD and the autoencoder approach Fulton al. (2019), respectively.
Memory footprint. Most prior ROM approaches are discretizationndeeeedeet. As discretization resolution (P) increases, so does memory consumption for g. By contrast, our discretization- independent architecture does not "see" the parameter P and in practice memory footprint does not scale with P. Rather, as with all neural fields (Xie et al. 202T). memory grows with the intrinsic complexity of the training data.
Figure compares memory footprint and reconstruction error for CROM, POD (Barbič 2012). and based autoencoder approaches (Fulton et al. 2019 Shen et al. (021). as a function of latent space dimension; Figure
Figure@Jeeiccss depicts the corresponding large deformation simulations. Figure compares CROM's 's memory and accuracy another large deformation example. For these examples, CROM offers simultancously lower memory consumption and reconstruction error. Furthermore, due to hyper-reduction (from P = 66, 608 to [M| 40), CROM also obtains significant wall-clock speedups on CPUs, consumer GPUs, and data center GPUs over the unreduced model (see Figure 9
Different discretizations, same architecture. Because CROM is not constructed for fixed dis- cretization, we can adopt an identical manifold-paramette parameterization function architecture to model- reduce simulations of both tetrahedral meshes and point clouds (see Figure (g). Moreover, even though these examples (Figures g and leverage completely different meshes, they use the same network architecture. This opens the door for future work transfer learning (Weiss et al. 2016) and weight sharing among various discretizations and PDEs.
Additional comparisons between CROM, previous ROM, and the full-order models are discussed in Appendices M andd] N
6 DISCUSSION AND CONCLUSION
CROM a model reduction framework for PDEs featuring a discretization dissretization-nndeeeeeen architecture. CROM disentangles the low-dimensional embedding from the discretization by reformulating the reduced manifold as a map accepting not only q but also æ as an input. CROM outperforms discretization-deeendee ROM approaches, such as POD and autoencoders, in terms of accuracy and memory consumption for equal latent space dimension, and in the ability to dynamically adapt the discretization during time integration.
While offering key advantages over prior ROM methods, CROM also inherits limitation of ROM. can only treat PDE solutions in the space spanned by the manifold-parameetrizaaoo function (determined by the training data) and does not generalize to arbitrary unseen scenarios. Such limitation also commonly found among other implicit neural representation works (Xie et al. 202T). Future research may consider improving generalizability and data efficiency via meta-learning and integration of stronger priors (Sitzmann al. 2020a).
Compared to end-to-end ML solutions to PDEs "Sanchez-Gonzaleee et al. 2020). CROM employs the neural network strictly as a spatial representation (see Sections 4.] and 4.3) for the kinematics and solves the PDE using classical PDE-integration numerical methods (see Section 4.2). As such, we believe CROM will open doors for more forthcoming hybrid ML-PDE solutions. As shown in work, these solutions integration with existing PDE solvers (see Section 4.2), and obtain practical computational savings that can be directly employed in production (see Figure
retain the PDE's physical invariants (see Figure 4 allow for easy

---

Published as a conference paper at ICLR 2023
ACKNOWLEDGMENTS
This work was supported in part by the National Science Foundation (Grant CBET-17-06689,, Meta, Natural Sciences and Engineering Research Council of Canada (Discovery Grant), and SideFX. We thank Honglin Chen and Rundi Wu for their insightful discussions. We thank Rundi Wu for sharing his implementation of Stable Fluids. We thank Raymond Yun Fei for rendering advice. We thank Keenan Crane for the "Origins of the Pig" mesh.
10

---

Published as a conference paper at ICLR 2023
REFERENCES
Steven S An, Theodore Kim, and Doug L. James. Optimizing cubature for efficient integration of subspace deformations. ACM transactions on graphics (TOG), 27(5):1-10, 2008.
C Bach, Linze Song, Tatjana Erhart, and Fabian Duddeck. Stability conditions for the explicit integration of projection based nonlinear reduced-order and hyper reduced structural mechanics finite element models. arXiv preprint arXiv: 1806. / 1404, 2018.
Jernej Barbič. FEM simulation of 3D deformable solids: A practitioner's guide to theory, discretiza- tion and model reduction. Part 2: Model reduction. In SIGGRAPH 2012 Course Notes. Citeseer, 2012.
Jernej Barbič and Doug L James. Real-time subspace integration for St. Venant-Kirchhoff deformable models. ACM transactions on graphics (TOG), 24(3)9882-99,, 2005.
Ulrike Baur, Christopher Beattic, Peter Benner, and Serkan Gugercin. Interpolatory projection methods for parameterized model reduction. SIAM Journal on Scientific Computing, 33(5): 2489-2518, 2011.
Peter Benner, Serkan Gugercin, and Karen Willcox. A survey of projection-based model reduction methods for parametric dynamical systems. SIAM review, 57(4):483-531, 2015.
Gal Berkooz, Philip Holmes, and John Lumley. The proper orthogonal decomposition in the analysis of turbulent flows. Annual review of fluid mechanics, 25(1)5539-575 1993.
Don O Brush, Bo O Almroth, and JW Hutchinson. Buckling of bars, plates, and shells. 1975.
John Rozier Cannon. The mie-dimensional heat equation. Number 23, Cambridge University Press, 1984.
Kevin Carlberg, Charbel Farhat, Julien Cortial, and David Amsallem. The GNAT method for nonlinear model reduction: effective implementation and application to computational fluid dynamics and turbulent flows. Journal of Computational Physics, 242:623-644 647,2013. 2013.
Kevin Carlberg, Matthew Barone, and Harbir Antil. Galerkin least-squares Petrov-Galerkin projection in nonlinear model reduction. Journal of Computational Physics, 330:693-734, 2017.
Kevin Thomas Carlberg. Model reduction of nonlinear mechanical systems via optimal projection and tensor upproximation. PhD thesis, Stanford University, 2011.
Peter Yichen Chen, Maurizio M Chiaramonte, Eitan Grinspun, and Kevin Carlberg. Model reduction for the material point method via an implicit neural representation of the deformation map. Journal of Computational Physics, 478:1 1908, 2023.
Zhiqin Chen and Hao Zhang. Learning implicit fields for generative shape modeling. Proceedings of the IEEF/CVF Conference on Computer Vision and Pattern Recognition, pp. 9939-5948, 2019.
Alexandre Joel Chorin. Numerical solution of the Navier-Stokes equations. Mathematics of computa- tion, 22(104):745-762, 1968.
Pi-Yueh Chuang and Lorena A Barba. Experience report of physics-informed neural networks in fluid simulations: pitfalls and frustration. arXiv preprint arXiv:2205. 2022.
Djork-Arné Clevert, Thomas Unterthiner, and Sepp Hochreiter. Fast and accurate deep network learning by exponential linear units (elus). arXiv preprint 1.11::25 2015.
Dylan Matthew Copeland, Siu Wun Cheung, Kevin Huynh, and Youngsoo Choi. Reduced order models for Lagrangian hydrodynamics. Computer Methods in Applied Mechanics and Engineering, 388: 1 14259, 2022.
Roy R Craig Jr and Mervyn CC Bampton. Coupling of substructures for dynamic analyses. AIAA journal, 6(7):1313-319, 1968.
11

---

Published as a conference paper at ICLR 2023
John R Dormand and Peter J Prince. A family of embedded Runge- Kutta formulae. Journal of computational and applied mathematics, 6(1):19-28 19:26,,9,6,11 1980.
Yilun Du, Katie Collins, Josh Tenenbaum, and Vincent Sitzmann. Learning signal-agnostic manifolds of neural fields. Advances in Neural Information Processing Systems, 34:8320 8331, 2021.
Virginie Ehrlacher, Damiano Lombardi, Olga Mula, and François-Xavice Xavier Vialard. Nonlinear model reduction on metric spaces. application to one-dimensional conservative pdes in wasserstein spaces. ESAIM: Mathematical Modelling and Numerical Analysis, 54(6):2159-25 2197, 2020.
N Benjamin Erichson, Michael Muchlebach, and Michael W Mahoney. Physics- Physics-informen informed autoen- coders for Lyapunov-stable fluid flow prediction. arXiv preprint arwiv: 10866, 2019.
William Falcon et al. PyTorch Lightning. GitHuu. Note: https///ghhuu.coo/////fftttthhhs lighriing, 3, 2019.
Lawson Fulton, Vismay Modi, David Duvenaud, David IW Levin, and Alec Jacobson. Latent- space dynamics for reduced deformable simulation. Computer graphics forum, volume 38, pp. 379-391. Wiley Online Library, 2019.
Sebastian Grimberg, Charbel Farhat, Radek Tezaur, and Charbel Bou-Mosleh. Mesh sampling and weighting for the hyperreduction of nonlinear Petrov-Galerkin reduced-order models with local reduced-order bases. International Journal for Numerical Methods in Engineering. 122(7): 1846-1874, 2021.
Chenjie Gu. Model order reduction of nonlinear dynamical systems. Berkeley, 2011.
University of California,
Serkan Gugercin, Athanasios C Antoulas, and Christopher Beattie. H2 model reduction for large- scale linear dynamical systems. SIAM journal on matrix analysis and applications, 30(2):609-66 638, 2008.
David Hartman and Lalit K. Mestha. A deep learning framework for model reduction of dynamical systems. In 2017 IEEE Conference on Control Technology and Applications (CCTA), pp. 1917- 1922, 2017. doi: 10.1 11119//CCAA22177866775333
Amir Hertz, Or Perel, Raja Giryes, Olga Sorkine-Hornung, and Daniel Cohen-Or. Sape: Spatially- adaptive progressive encoding for neural optimization. Advances in Neural Information Processing Systems, 34, 2021.
Philip Holmes, John I Lumley, Gahl Berkooz, and Clarence W Rowley. Tarbulence, coherent structures, dynamical systems and symmetry. Cambridge university press, 2012.
Yixin Hu, Qingnan Zhou, Xifeng Gao, Alec Jacobson, Denis Zorin, and Daniele Panozzo. Tetrahedral meshing in the wild. ACM Trans. Graph., 37(4))60-1, 2018.
Patrick Huerre and Peter A Monkewitz. Local and global instabilities in spatially developing flows. Annual review of fluid mechanics, 22(1):473-537, 1990.
Thomas JR Hughes. The finite element method: linear static and dynamic finite element analysis. Courier Corporation, 2012.
George Em Karniadakis, Ioannis G Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang. Physics-informet machine learning. Nature Reviews Physics, ((6):422-44)) 2021.
Kenji Kashima. Nonlinear model reduction by deep autoencoder of noise response data. 2016 IEFE 55th Conference Decision and Control (CDC), pp. 5750 5755, 2016. doi: 10.1 2016.7799153.
Byungsoo Kim, Vinicius Azevedo, Nils Thuerey, Theodore Kim, Markus Gross, and Barbara Solenthaler. Deep fluids: A generative network for parameterized fluid simulations. In Computer Graphics Forum, volume 38, pp. 59-70. Wiley Online Library, 2019.
Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv: 148806980,... 2014.
12

---

Published as a conference paper at ICLR 2023
Karl Kunisch and Stefan Volkwein. Galerkin proper orthogonal decomposition methods for general equation in fluid dynamics. SIAM Journal on Numerical analysis, 40(2):492-515, 2002.
Kookjin Lee and Kevin Carlberg. Deep conservation: A latent dynamics model for exact satisfaction of physical conservation laws. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 35, pp. 277-285, 2021.
Kookjin Lee and Kevin T Carlberg. Model reduction of dynamical systems on nonlinear manifolds using deep convolutional autoencoders. Journal of Computational Physics, 404: 108933, 2020.
David B Lindell, Julien NP Martel, and Gordon Wetzstein. Autoint: Automatic integration for fast neural volume rendering. In Proceedings of the IEEE/CVE Conference on Computer Vision and Pattern Recognition, pp. 14556-14563 2021.
Hsuch-Ti Ti Derek Liu, Francis Williams, Alec Jacobson, Sanja Fidler, and Or Litany, Learning smooth neural functions via lipschitz regularization. In ACM SIGGRAPH 2022 Conference Proceedings, pp. 1-13, 2022.
Lingjie Liu, Jiatao Gu, Kyaw Zaw Lin, Tat-Seng Chua, and Christian Theobalt. Neural sparse voxel fields. Advances in Neural Information Processing Systems,
2020.
Bethany Lusch, Nathan Kutz, and Steven I. Brunton. Deep learning for universal linear embeddings of nonlinear dynamics. Nature communications, (1):490,, 2018.
Julien NP Martel, David B Lindell, Connor Z. Lin, Eric Chan, Marco Monteiro, and Gordon Wetzstein. Acorn: Adaptive coordinate networks for neural scene representation. arXiv preprint arXiv:2105.02788. 2021.
Romit Maulik, Arvind Mohan, Bethany Lusch, Sandeep Madireddy, Prasanna Balaprakash, and Daniel Livescu. Time-series learning of latent-space dynamics for reduced-order model closure. Physica D: Nonlinear Phenomena, 403618 2020.
Romit Maulik, Bethany Lusch, and Prasanna Balaprakash. Reduced-orde modeling of advection- dominated systems with recurrent neural networks and convolutional autoencoders. Physics of Fluids, 33(3):037106 2021.
Lars Mescheder, Michael Oechsle, Michael Niemeyer, Sebastian Nowozin, and Andreas Geiger. Occupancy networks: Learning 3d reconstruction function space. In Proceedings of the IEEE/CVF Conference Computer Vision and Pattern Recognition, pp. 4460-4470, 2019.
Ben Mildenhall, Pratul P Srinivasan, Matthew Tancik, Jonathan T Barron, Ravi Ramamoorthi, and Ren Ng. Nerf: Representing scenes neural radiance fields for view synthesis. In European conference on computer vision, pp. 405-421. Springer, 2020.
B.C. Moore. Principal component analysis in linear systems: Controllability, observability. and model reduction. IEÉE Transactions on Automatic Control, 26(1):17-32, 1981. doi: 10. 1211099TAC.1999111105555 httts:/t::/// //www. scopus com httmm//wwawddddemmmm wtw.://...dree.......... 2io=22-22.00000155333-
11w.:1......999......... .10.c..9819 1981. 11 w11055684aatterriiiii 40 400m5622266377022221222222288833
cited By 3526.
Thomas Müller, Alex Evans, Christoph Schied, and Alexander Keller. Instant neural graphics
primitives with a multiresolution hash encoding. ACM Transactions on Graphics (ToG), 41(4): 1-15, 2022.
Mohammad Sina Nabizadeh, Stephanie Wang, Ravi Ramamoorthi, and Albert Chern. Covector fluids. ACM Transactions on Graphics (TOG), 4164):1-16, 2022.
Jorge Nocedal and Stephen Wright. Numerical optimization. Springer Science & Business Media, 2006.
Mario Ohlberger and Stephan Rave. Nonlinear reduced basis approximation of parameterized evolution equations via the method of freezing. Comptes Rendus Mathematique, 351(23-24): 901-906, 2013.
13

---

Published as a conference paper at ICLR 2023
Shaowu Pan and Karthik Duraisamy. Data-driven discovery of closure models. SIAM Journal Applied Dynamical Systems, 17(4):2381-2413, 2018.
Shaowu Pan, Steven L Brunton, and J Nathan Kutz. Neural implicit flow: mesh-agnostic dimen- sionality reduction paradigm of spatio-tempores data. Journal of Machine Learning Research, 24 (41):1-60, 2023.
Jeong Joon Park, Peter Florence, Julian Straub, Richard Newcombe, and Steven Lovegrove. Deepsdf: Learning continuous signed distance functions for shape representation. Proceedings of the IEEE/CVF Conference Computer Vision and Pattern Recognition, pp. 165-174, 2019.
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al. Pytorch: An imperative style, high-performance deep learning library. Advances in neural information processing systems, 32, 2019.
Benjamin Peherstorfer. Model reduction for transport-dominatee problems via online adaptive bases and adaptive sampling. SIAM Journal on Scientific Computing, 42(5): A2836, 2020.
Benjamin Peherstorfer. Breaking the kolmogorov barrier with nonlinear model reduction. Notices of the American Mathematical Society, 69(5):725-733 2022.
Benjamin Peherstorfer and Karen Willcox. Online adaptive model reduction for nonlinear systems via low-rank updates. SIAM Journal on Scientific Computing, 37(4) 2015.
J Peraire, M Vahdati, K Morgan, and O.C Zienkiewicz. Adaptive remeshing for compressible
flow computations. Journal Cf Computational Physics, 72(2):449-466, 1987. ISSN 0021-9991. doi: httpl///ceeeeeerreeee//////////0 com/
URL BtLhttts:/// www. sc wwt.://w.....n.cc.cmmcom. enced rect
Pietro Perona and Jitendra Malik. Scale-space and edge detection using anisotropic diffusion. IEEE Transactions on pattern analysis and machine intelligence, 12(7):629-66 1990.
C. Prud'hone homme, D.V. Rovas, K. Veroy, L. Machiels, Y. Maday, A.T. Patera, and G. Turinici. Reliable real real-time solution of parametrized partial differential equations: Reduced-basis output bound methods. Journal of Fluids Engineering, Transactions the ASME, 124(1):70-80, 2002. doi: 10.1115/1. 1448332. URL https: wwwwwwwwwooooo scwpus.cuuuuuoo com wccoppas.commmwwaarrreer........
1115$2f1 1a0.55221111122222aaarrrrr0-0000 2311003a000003222333322222222220
eid= 00002222200000332200 00000000030000000000000000000000
cited By 338.
Elizabeth Qian, Boris Kramer, Benjamin Peherstorfer, and Karen Willcox. Lift & learn: Physics- informed machine learning for large-scale nonlinear dynamical systems. Physica D: Nonlinear Phenomena, 406: 132401 2020.
Alfio Quarteroni, Gianluigi Rozza, et al. Reduced order methods for modeling and computational reduction, volume 9. Springer, 2014.
Maziar Raissi, Paris Perdikaris, and George Karniadakis. Physics-informed neural networks: deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378:686-707, 2019.
Francesco Regazzoni, Luca Dede, and Alfio Quarteroni. Machine learning for fast and reliable solution of time-dependent differential equations. Journal of Computational physics, 197:10 2019.
Cristian Romero, Dan Casas, Jesús Pérez, and Miguel Otaduy. Learning contact corrections for handle-based subspace dynamics. ACM Transactions on Graphics (TOG), 40(4):1 2021.
G. Rozza, D.B.P. Huynh, and A.T. Patera. Reduced basis approximation and posteriori error estimation for affinely parametrized elliptic coercive partial differential equations. Arch. Comput. Methods Eng., 15(3):1-47, 2007. cited By 27.
14

---

Published as a conference paper at ICLR 2023
Alvaro Manchez-Gonzalee Jonathan Godwin, Tobias Pfaff, Rex Ying, Jure Leskovec, and Peter Battaglia. Learning to simulate complex physics with graph networks. International Conference on Machine Learning. pp. 8459 8468. PMLR, 2020.
Siyuan Shen, Yin Yang, Tianjia Shao, He Wang. Chenfanfu Jiang, Lei Lan, and Kun Zhou. High-order differentiable autoencoder for nonlinear model reduction. ACM Trans. Graph., 40(4), jul 2021. ISSN 0730-0301. doi: 10.114//3550626345557555 URL https:://oi..oooeeeeeeo http://///oooooo1111///10 041006666.459755
Eftychios Sifakis and Jemnj Barbic. Fem simulation of 3d deformable solids: A practitioner's guide to theory, discretization and model reduction. In ACM SIGGRAPH 2012 Courses, SIGGRAPH '2, New York, NY, USA, 2012. Association for Computing Machinery. ISBN 9781450316781. doi: 10.1145/233383.2243300 URL (RLLttts:///oo.. htt.:/......//111111222224433332 2343500
Vincent Sitzmann, Eric Chan, Richard Tucker, Noah Snavely, and Gordon Wetzstein. Metasdf: Meta-learning signed distance functions. Advances in Neural Information Processing Systems, 33: 10136-10147, 2020a.
Vincent Sitzmann, Julien Martel, Alexander Bergman, David Lindell, and Gordon Wetzstein. Im- plicit neural representations with periodic activation functions. Advances in Neural Information Processing Systems, 33, 2020b.
Olga Sorkine and Marc Alexa. arrigid-as-possibee surface modeling. In Symposium Geometry processing, volume 4, pp. 109-116, 2007.
Jos Stam. Stable fluids. In Proceedings of the 26th annual conference on Computer graphics and interactive techniques, pp. 121-128, 1999.
Tommaso Taddei, Simona Perotto, and ALFIO Quarteroni. Reduced basis techniques for nonlinear conservation laws. ESAIM: Mathematical Modelling and Numerical Analysis, 49(3):787-814, 2015.
Towaki Takikawa, Joey Litalien, Kangxue Yin, Karsten Kreis, Charles Loop. Derek Nowrouzezzarrai Alec Jacobson, Morgan McGuire, and Sanja Fidler. Neural geometric level of detail: Real-time rendering with implicit 3D shapes. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 11358-11367, 2021.
Han Vanholder. Efficient inference with tensorrt. In GPU Technology Conference, volume 1, pp. 2016.
Karl Weiss, Taghi M Khoshgoftaar, and DingDing Wang. A survey of transfer learning. Journal of Big data, 3(1):1-40, 2016.
Karen Willcox and Jaime Peraire. Balanced model reduction via the proper orthogonal decomposition. AIAA journal, 40(1 (0(11):2323-2330 2002.
Yiheng Xie, Towaki Takikawa, Shunsuke Saito, Or Litany, Shiqin Yan, Numair Khan, Federico Tombari, James Tompkin, Vincent Sitzmann, and Srinath Sridhar. Neural fields in visual computing and beyond. arXiv preprint arXiv:2111.11426, 2021.
Guandao Yang, Serge Belongie, Bharath Hariharan, and Vladlen Koltun. Geometry processing with neural fields. Advances in Neural Information Processing Systems, 34, 2021.
Yuan Yin, Matthieu Kirchmeyer, Yves Franceschi, Alain Rakotomamonjy, and patrick gallinari. Continuous PDE dynamics forecasting with implicit neural representations. In The Eleventh International Conference on Learning Representations, 2023, URL /ttl:tts///peereeieee net nett/Eorum220 id- 373
Jonas Zehnder, Yue Li, Stelian Coros, and Bernhard Thomaszewski. Ntopo: Mesh-free topology optimization using implicit neural representations. Advances Neural Information Processing Systems,
2021.
15

---

Published as a conference paper at ICLR 2023
Appendix
We highly encourage viewing the supplementary video where our results are best illustrated.
A ADDITIONAL LITERATURE REVIEW
A.1 REDUCED-ORDER MODELING
Accelerating PDE solutions via reduced-order modeling (ROM) has a rich history. Two common building blocks of ROM are (1) identifying a low-dimensional latent space of the original complex system and (2) solving PDEs via latent space dynamics, i.e., evolving the latent space vector over time.
Early works on latent space construction started with linear methods where the the latent space and the full space has linear relationship. Notable works include proper orthogonal decomposition
(Berkooz et al. 1993] Holmes et al. 2012 Barbič & James 2005 Barbič 2012). the reduced-basse
technique (Prud homme et al. 2002, Rozza al. 2007]. balanced truncation (Moore 1981), rational
interpolation (Baur et al. 2011 Gugercin et al. 2008), and
model reduction CCrigg Jr
& Bampton 1968). Recently. nonlinear methods (Hartman & Mesthal 2017 Kashimal 2016 Gu
2011 Brichson et al. 2019: Maulik et al.
2020 Regazzoni et 2019 Fulton et al. 2019 Maulik
al. |021] Romero al. 2021 [Shen et al.
2022 often constructed via neural networks, have gained
significant attentions due to their ability to more accurately construct the wwwpffffddparammterizzttiie
function. In particular. nonlinear manifolds significantly outperform their linear counterparts on
slowly decaying Kolmogorov n-width problems (c.g., advection) (Peherstorfer, 2022 Ohlberger &
Rave 2013: Peherstorfer & Willcox 2015
Taddei al. 2015 Peherstorfer 20200] Ehrlacher et al.
2020
Lee & Carlberg 2021).
dicretized vector fields of interest. [Cen et al. (2023); ; Pan et al. (202)) and our work are the first attempts at constructing the latent space for continuous vectors themselves, thereby allowing the training data to come from any discretizations.
Most of these prior works exclusively focus on building latent space for the
After the latent space is identified, solving PDEs with parameters unseen during training requires evolving the latent space vector over time. Two families of approaches have been proposed for latent space dynamics. The first family of approaches (Kim et al. 2019) computes latent space dynamics via a data-driven operator that is learned from the training data. A key benefit from learning the operator from data is that the user requires no knowledge of the system's governing equations and does not need access the training solver's source code. However, the downside of this approach is that additional assumptions (e.g., the availability of a Koopman invariant subspace as in Lusch finite-dimensional quadratic representation as is in the lift and learn approach by Qian et al. (2020)) are made. Without these assumptions, the learned operator prone violation of physical laws. The second family of approaches (Lee & Carlberg 2020 Fulton et al. 2019). to which our method belongs, computes latent space dynamics according to the PDE used for training data generation. The advantage of this approach is that it does not require learning an additional time-evolution operator. It works with arbitrary dynamical systems with known PDE governing equations and conserves physical laws by construction since explicit PDE model employed. However, the disadvantage is that the PDE model must be known.
al. (2018) and the existence
B ROBUST SAMPLING
To select the integration samples robustly, we aim to balance the computation speed and accuracy. CROM employs [M] spatial samples, and the computation speed roughly scales linearly with the number of samples. To maximize computation speed, we aspire to select the fewest number of samples possible.
16

---

Published as a conference paper at ICLR 2023
Given target accuracy, our greedy algorithm selects the fewest number of integration samples from the discrete samples {a [x')== of the full-order PDE solution order to achieve the target accuracy.
Algorithm 1: Robust sampling
M ={**) // randomly initialize integration sample set, l E (1,...,
while true do
res_vec= = Calculate, Cclcuuaae_Resiitudl....
-Metric(res__eer _ec) <
then
return .M
else
indices Max aaa__Indicee(rssserrr Q) find Q indices with the largest individual
residuals
_ist_littt 0 Rº for in indices do res_ res_list[i Metric(Calluddteeeesssddd ResiiuulleeeiiiaaMMMMMM U (x')))
end
j = Min Inddinnccccess_saaeee _ist, 1) // find the index with the smallest global residual
12
18
M (
14
end
15
end
Function Colculate Colextole_Resiiull(MM.. R
16
17
res_VeC
for Dirain do Obtain the final latent space vector q(tr: (4) with the integration samples (M) for i 1,..... do
19
res res_vec[i]++ 19s, (x', qqt:::::t)) f(x', tr
end end
23
24
return res_vec
25
Function Metric(res_ :
26
return mean(res_vent + mar(res__ vec)
In every iteration, the greedy algorithm adds one spatial sample to the sample M and lowers the error of latent space dynamics. Specifically, the algorithm loops over Q spatial samples with the largest individual residuals. We test adding each sample the integration sample set and compute a new global residual. From these Q samples, we select the one leading to the smallest new global residual and add it to the actual sample set. The algorithm repeats until target accuracy met. In practice, we find Q = 10 gives sufficient results. Algorithm|] [ describes the details. Compared to the naive uniform sampling approach, our robust sampling scheme yields significantly more accurate results (see Figure 4.. We also highlight that Algorithm happens before any simulation occurs (latent space dynamics). strictly a pre-computation.
C NETWORK INVERSION
In order to generate latent space dynamics, need invert the neural network employed for the manifold-parameteriaen function. We do so by solving the minimization problem:
min a ER
Ilo,, (n+1) f (n))2..
DEM
Since this is a nonlinear least- squatesquares problem, we solve with the classic iterative Gauss-Newtar algorithm (Nocedal Wright 2006). Unlike Newton's method, Gauss-Newton does require the
Hessian of the objective function. Consequently, in addition to the neural network itself, we only need
0ge, to evaluate the gradient of the neural network with respect to the latent space vector Bq (y. 9-+1)- We also equip our Gauss-Newton implementation with a standard backtracking line-seardh scheme to improve the convergence rate.
17

---

Published as a conference paper at ICLR 2023
C.1 INITIAL GUESS FOR THE NONLINEAR SOLVER
An initial guess close the minimum is necessary for the convergence of Gauss-Newton. The previous time step latent space vector q, serves as an ideal initial guess since does not vary much over time step. However, in the extreme cases where f deviates significantly over time step, find it helpful to compute the initial guess through an encoder network. For convenience, we employ the encoder network from training, which requires computing the dynamics for P samples. Future work should consider building this encoder network strictly for the degrees of freedom in M. We can further facilitate nonlinear solver convergence by encouraging smoothness of the latent space, such as employing Jacobian penalty term (Chen et al., 2023) / local isometry term (Du et al.) (202T) and enforcing Lipschitz constraints (Liu et al. 2022). These smoothness regularizations would be particularly important for challenging PDE cases involving bifurcations and buckling (Brush et al. 1975).
C.2 LINEARIZATION
We can further obtain computation-san save by bypassing the iterative nonlinear solver and linearizing the least-squares problem.
We define f as,
Af(y,bat1) = fly (n+1) - 5(y, (n) = fy, (n+1) Iu, (y, an)
(9)
(10)
Taylor expansion 90, at qn we have go, (3+ 9n+1) Oges (y, 9n An+11 where an+1 In + Aqn+1- Therefore, the original objective function can be approximated as,
Employing
go, (y. q,,)
ge, (y, an+1) ((y,knt) = Ie, (y, (n+1) ) 90, (y: 9n)) - A5 bn+1)
(11)
0ga, (y. q., Am)4 tn41)))))
(12)
Ōq
Therefore, we can obtain the latent space dynamics by solving the linear least-squares problem,
Σ 11 Oge, (y. qm)))amT1 Af(g)t41111
min MER"
(13)
0q
DEM
This is the normal equation and can solved in the closed form:
Aq (J J)
(14)
agaz oq
where J is the (IM) d) by r Jacobian matrix that contains by 1 residual vector containing Af fn41), Vy e M.
(y. a..,, E M; the (IM)-
Since our work features both small |M| and small r, the closed-form network inversion costs just few small dense matrix multiplications and thereby introduces minimal overhead. In fact, the solver itself is extremely efficient, and the majority the computation cost forming the Jacobian matrix J itself, whose optimization will be discussed in Appendix E.I.I Additional performance details listed Appendix N
Although network inversion step follows the standard optimal-projection-besss ROM literature (Carlberg et al. (013)), future work may consider directly learning an inverse during training time. The encoder introduced in the manifold construction process may be step towards learning inverse. However, the encoder its current form does not have optimal accuracy and is bounded to a particular discretization. Future work should consider lifting these limitations.
18

---

Published as a conference paper at ICLR 2023
D NETWORK DETAILS
Consistent with the implicit neural representation literature, we parameterize the manifold- parameterization function manifold with MLP network 98, The input dimension of the network is m while the output dimension of the network is d, where is the dimension of the input spatial vector while d is the output dimension of the vector field of interest. Our MLP network contains hidden layers, each of which has width of (3 d), where 3 the hyperparameter that defines the learning capacity of the network. Essentially, our network has two tunable hyperparameters, r and B. detailed study on these hyperparameters will be discussed in Appendix F
Since we require the network to be continuously differentiable with respect to both the spatial coordinates æ and the latent space vector q. we adopt continuously differentiable activation functions. In practice, either ELU (Clevert et al., 2015) or SIREN (Sitzmann et al. (020b) serves this purpose.
For the encoder network Co., the input vector of dimension r. Using MLPs would lead to enormous network size when P is large. We therefore design an encoder network first to apply multiple ID convolution layers kernel size stride size 4, and output channel size d until the output ID vector's length is the closest to 32/d /d but no smaller. Afterward, the encoder network reshapes the vector into I channel and applies an MLP layer to reduce the dimension of the vector to 32. The last MLP layer then transforms the previous 12-dimensional vector into dimension r. Discretization-speaiii e.g., PointNet for point clouds and convolutional neural network for grid data.
ID vector of length P with d channels. The output is
encoder networks can also be used,
D.1 TRAINING DETAILS
We the Adam optimizer (Kingma & Ba 2014) for stochastic gradient descent. We the Xavier initialization for ELU layers and the default initialization (wb 30) for SIREN layers. Unless otherwise noted, we train with base learning rate of Ir le 4 and adopt a learning rate decay strategy (10 lr + 5 Ir + 2 Ir + 1. Ir + 0.5 Ir + 0.2 Ir). For each aforementioned learning rate, we train for 30, 000 epochs. We adopt batch size 16 (i.e.. 16 simulation snapshots) for the encoder and therefore batch size of 16 P for the implicit-neus neurallrepressentaion based manifold- parameterization function. We implement the entire training pipeline in Py Lightning et al. 2019) which facilitates distributed training across multiple GPUs. The training vector fields are standardized have zero mean and unit variance. The training spatial coordinates are also standardized for the ELU implementation, but they are preprocessed to be between [-1, 1] for the SIREN implementation.
E NETWORK GRADIENTS
Gradients of the implicit neural representation (va9,, and öge, are crucial for our reduced- aq order pipeline. We use them to compute the spatial and temporal gradients of the vector fields,
öge,
vf V a (x, q) and f(x,t) = ōq (x, q)q. Furthermore, the gradient with respect to öge, is also key ingredient for the network inversion (Appendix Co.
the latent space 0q
E.1 NETWORK GRADIENTS VIA DIRECT DIFFERENTIATION
One way to compute the gradients is through directly differentiating the continuously differentiable network.
19

---

Published as a conference paper at ICLR 2023
O. 0.182740169
Relwwckk (ELUU mwan.
0.06446 1003
0.064905781
Cadeent Ento
85E:006
Sample
numerical (BIEEN)
Taboorkk (SIREN)
nummercaa (ELU)
Detwickk (ELU)
Figure 11: Gradient errors evaluated at the centers of tetrahedra. Numerical differentiation (specifi-
cally, finite element with linear basis) significantly improves the gradient accuracy of SIREN.
E.1.1
EFFICIENT IMPLEMENTATION
Algorithm 2: Network gradient
Input: a
8ge, duz
(x),, Io, (x) Let 90,- In.
Output:
91 be n-llaeer MLP network of interest.
Initialization:
dy
Bu
I. y
x
for
1...., n do By
0g. (y)
By
8m
(v)
end
By Ow
return
Unfortunately, computing these differentiations via auto-diff (computational graph tracking) is too slow for the nigh-performance application explored in this work. To ensure maximum efficiency, we implement the gradients analytically through the chain rules of each network layer, similar to the grad net approach by (Lindell et al. 2021). Since only use fully-connected layers, the implementation straightforward. In addition, along the way of computing the gradient via chain rules, we also obtain the function value of the neural network itself. Therefore, obtain both the gradient and the function value in a single forward pass. Algorithm 2 describss the details.
1.2 FAILURE MODES
While direct differentiation is generally accurate, we observe failure modes when the training samples sparse. Consider the where 9a, (x, q) æ, e.g., the undeformed, reference configuration of elastodynamics. In this case, the ground truth spatial gradient is V v9
We train the neural network with spatial samples discussed in Appendix L.4.1 After training, evaluate the spatial gradients at the center of each tetrahedron and compute its error (from the identity matrix) in the L2 norm. As shown in Figure []] errors observed in both ELU and SIREN networks while SIREN's s error is significantly larger, c.f., (Yang et al. 202T). This discrepancy be understood as the high-frequency prior by SIREN being more suitable for approximating
20

---

Published as a conference paper at ICLR 2023
0.00012
0.0001
Reconstruction Error (mean squared)
0.00008
0.00004
0.00002
25
35
65
MLP Layer Width
Figure 12:
Reconstruction accuracy vs. the MLP layer width (3). For each setup, we repeat the training for 8 times. The mean values and the error bars are shown. Increasing MLP layer width leads to higher accuracy.
high-frequency functions. By contrast, the low-frequency prior by ELU is better approximating low-frequency functions (Hertz et al. (02T).
We find training with more samples and/or gradient supervision (Chen et al. 2023) reduces the
above at the cost of additional training resources but does not resolve the issue completely. Future work may consider more advanced gradient regularization methods (Liu al. 2022
E.2 NETWORK GRADIENTS VIA NUMERICAL DIFFERENTIATION
Another approach to computing neural network gradients is numerical approximation, such as the finite difference method and the finite element method. Figure 11 demonstrates the drastically improved gradient accuracy with the finite element method (linear basis function) on networks trained with SIREN. We find such a hybrid approach be ideal middle ground, i.e., represent the function itself with the discrreizaaiionnnddeeeeeen neural network and compute the gradient via discretization-depal dependenn numerical differentiation.
With numerical gradients, users still benefit from our continuous representation that allows them use any discretization for gradient computation (c.g., ground impact elasticity deformation with remeshing). Prior ROM approaches require these numerical gradients to be computed by fixed mesh.
F HYPERPARAMETERS
During the offline training stage, we can modify two network architecture hyperparameters: the size of the latent space (r) and the width of the MLP (8). After training and during the online deployment stage of the model, the key hyperparameter is the number of integration samples (.M).. In this section,
sensitivity analysis these hyperparameters.
Figure 12 shows that increasing the width each MLP layer improves the reconstruction accuracy. Unless otherwise noticed, all reconstruction errors reported are mean squared errors (MSE). While Figure 12 demonstrates an empirical convergence, future work should consider theoretical analysis
(Kunisch & Volkwein 2002).
21

---

Published as a conference paper at ICLR 2023
Inference Cost
MLP Layer Width B
-+64X Intel(F)) Xeon(R) CPU E5-2686 2.30GHz
NVIDIA Tesla V100
Figure 13: Network inference cost (CPU and GPU) vs. the MLP layer width (3). Inference cost (over 100 runs) increases as the network size increases. CPU suffers from a more severe performance drop than GPU.
However, large layer widths also entail higher computational cost (see Figure 13 the highly parallelized GPUs, CPUs suffer from more severe performance drop due to the larger network size.
Compared to
Unlike traditional PCA-based model reduction approaches (e.g., POD), increasing the dimension of the latent space vector does not yield higher accuracy for our approach (see Figure (0]. In prior ROM approaches, augmenting the dimension of the latent space vector directly increases the number of entries in the manifold- paramttttttatmniitttionn function matrix, which, in turn, improves the approximation capacity of the manifold pariieerraaaamttrizattoooon function. However, as shown Figure [2] control the learning capacity of the network
in the implicit neural representation framework, architecture by explicitly modifying the width of the MLP layer. Therefore, the importance of the latent space dimension diminishes. This experimental observation is consistent with prior work's 's theoretical analysis (Remark 2.1 by Lee & Carlberg (2020) which says the lower bound the latent space dimension is the dimension of the problem parameters (e.g., material properties) plus one. Consequently, one should consider increasing the learning capacity of the canifold-paramettrriattion function instead of increasing the latent space dimension.
During the online stage of the reduced-order model, we select integration samples via the greedy sampling scheme (Appendix B). Increasing the number of integration samples leads higher accuracy (Figure 14p.
G
THERMODYNAMICS
G.1 CONTINUOUS PDE
Du Øt
82u O.
v(x)
(15)
In thermodynamics, we study the 1D heat equation of the spatiotemporal dependent temperature u. We assume zero-Dirichlet boundary condition, though other boundary conditions can also be incorporated.
G.2 FULL-ORDER MODEL
The full-order model discretizes the spatial vector field with regular Eulerian grid using 501 equally spaced samples (P 501).
22

---

Published as a conference paper at ICLR 2023
10
Enror (%)
0.01
0.001
0.0001
0.00001
0
100
200
300
400
500
Number of Samples
Training data
Testing data
Figure 14: PDE solution accuracy vs. the number of samples. Increasing integration sample counts yields higher accuracies on both training and testing data.
We then approximate the spatial gradient using the finite difference method,
2uu u (x" tn)
u(x² Ar, t.) + u(x* + Ar, t..) 2u(x", tx))
(16)
dr2
2
where Ai is the grid spacing.
We use this to compute the next-time step velocity,
d = v(x²) 02
(x tn).
(17)
dz
We assume first-order explicit time-stepping scheme,
" n/ Ati
(18)
G.3
REDUCED-ORDER MODEL.
Instead of a finite-difference treatment of the spatial gradient, we compute the gradient via direct differentiation of the network. The latent space dimension is r 16 and the width of the MLP B 128.
G.4 TRAINING AND TESTING DATA
We generate training data by setting v to different piecewise constant functions of three regions, i.e., the parameter vector μa D [0.2, 1.0 C R Figure 15 displays the initial and the final states of training data. Notice how different regions diffuse with different speeds due to the spatially varying v. In total, 8 PDE temporal sequences (of 100 time steps) with different v's generated for training. We sample another v's for testing purposes. All training and testing data adopt the same initial condition (Figure 15
G.5 RESULTS
Figure 16 demonstrates the performance of the reduced-order model the testing dataset. With the robust sampling scheme, we can achieve 0.43% етrоr using just 22 spatial samples (see Figure 4 In addition, by incorporating different initial conditions (e.g., Gaussian profiles centered at different horizontal locations) in the training data, our approach can capture diffusion with a wide range of initial conditions (Figure
23

---

Published as a conference paper at ICLR 2023
~
1.0
0.2
1.0
1.0
1.0
1.0
0.2
0.3
0.7
(a) Initial state
(b) Training case
(c) Training case
(d) Training case
Figure 15: Thermodynamics: training data.
(a) Testing case 1
(b) Testing case 2
(c) Testing case 3
(d) Testing case 4
Figure 16: Thermodynamics: reduced-order simulation on the testing dataset.
1.9
0.5
0.4
X
Figure 17: Our reduced-order diffusion solver with various initial conditions (transparent) and the final diffused states. Different colors correspond to different initial conditions.
24

---

Published as a conference paper at ICLR 2023
To gauge approach's 's ability to conserve physical laws, time the region between æ1 = 0.396 and 22 = 0.526 (Figure T6p. As shown by Cannon (1984), th the escaped energy (EE) is measured by the flux difference &e , t) 8e (22, t.) the left and right boundaries integrated over time; the stored energy (SE) is measured by the heat energy density, proportional to the temperature (u), integrated over space; the total energy (TE) the sum of the prior two energies, i.c. TE=EE+SE. Figure 4 demonstrates that the reduced simulation preserves the total energy over time.
the energy evolution
du
H
IMAGE PROCESSING
H.1 CONTINUOUS PDE
v
(19)
In image processing, we study the Perona Malik diffusion equation of the gray-scale image u (Perona)
& Malik 1990p.
H.2 FULL-ORDER MODEL
The full-order model discretizes the spatial 2D pixel field with a regular Eulerian grid using 256 x 256 equally spaced samples (P = 65, 536).
then approximate spatial gradient using the finite difference method,
8au Ox (x', v. ,tn) u(x² Ax, y², tn) + u(x² + Ar, y², tx) 2u(x', , ,t,tt)
(20)
Ar
Su ay² (x', (I',v', ,t,)= u(xi,y² Ay. tn) + ula", yf + Ay. fn) "u(x",,".
(21)
Ay²
We use these to compute the next-time step velocity,
d v(x', y²)( 82 Oz
Ou
(x', v',tna kn)
Oyz (x',8', tt,)))
(22)
We assume first-order explicit time-stepping scheme,
"h = u/, + Ali
(23)
where Ar = Ay is the grid spacing.
H.3 REDUCED-ORDER MODEL
The reduced-order model follows the same treatments as Appendix G.3
H.4 TRAINING AND TESTING DATA
Similar to Appendix G.4 we generate training data by setting / to different piecewise constant functions four squally-spaces regions, i.e., the parameter vector μ E D {0, 0.2)* 4 R Diraiu contains 11 parameter vectors (Figure 18), where 0, 2 or 3 regions have zero diffusion coefficients, ie., no blurring. Diest contains parameter vectors, where only region is unblurred. For each training and testing parameter vector, run a simulation of 20 time steps. All training and testing data adopt the same initial condition (Figure 18)).
H.5 RESULTS
Figure 19 demonstrates the reduced simulation's 's performance the testing dataset. Furthermore, with the robust sampling scheme, we can achieve 36.1 PSNR using just 63 spatial samples (see Figure 20p. The dimension reduction is 99.90%.
25

---

Published as a conference paper at ICLR 2023
2
2
(a) Initial state
(b) Training case 1
(c) Training case 2
(d) Training case 3
Figure 18: Image Processing dataset. The blurred region is numbered.
2
2
2
4
(a) Testing case
(b) Testing case 2
(c) Testing case 3
(d) Testing case 4
Figure 19: Image Processing: our reduced-order simulation successfully blurs the target regions in the testing dataset. The blurred region is numbered.
(a) Initial state
(b) Ground truth
(c) Optimal sampling (ours)
(d) Uniform sampling
Figure 20: Blurring of portions of the image (20 time steps). (b) Ground truth solution uses all
P = 65, 536 pixels. (c) Our approach uses very few integration samples (IM) 63, blue circles) and obtains a similar result as the ground truth (PSNR 36.1). (d) Naive uniform sampling (MM] = 64, blue circles) leads to poor agreement with the ground truth (PSNR 23.4).
26

---

Published as a conference paper at ICLR 2023
I ADVECTION EQUATION
Nonlinear ROM techniques have been shown to significantly outperform linear ROM techniques on problem with slowly decaying Kolmogorov n-widths (Peherstorfer 2022 Lee & Carlberg, 2020). In this section, we test CROM, also nonlinear ROM method, on transport-dominated problem, given by the advection equation.
1.1 CONTINUOUS PDE
au 0t + (a V u = 0,
(24)
where u is the advected quantity and is the advection velocity,
1.2 FULL-ORDER MODEL
The full-order model discretizes the ID spatial field with regular grid using 100 equally spaced samples. The PDE is temporally discretized via forward Euler.
1.3
REDUCED-ORDER MODEL
We vary the latent space dimension and fixed width of the MLP 8 20.
1.4 TRAINING AND TESTING DATA
In this example, consider reproductive case where the training and the testing data are the same. In particular, the initial condition is given by Gaussian profile and is then advected under a constant velocity.
I.5 RESULTS
Since this example is reproductive and the only parameter of the system is time f, the intrinsic solution-manifold dimension (Lee & Carlberg 2020) is 1. Figure [2] demonstrates that CROM indeed accurately captures the transport behavior of the Gaussian profile even at the the intrinsic colution-manifoll dimension (r 1). By contrast, POD suffers from serious artifacts with lower dimensional latent space and requires significantly higher dimensional latent space (r 16) in order model the advection behavior. Furthermore, Figure showss that our approach is also more accurate than the convolutional autoencoder approach by Lee & Carlberg (2000) when =
BURGERS' EQUATION
In this section, we continue to investigate problems with slowly decaying Kolmogorov n-widths. In particular, obtain the setup Lee & Carlberg 2020) (see Section 7.1 of their paper) and compare with their dnisxrrtizzaioo-deeeeddee convolutional autoencoder manifold.
J.1 CONTINUOUS PDE
Ow 30.5w²
(25)
Øt
dz
where w is the variable of interest. The initial condition is as follows.
(0 0) = 4.25
(26)
w(I 0) We further assume Neumann boundary condition.
Yx E (0, 100
(27)
27

---

Published as a conference paper at ICLR 2023
900 90
08 96
x=4
40
7=8
16
fruth
menifold)
7: space t:
Figure 21: Advection. Our CROM employs nonlinear manifold oarifodddprrammeerizaaiioo function and accurately captures the advection phenomenon with just one-dimensional latent space. By contrast, linear manifolds, such POD, require significantly larger latent space dimension. Lower dimen- sional latent spaces on linear manifolds yield significant artifacts due to the inability to model slowly decaying Kolmogorov n- n-widths problems.
28

---

Published as a conference paper at ICLR 2023
Intrinsic solution-manifold dimension
0.1
Relative Emro
0.01
0.001
7
13
16
Latent S pace Dimension
C O M (Ours)
- [Lee and Carlberg 2020]
Figure 22: Reconstruction errors. Our approach yields higher accuracies than the convolutional autoencoder approach by Lee & Carlberg (2020) across different latent space dimensions. Our method's accuracy advantage is particularly noticeable at the intrinsic solution-manifol dimension (see the vertical dashed line).
J.2 FULL-ORDER MODEL.
The full-order model discretizes the ID spatial field with a regular grid using 256 equally spaced samples. The PDE is temporally discretized via forward Fuler.
J.3 REDUCED-ORDER MODEL
We vary the latent space dimension and use fixed width the MLP 3 64.
J.4 TRAINING AND TESTING DATA
We generate training and testing data by varying the parameter jap. We set the training parameter to D.Div (0.015 + (0.015/7)))5=0..... resulting in training cases. For online testing,
the parameter to JDD..t = 0.021.
J.5 RESULTS
Since we can uniquely parameterize the family of training and testing data with two variables / and μDD the intrinsic solution-manifole dimension (Lee & Carlberg 2020 is p* 2, which the lower bound for the latent space dimension
Figure 6p reports the result for p = 2. While both POD and convolutional autoencoder (CAE) approach from Lee & Carlberg (2020) display clear artifacts, CROM agrees well with the ground truth. Such an agreement at the latent space dimension lower bound (i.e., r = p²) shows that CROM very close achieving the optimal performance of any nonlinear trial manifold.
Furthermore, Figure 2auure222eene demonstratess the accuracy advantage our approach over the CAE approach across different latent space dimensions.
29

---

Published as a conference paper at ICLR 2023
K FLUID DYNAMICS
K.1
CONTINUOUS PDE
uu V u Vp + V eet: V =0
(28)
We study the incompressible Navier-Stoker equations, where u is the velocity field. is the pressure field, and is the viscosity. We assume there is external force, i.e., Sent O. We assume no-penetration boundary condition for the velocity and a zero-Dirichlet boundary condition for the pressure.
K.2 FULL-ORDER MODEL
The full-order model assumes an operator-spliting scheme following the classic Chorin's projection method (Chorin 1968). In particular, follow the implementation by Stam 1999). We summarize the major ingredients and refer to his paper for details.
We sequentially apply 3 linear operators the velocity field: diffusion, advection, and projection. From the previous time step velocity field un, we apply the diffusion operator according to the viscosity and obtain the diffused velocity u Next, we employ sami-Lagrnngian to obtain the advected velocity field u.... For the last projection step, we first compute the pressure Pn+1 by solving the poisson equation, V2pn+1 obtain the divergence-free free velocity un+1- All operations are done on 2D 16 by 16 Eulerian grid.
u. Afterwards, we apply the pressure gradient to
K.3 REDUCED-ORDE MODEL
We compute latent space dynamics (4n+1) of the velocity field by projecting un+1 onto the manifold- parameterization function 98,- For network hyperparameters, we adopt latent space dimension
r = 6 and a MI.P width | = 10.
K.4 TRAINING AND TESTING DATA
We generate training data by varying the spatially-consten viscosity: μ E D 0.00. R. We sample 3 different viscosities from D. For each viscosity v. we run a simulation of 50 time steps. In total, 150 time steps are used for training. For testing, first sample 3 viscosity values from D [0, 0.02). Furthermore, we also sample 3 values from the extrapolated region Dextrapolate [0.02. 0.1]. Notice how the extrapolated region has larger range (0.08) than the interpolated region (0.02). All training and testing data adopt the same initial condition (see Figure 23
K.S RESULTS
Figure 23 demonstrates method's ability qualitatively capture fluid dynamics with different viscosities on both interpolated data as well as extrapolated data. Surrogate models that predict velocity fields directly from viscosities (e.g., using a neural network) typically cannot handle extrapo- lation. By contrast, our approach handles extrapolated scenarios because we only approximate the spatial representation (kinematics) with neural network, not the PDE itself. Therefore, as long as the extrapolated cases involve vector fields that can be represented by the low-dimensional implicit neural representation, our reduced-order solver is able evolve the latent space vector capture them according the PDE. However, if the extrapolated quantity lies outside of the space spanned by the manifold-parameterre CROM, like most ROM approaches, does not handle arbitrary extrapolation scenarios. See more discussion on generalizability at the end of the main text.
function, our approach would not be able to capture the dynamics.
To further test the proposed method, we set up Karman vortex street experiment. The Karman vortex street is particularly challenging as the flow bifurcates from a steady state to the periodic vortex shedding regime. Mathematically speaking, the fixed point solution undergoes Hopf bifurcation and becomes a limit cycle (Huerre & Monkewitz/ 1990). To capture such an effect, the latent space dynamics needs to be robust under extreme nonlinearity. Many neural-network- methods, such
30

---

We follow the geometry and material parameters outlined by Nabizadeh al. (2022). No hyper- reduction considered. Figure shows that our latent space dynamics method faithfully represents the ground truth flow: from the initial steady state to the bifurcated vortex shedding regime. The ground truth solver utilizes a 100 by 200 grid, with a total 40, 000 degrees of freedom. By contrast, our method captures the same level of details while employing a latent space dimension of only nature, our approach also generalize beyond Published as a conference paper at ICLR 2023
Interpolation
0.0000
0.0075
0.0125
Initial condition
16
Extrapolation
0.0200
0.0600
V 0.1000
Fluid dynamics: our method captures the velocity fields of a wide range of viscous fluids. Large viscosity severely dissipates the velocity. All frames shown are captured from time step 25. The first row consists of the interpolated result while the second row consists of the extrapolated result well outside the training range.
Figure 23:
physics informed neural network (Raissi al. 2019). well-known fail in this case (Chuang & Barba 2022).
r = 4. Due to CROM's ainematies-pppooximttionn the training temporal range.
Future work should test CROM on more challenging boundary conditions (Copeland et al. 2022 and chaotic systems involving Markovaoo effects (Pan & Duraisamy, 2018).
L SOLID MECHANICS
L.1 CONTINUOUS PDE
We solve the clastodynamics equation arisen from solid mechanics for deformable soft body modeling:
P(V + B.
(29)
where ((X, maps the undeformed (reference) position (X) of an arbitrary material point from the reference configuration to its deformed (current) position (x) at time Proper Dirichlet and Neumann boundary conditions are applied.
We assume the first FRla-KKirhhhe stress is strictly a function of the deformation gradient F VΦ. In particular, adopt the corotated linear elasticity constitutive law,
P(F) 2F(F - R. Atr ( F I)R,
(30)
where the rotation tensor from the polar decomposition of the deformation gradient F R.S. and μ are the first lame parameter (Pa) and the second lame parameter shear modulus (Pa). The proposed framework also works with other elasticity constitutive laws and can be easily extended to capture hysteresis and inelasticity.
31

---

Published as a conference paper at ICLR 2023
L.2 FULL-ORDER MODEL
We adopt the retrahedral-mest bassh linear finite element method (FEM) for the full-order model. We closely follow the course by Sifakis & Barbic (2012). Here we review only the salient features and refer to their work for theoretical and practical details.
We spatially discretize the domain of interest using linear tetrahedra. For each tetrahedron, we have four undeformed vertex positions (X I, I. X 4) and four deformed vertex positions (0),' 244). For each vertex, we have x; FX; + b, where is the translation.
Consequently, we can obtain the deformation gradient F of each tetrahedron by solving D. = FD=
where D2 = [ 324 $2 2'4 3333 24 and Dm = =X - X X² Xx X3- Xa].
Furthermore, the internal elastic forces the vertices (f1 (f1.... .fa) can be computed though det m. (mF)) -T and f. f: f2 f3.
[f1 f2 fzl det
L.3 REDUCED-ORDER MODEL
Since the initial condition known for a given problem, we opt to construct manifold- parameterization function for the displacement field uXX,t) Φ X. The deformation map can then be computed via Φ = u + X. The hyperparameters of the manifold-parammettiiziiii function are and 20.
Unlike previously discussed PDEs, the elastodynamics equation second-order in time. Conse- quently, in addition to the function value of the manifold, also need to infer the temporal derivative during latent space dynamics. We do so by employing the tangent space of the manifold,
0ga= q. 0q
ल(x, tn)
(31)
where 9n = (q., - 1n and 90 = 0.
To compute the acceleration (and increment the velocity) of vertex belonging to the integration samples set M, we need to evaluate the internal force at this particular vertex. We can calculate the vertex force by accumulating the internal forces of tetrahedra incident on this particular vertex. Tetrahedron forces are computed using the formula described in the previous section. To attain optimal gradient accuracy (Appendix E.2). we opt use finite element's linear basis function to compute the deformation gradient F of each tetrahedron.
To compute these internal forces and deformation gradients, we also need to obtain the position information (via the manifold-parre parameterization function) of all vertices in the one-ring neighborhoods of the integration samples. We refer to the set of these neighboring vertices as N. Therefore, the total number of spatial samples in the reduced-order model now becomes [M] + [W]. Since one-ring neighbors only entail a small subset of the original degrees freedom, we still offer significant spatial sample reduction from the full order model (P).
L.4 TRAINING AND TESTING DATA
L.4.1 GRAVITY-INDUCED IMPACT
An initially-static rectangular-shang shaped deformable object (P = 2, 065 vertices and 9, 346 tetrahedra) accelerates due to downward gravity. After impacting the collision objects underneath, the material undergoes intense deformation.
We generate training and testing data by sampling the shear modulus μ = [60, (00, 70, 000] C R.. The first lame parameter À is fixed to be zero, in which case the coro- tated linear elasticity constitutive law simplifies to the As-Rigid-As --sPoiblAsee energy (Sorkine & Alexa 2007). Specifically,
D
uniformly generate 4 samples of the shear modulus for training data. We simulate 200 time steps for each shear modulus. Therefore, total of 800 time steps are used for training. We then randomly generate 3 samples of the shear modulus from D. We also simulate 200 time steps for these shear moduli and use them for testing. We use the same initial condition for training and testing.
32

---

Published as a conference paper at ICLR 2023
6200, 46667; 1.28% error
# 10800, 10000; 0.91% error
# = 15400, = 65000; 2.24% eIrOr
= 20000, = 28333; 2.79% error
Figure 24: Our approach captures a wide range elastic moduli with modest errors compared to the full-order ground truth.
L.4.2 TORSION AND TENSION
We study another common solid mechanics application where the material undergoes torsion and tension. The left and the right boundaries kinematically move constant translational and angular velocities of the same magnitude but in opposite directions.
We solve the parameterized family of problem where μ (H, E D [6. 200, 20, 000] x [10, 000, 65, 000] C R2. We obtain 9 training samples via uniform full-factorial sampling of the parameter space. We simulate 100 time steps for each training sample. We then use Latin hypercube method to randomly generate another 4 testing samples. We also simulate 100 time steps for these testing samples. Training and testing data adopt the same initial condition.
L.5 RESULTS
Figure demonstrates testing (μ 62500, 2.68% emror from the full-order ground truth) for the experiment from Appendix [.4.1 Since our method has ilisereizztiionagnnsst architecture, we can adapt the simulation resolution improve efficiency. The adaptive simulation adopts contact-based oracle that switches from the low-resolution mesh to the high-resolution mesh once contact event is detected. Before the start the simulation, we generated both the low-resolution and the high-resolution meshes using TetWild (Hu et al. 2018) with small and large ideal edge lengths.
After training on the experiment from Appendix L.4.2] our method captures a wide range of shearing and volume-preserrvng behaviors (see Figure
Furthermore, to showcase CROM's ability to handle large deformations, we reproduce the dinosaur example Syen et al. (2021). Figure[5]shoove that CROM accurately captures the large deformations of the dinosaur. By contrast, under the same latent space dimension, POD has a limited expressivity and suffers from visual artifacts: the head and the tail's motions are constrained. This observation is consistent with the results reported by Shen et al. (2021) where nonlinear approaches outperform linear approaches. In terms of memory consumptions, CROM is offlers---mmannuuude more efficient than these prior discretization -eepiiiatiion-dddeentteee ROM approaches (Barbič 2012 Fuulon et al. 2019 Shen
al. 2021
33

---

Published as a conference paper at ICLR 2023
11X memory usage than ours
32X memory usage than ours
542X memory usage than ours
Ground truth
POD
Autoencoder [Fulton et al. 2019]
Autoencoder [Shen et al. 2021]
Ours
[Barbic 2012]
Figure 25: CROM agrees well with the ground truth simulation while consuming far less memory than prior discretization-deeendde network-bass based autoencoder approaches (Fulton et al., 2019 Shen et al. 2021). The experimental setup is taken from Shen et al. (2021).
ROM approaches, including POD (Barbič (0022) and neural-
M COMPARISON WITH PRIOR APPROACHES
We compare our approach with prior reduced-order approaches that generate fast PDE solutions via latent space dynamics. In total, we compare with four baselines: (1) the classic widely-adopte
proper orthogonal decomposition (POD) method (Berkooz al.) 1993; Holmes et al., 2012). Our POD implementation follows the PCA approach by Barbič & James (2005): Barbic 2012); (2) the
neural-network- network-basedd autoencoder approach by Fulton et al. (2019); (3) the neural-network-hasee
autoencoder approach by [Shen et al. (202T); (4) the convolutional autoencoder (CAE) approach by
Lee & Carlberg (2020). All these four baselines are dbscretization-deeeeddeee
SIREN, PINN, or GNN are related to our work but are very different in nature. While their goal is to directly solve PDEs with neural network architectures, our goal is to leverage neural networks for dimension reduction. Therefore, we do not compare with them and focus on comparison with other dimension reduction methods.
POD approximates the discretized vector field with a linear basis U,
7(1) (f(x', (f((x")).... .).... f( (), - f(a .D))² Uq(1),
(32)
where U is matrix of size Pd U is most easily constructed via a singular value decomposition of the discretized training data (Sifakis Barbic 2012).
Autoencoders have the same input and output POD. However, the linear basis is replaced with nonlinear neural networks,
7 ()( f(xt, ....)....a... ,0) h(a()),
(33)
where h is the decoder portion of the autoencoder structure. Fulton et al. (219 Shen et al. 2021): Lee & Carlberg (2020) all follow this formulation but differ their particular architecture choices.
We implement both POD and the autoencoder approaches the same Py framework as CROM. These approaches share exactly the same PDE time-stepping stepping algorithm (step 2) during latent space dynamics. For POD and the autoencoder approach, the network inference step (step 1) and the network inversion step (step 3) employ the linear basis U and the nonlinear decoder function h, respectively.
M.1 SPEED
We measure the wall-clock computation time during the falling stage (before contact) of the gravity- induced impact experiment (see Appendix Ind Figure In particular, since all model reduction methods (including our approach, POD, and the autoencoder approaches) share the same PDE time- stepping algorithm (step 2), we measure its wall-clock computation cost to ensure the generalizability of the comparison.
Since POD and the autoencoder approach are baked into particular discretization, it does not allow for adaptive discretization. Therefore, the high-resolution mesh (2, 065 vertices, 9, 346 tetrahedra) has to be used throughout the falling stage even though there is minimal deformation. Consequently, its computation cost is significantly higher than our approach, which employs a low-resolution mesh (8 vertices, 5 tetrahedra) during the falling stage. Figure shows that approach is faster
34

---

Published as a conference paper at ICLR 2023
SIMULATION EFFICIENCY: OURS VS. PRIOR CROM (Ours)
SIMULATION EFFICIENCY: OURS VS. PRIOR ROM
SIMULATION EFFICIENCY: OURS VS. PRIOR
=CROM (Oura)
CROM (Ours)
Prirr ROM
667
5.923
Wall (8)
Wall clock (8)
Wall clock (s).
00:377
0:136
0:156
1 CPU thread
4 CPU threads
32 CPU threads
Figure 26: Wall clock time comparison with prior ROM approach. Our approach obtains considerable
speedups across different computational resources.
than the prior ROM approach with different CPU thread counts, Our method's advantage is most obvious when there is extremely limited computational resource (1 CPU thread). Since this speedup comes from the discretization independent nature CROM, we expect it to generalize to other prior discretization-
ROM approaches as well. However, for application where remeshing is unnecessary, we do not expect this kind of speedup.
Hyper-reduction techniques (An al. 2008), such as the integration samples introduced in approach, can also be employed to speed up the high-resolution mesh solution by reducing spatial samples. However, this is orthogonal to our adaptive discretization contribution, and the hyper- reduction techniques can be employed to further accelerate the low-resolution mesh solution as well.
M.2 MEMORY
We compare the memory consumptions of our method, POD, and the autoencoder approaches with the same latent space dimension (r), both trained on data from the torsion and tension experiment (see Appendix[L.4.2 For all implicit neural representation networks, we set the MLP width to be 3 = 20. Independent of latent space dimension (r), network has near-constant memory consumption (see Figure TOP. By contrast, the memory consumptions scale linearly with More importantly, since the output dimension of POD and the autoencoder is Pd, their memory consumptions scale linearly with the number of discretized positions P. By contrast, our approach's memory consumption independent of the number of spatial samples. Consequently. in this large-scall scale example that features P = 66, 608 vertices, we observe more than ten-fold advantages with our method across all latent space dimensions (see Figure TO).
POD and the autoencoder approaches
In addition, we also compare our method with POD on the image processing experiment (see Appendix H). CROM uses an order-of-magnn -aagnit-- less memory than POD (see Figure (5t). In Figure 27 we also demonstrate CROM's memory advantage over the convolutional autoencoder approach by Lee & Carlberg 2220)).
Even though the overall memory consumption of our approach is significantly lower, our approach's 's per-sample evaluation cost is higher. For any given discretized spatial sample, POD only requires a vector-mas mittor--attixxxoo multiplication of d by matrix. By contrast, our approach needs to go through the entire MLP that employs d3 by d3 matrices. The autoencoder approaches share similar computation overheads.
M.3 ACCURACY
Figure 10 demonstrates the accuracy the networks studied in the previous memory section. Our method consistently offers onlers-o--magiiude smaller training accuracies. This accuracy advantage especially noticeable with small latent space dimensions (e.g.. r = 2). Further increasing the latent space dimension leads to higher accuracies with POD and the autoencoder approaches but not our approach. As discussed Appendix F our approach's 'ss accuracy depends more saliently on the MLP width as opposed to the latent space dimension. Furthermore, unlike these discreizaation-deenddee approaches, our continuous manifold manifold-paramettee parameterization function also facilitates accurate handling of adaptive spatiotemporal data (Pan et al. 2023).
35

---

Published as a conference paper at ICLR 2023
1.00E +06
Memory (Number of Parameters
1.00E +05
+04
3
10
Latent pace Dimension
CROM (Ours)
[Lee and Carlberg 2020]
Figure 27: Memory comparison between our CROM and Lee Carlberg 2020
POD
Ours
[Fulton et al. 2019]
[Shen et al. 2021]
40
35
30.1
Training Cost (min)
30
25
20
15
14.1
9.6
10
0.1
0
Figure 28: Training time for the elasticity example presented in Figure significantly less time than neural-network-baeee methods, including Fulton et al. 2019 2021) and ours.
25
POD Barbic
2012) uses Shen et al.
Furthermore, 5 shows that method is more accurate than POD the image processing experiment, both visually and quantitatively. Figure 22 displays our accuracy advantage over CAE on Burgers' equations.
M.4 TRAINING TIME
Figure 28 reports the offline training time for different ROM approaches. All neuural-ettorr-baare
methods, including ours, take longer time to train than the linear POD approach. We emphasize that the timing data is purely informative and calculated post-hod. The main target of our work is not
36

---

Published as a conference paper at ICLR 2023
SIMULATION EFFICIENCY: REDUCED VS. FULL Full
SIMULATION EFFICIENCY: REDUCED VS. FULL
SIMUL EFFICIENCY: REDUCED VS. FULL FFll 24.086
a cCRM (Ours)
=CRMM (Ours)
aCROM (Ours)
139.80
Wail clock (s)
26.631
Wall clock (9)
Wall clock
1:008
00:099
0:585
64X Intel(R) Reon CPU ES- 2686 v4 @ 2.30GHz
NVIDIA GeForce RTX 3080
NVIDIA Tesla V100
Figure 29: Wall clock time comparison with the full-order method. Our approach obtains significant speedups across different computing platforms.
COMPUTATION TIME BREAKDOWN OF THE REDUCED MODEL 5.05%
18.18%
76.77%
Network inference (step 1)
PDE time-stepping (step 2)
Network inversion (step 3)
Figure 30: CROM timing breakdown. Neural network operations introduce minimal overhead. The computation bottleneck remains to be PDE time-stepping.
optimize training time but to achieve efficient latent space dynamics. Nevertheless, specifically for implicit neural representation approach, future work may consider expediting training time via more advanced data structures and optimization (Müller et al. 2022: Liu et al.) 2020 Martel et al.
1 Takikawa et al. 2021).
N COMPARISON WITH THE FULL-ORDER MODEL.
To gauge the practical performance of our approach over the (unreduced) full-order model, we measure the sheer wall-clock performances of our reduced-order model and the full-order model.
To facilitate fair comparison, we implement, optimize, and (fully) parallelize the full-order PDE solution and the reduced-order approach within the PyTorch framework (Paszke et al. 2019) without any other dependency. The full-order model and the reduced model share the same (parallelized) PDE solver code (step 2), and the only difference between them is the neural network evaluation in the reduced model. To encourage fair comparison, for operations involving neural networks (step I and step 3), do not use customized SIMD vectorizations, CUDA kernels, or highly optimized inference libraries such as TensorRT (Vanholder 2016). though these optimizations should lead to further performance gain of algorithm. In addition, having both models implemented in the same framework allows us to compare them under the same computing environments, avoiding biased comparison where the full-order model runs on significantly limited computing resource (e.g.. 8-core CPU) while the reduced-order model runs on nigh-performance computing device (e.g.. 5120-core GPU).
Figure 29 reports the wall clock time of our reduced-order method and the full-order method across different computing platforms. While our approach is significantly faster on all computing platforms, its strongest speedup is obtained on the CPUs. This is unsurprising since the computationally expensive full-order model benefits more from high-end data-center GPUs (e.g.. V100) than the reduced-order model. In addition, because the neural network in the reduced-order model only employs a few 60 by 60 matrices (d = 3, 3 = 20). inference on CPU is also efficient. This opens doors for employing our models limited computing platforms such mobile and VR devices.
37

---

Published as a conference paper at ICLR 2023
100.00%
50.00%
12.50%
Enor
8.29%
3.13%
1.58%
0.78%
0.5
1.5
2 Timing
2.5
3.5
8-5, LLS, GPUs, IM-10 BB-,, LLS, GPUs, IMI40 4B-5, NLS. GPUs, IMi-10 xB-6, NLS, GPUs, IMI-40
B-20, LLS, GPUs, IMi-00 -8=20, LLS, GPUs, IMI-00 =B=20, NLS, GPUs, IMI=10 B-20, NLS, GPUs, IMI+00
8-5, LLS, CPUs, IMI10 48-5, LLS, CPUs, IM-40 4B-5 NLS. CPUs, IMI10 xB-6, NLS, CPUs, IMI400
LLS, CPUs, IM=10
CPUs, IMI=40
CPUs, IMI110
CPUs, IMI=40
Figure 31: Accuracy vs. Speed: our approach under different setups. In comparison, the ground truth full-order takes 139.855s on CPUs and 24. 24.086s on GPUs (see Figure 29).. Our approach obtains significant speedups with less than 5% error under various settings.
Figure 30 further breaks down the computation time spent on each component the reduced-order model. While the overhead of the neural network operations (step I and step is non-negligibl,, the bottleneck of the algorithm remains to be PDE time-stepping (step 2). In particular, the network inversion cost is very low, thanks to the linearized inversion solver discussed in Appendix C
Figure 31 reports the accuracy and the speed of our approach using different MLP layer widths (8). neural network inversion methods (nonlinear least squares (NLS) vs. linear least squares (LLS)). computing platforms (CPU vs. GPU), and integration sample sizes ([M)). In general, we observe similar accuracies between NL.S and L.SS (e.g.. [ 20, L.LS, GPUs, | = 40] vs. [3 20, NL.S, GPUs, |M| 40]), while LLS is significantly faster. Increasing sample size lead to higher accuracies but also require longer computation time (e.g., [3 = 5, NLS, CPUs, [M| = 10] vs. [3 = 5. NLS, CPUs, [M| 40]). Increasing the MLP layer width yields higher accuracies but also increases the computation time. In particular, the computation time increased more on CPUs (e.g.. [ = 5, LLS,
CPUs, M 40] vs. [5 20, LLS, CPUs, [M| 40]) than on GPUs.
O LATENT-SPACE TRAJECTORY
A key element of CROM is latent space dynamics, in which the latent space vector evolves nonlinearly over time under explicit PDE constraints.
Figure 32 demonstrates nonlinear latent space trajectory from the image smoothing experiment. Guided by the PDE, the latent space vector traverses through the latent space and obtains the desirable smoothing effects in the image (see Figure 33
38

---

Published as a conference paper at ICLR 2023
0.008
Principal Component
0.006
0.004
320,.
320
.200
20
0.000
80
40
-0.002
Latent trajectory
-0.004
Latent interpolation extrapolation
-0.02
0.00
0.02
0.04
Principal Component 2
Figure 32: Latent-space trajectory. We visualize the first two principal components of the latent space vector. The red numbers displayed on top of the trajectory are time step numbers. The simulation starts from 0 and ends at 320. The corresponding images of these time steps are shown in Figure 33 Overall, we observe smooth nonlinear trajectory over time. In addition, we also linearly interpolate and extrapolate inside the latent space. The corresponding images are displayed in Figure 34
36
t=
t= 107
t= 142
t= 178
213
t 249
== 284
t 320
Reduce-order simulation of image smoothing at different time steps t. The corresponding latent space trajectory visualized Figure
Figure 33:
After the simulation is finished, a postprocessing experiment, we linearly interpolate and extrapolate inside the latent space (see Figure 32 The interpolation end points are the latent vector for the unmodified image and the latent vector for the final image determined by the latent space dynamics. Figure 34 demonstrates the corresponding smooth transitions in the image space.
P STABILITY ANALYSIS
systtmmttiaalyy studies the stability of our method. CROM exhibits much larger stable time step sizes than unreduced counterpart. This is consistent with the theoretical analysis reported in the literature (Bach et al. 2018
Figure 35
extrapolation
t=
interpolation
t= 320
extrapolation
Figure 34: Linear interpolation and extrapolation of the latent space vector. We observe smooth
transitions. The corresponding latent space vector values shown in Figure 32
39

---

Published as a conference paper at ICLR 2023
10a
104
Overall Error 196)
10a
10-4
Reducedi order Model (Ours)
10 2
16
64
Time Step Size (Multipliers of the Base Time Step Size)
Figure 35: Stability analysis. We systematically increase the time step size (At). Our method remains stable in a wide range of time step size and maintains an error close to 1%. By contrast, the full-order model (that was used for training data generation) goes unstable the time step increases and suffers from major visual artifacts.
Q LATENT-SPACE DYNAMICS PSEUDOCODE
Algorithm 3: Latent-space Dynamics
Input: Latent space vector 9T,
Output: Next-time-step latent space vector q,+1
Network inference: gather full-space information (Algorithm 4)
PDE time-stepping with boundary conditions: update full-space information via differential equations, e.g., the elastodynamics equation (more details available from Appendix G to Appendix D.)
Network inversion: optimally project the updated full-space information onto the low-dimensional embedding (Algorithm 5).
Algorithm 4: Network inference: gather full-space information
Input: Latent space vector Qn.
Output:
Full space information f f.
for E M do
Compute the vector field value: fly tn.) = 90, q.).
8ge, oq In.
Compute the temporal gradient either by differentiating the network f(y, (n)
where q.. (q., 4 (n-1)/At. or by numerical approximation (Appendix Compute the spatial gradient either by differentiating the network, Vf(y. 6n) V "y9,, or by numerical approximation (Appendix
end
Algorithm 5: Network inversion: optimally project the updated full-space information onto the low-dimensional embedding
Input: Updated full space information fa+1
Output: Updated latent space vector qn+11
Find the latent space vector an+1 that best matches the evolved configuration fn+1 by solving
the least- squares:
min 4o+t ERr
Σ DEM
199, (n41) fly, (1,,4411113
40

---

Published as a conference paper at ICLR 2023
PDE
Full-order degrees of freedom
Latent space dimension
Dimension reduction
Spatial sample count
Spatial sample reduction
ID Thermo
501
16
31x
22
23x
2D Image
65,536
16
4,096x
63
1,040
ID Advection
100
100x
n/a
n/a
ID Burger
256
2
128x
n/a
n/a
2D fluid
512
6
86x
n/a
n/a
3D solid impact
6,195
3.098
n/a
n/a
3D solid piggy
199,824
99,912x
40
1,665
3D solid dinosaur
11,565
1
11,565x
n/a
n/a
Table 1: Dimension reduction statistics of our reduced-order model
R
REDUCED-ORDER MODEL'S STATISTICS
CROM obtains considerable dimension reductions across different PDEs (Table
41


Reviews:
[{"review_id": "wf0HTCiCvM", "reviewer": "Reviewer_YNJM", "summary": "", "questions": "", "limitations": "", "rating": 8, "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": "", "summary_of_the_paper": "The CROM approach proposes to solve the driving high dimensional PDE in a low-dimensional manifold, without having to discretize the continuous high-dimensional vector field. The discretized high dimensional spatio-temporal solution is projected in the low dimensional manifold through an encoder. The encoder produces a latent representation of the continuous vector field. A positional-dependent decoder is then employed to reconstruct the continuous-vector field. In order to evolve the PDE in time few saples of the continuous vector field are required to minimize the residual of the driving PDE with a PINNs-like loss. The limited number of spatial samples y and the low dimension of the latent representation allow to boost the solution of the PDE with a considerable reduction of computational resources (e.g. time and memory) without degrading the accuracy of the solution.\n\nSeveral technical solutions are also proposed to effectively solve the dynamics of latent space (e.g. greedy sapling of the spatial points, numerical approximation of the gradients, solution of the least squares problem for the network-inversion).\n\n ", "summary_of_the_review": "The proposed approach is novel and interesting. The paper might have a coniderable impact in machine learning and numerical simulation community. For this reason I recommand acceptance of the paper. I suggest to the authors to be clearer on the limitations of the approach (i.e. discretized field needed to initialize the latent state, no guaranties on the smoothness of the latent space at least in the proposed version). ", "strength_and_weaknesses": "The approach can be seen as ROM-PINNs method, as also asserted by the authors. Strength And Weaknesses of the approach are then related to the PINNs approach plus the benefits introduced by the reduced representation of the manifold. The proposed approach can effectively accelerate the numerical calculation of PDE without being agnostic to the physics. The same tool can also be used to dinamically remesh the high dimensionla solution. Nevertheless, despite \"The ROM community has demonstrated that these three steps can yield strong long-time stability even on stiff and chaotic dynamical systems\" (cit.) the PINNs are well known to fail with chaotic systems. Therefore, it would have been interesting to understand whether or not model reduction is of benefit to the solution of chaotic dynamical systems also when a PINNs-like criteria is used to evolve the solution in time.\n\nDespite the evolution of the solution in time is mesh-free, the whole method is not. In the Network Inference Step the latent vector at time zero is assumed to be given. Since the latent vector is the result of an encoder network, a discretized vector field (or at least a point-cloud) is always necessary. For PINNs a discretized vector field is also necessary in the training stage to enforce the initial condition, but this is not needed in inference.\n\nThe main limitation of the approach concerns the projection of the high dimensional solution into a low dimensional space. Authors claim that the low dimensional manifold is smooth. The smoothness is proven empirically for a rather simple test case but is not enforced in the training stage. Any kind of regularization is introduced for the latent state. This might be a major issue when dealing with physical systems that can undergo transitions of bifurcations. I suspect that the experiments in the paper are rather sipmle to experience such issue.", "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": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works."}, "extra_text": {}}, {"review_id": "rHJmhX3S-f", "reviewer": "Reviewer_7uQY", "summary": "", "questions": "", "limitations": "", "rating": 8, "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": "", "summary_of_the_paper": "This paper proposes to use implicit neural representations to model and efficiently solve PDE dynamics. This allows to solve dynamics in the latent space instead of the whole parameter space (lower dimension), without discretizing it. The latent space and network are trained on a range of sample PDEs for a given task, and employed to solve PDEs with different parameters.\nSolving the temporal evolution is still done with an explicit PDE constraint at each time step, providing updates of latent space parameters, instead of the explicit parameters which have much higher dimensions.", "summary_of_the_review": "Clear motivation, elegant method and strong experimental section.", "strength_and_weaknesses": "Using neural networks to reduce the dimensionality of physics phenomenon in a learned and non linear but smooth manner is a very good motivation, and as noted by the authors it had already been explored. Here it is also combined (for the first time) with an IMPLICIT representation, which allows for quantization-free modelling.\nThe latter point has two benefits: the pipeline requires less memory, and can easily adapt to any resolution by simply querying the network more or less densely.\n\nThe experimental section is very rich, displaying the applicability of the method for many applications (image blurring, heat equation, advection, elastodynamic deformation, …). Comparing network outputs to a full fledged PDE solver (slow) and other network based solvers quantitatively and qualitatively favors the proposed pipeline: it is almost on par with the full PDE solver, and largely more accurate than other network based approaches.\nA runtime and memory analysis also demontrates that the proposed approach is much more efficient once trained.\n\nMy main doubt is that such a pipeline could only (locally) interpolate between seen behaviors from the training set. How does this prevent applicability? How much care must be taken in crafting a training set that spans a sufficiently large spectrum of possible dynamics?\n\nAnother concern is about runtime comparison with full order models: for a fair comparison, maybe the training set generation and network training time would have to be accounted for in some way. Indeed, the trained approach can be seen as a way to smartly “memorize” simulation results to do much faster predictions. But this has an upfront computational cost, which only gets amortized when the resulting network is used for many different simulations.", "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": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "iTkzKm3Z4n", "reviewer": "Reviewer_NfXC", "summary": "", "questions": "", "limitations": "", "rating": 8, "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": "", "summary_of_the_paper": "\nThis work proposes a novel reduced order modeling approach based on deep neural networks. Here, the solution of the problem is approximated by a neural network, rather than a basis of mesh elements as is done in most reduction approaches. Dimensionality reduction is conducted through a *context* vector which is non-linearly encoded from an initial full order solution. Despite using a non linear reduction scheme, the approach manages to use the original model equations to advance the reduced state in time, bypassing the need for the development of a surrogate forecasting system.", "summary_of_the_review": "The paper proposes a method able to combine non-linear dimensionality reduction and intrusive simulation. Despite the need for more complex test cases, the paper clearly demonstrates the ability of the approach to efficiently model dynamical problems.", "strength_and_weaknesses": "\n**Strengths :**\n\n* Combination of intrusive simulation based on a *PINNs-like* scheme and non-linear dimensionality reduction which is a key challenge in reduced order modeling.\n* \"Mesh free\" approach : Although a discretised solution is still necessary for the initialisation of the CROM, the model is able to predict the solution at arbitrary inference points during the simulation.\n* The method outperforms classical reduction methods on the proposed test cases\n\n**Weaknesses :**\n\n* The main weakness of the paper are the test cases which might be a bit too simple. 1D diffusion problems are used to showcase the method and the case of the burger's equation clearly presents the interest of using a non-linear dimensionality reduction method. However, the case of the Navier-Stokes equations could be extended to more complex conditions such as the cylinder flow, or the fluidic pinball to demonstrate the performance of the approach on a more realistic test case.\n* Tying back to the previous point, issues such as non-markovian effects in the dynamics might appear when dealing with the reduction of complex chaotic problems, which are not treated in this work.\n* The name of the method should be changed, as a paper on cluster based reduced order modeling ([1]) already uses the acronym.\n\n*[1] E. Kaiser, B. R. Noack, L. Cordier, A. Spohn, M. Segond, M. Abel, G.Daviller, J. Östh, S. Krajnović and R. K. Niven (2014). Cluster-based reduced-order modelling of a mixing layer. [*Journal of Fluid Mechanics* , **754**](http://journals.cambridge.org/action/displayIssue?jid=FLM&volumeId=754&seriesId=0&issueId=-1), pp 365-414.*", "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": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "R6IKdEUfWZ", "reviewer": "Reviewer_XNfR", "summary": "", "questions": "", "limitations": "", "rating": 6, "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": "", "summary_of_the_paper": "The paper proposes to use a reduced order model (ROM) to accelerate PDE resolution.\nThey differ from prior ROM approaches by directly modeling the continuous function of interest, and not an already discretized function.\nTo approximate the continuous PDE solution, CROM relies on a latent space and use implicit neural representations as decoder, which enable them to address any discretization. The PDE is solved in a lower dimensional space, enabling them lower memory consumption. Experiments are conducted on several PDE depicting various discretizations, and comparisons with several baselines are performed.", "summary_of_the_review": "The paper combines various recent methods (ROM, PINN, INR) to tackle an interesting problem, namely the acceleration in the resolution of PDEs. Some minor shortcuts hinder the comprehension of the method. To add clarifications on the first and last step of the implementation would enhance the overall quality.", "strength_and_weaknesses": "On overall, the paper is clear, well-written and easy to follow. The method has many advantages compared to traditional ROM approaches, among which applicability to any discretization (which appears as a great deal of current research in dynamical systems) and reduced memory footprint and consumption.\nHowever, some steps in the method are not mentioned in the main paper, only in the appendix, which hinders the overall comprehension of the implementation, as well as its reproducibility. Also, some statements are misleading (developed hereafter).", "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": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "29FtNH14rI", "reviewer": "Reviewer_yRQN", "summary": "", "questions": "", "limitations": "", "rating": 8, "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": "", "summary_of_the_paper": "This paper proposes a new way to solve PDEs using implicit neural representations (INRs) combined with the discrete time discretization/integration and neural network inversion. By leveraging lower dimensional manifolds or learned basis functions via INRs, the proposed method achieved efficient numerical schemes to temporally evolve solutions on a lower dimensional manifold. Also, the high-quality lower dimensional manifold that can be modulated by problem parameters allows more accurate solutions. The experiments demonstrate that the proposed method is efficient and more accurate than a baseline method ROM.\n", "summary_of_the_review": "Overall, this paper is well-written and studies a new combination to solve PDE more efficiently and accurately compared to other alternatives. By learning basis functions/manifolds with INR, the proposed method achieved impressive speed-ups and memory efficiency on top of higher approximate quality. This is a promising result and will inspire many other researches.", "strength_and_weaknesses": "Strengths:\n* The paper is well-written and explains complex concepts and the pipeline very clearly. Considering the numerous mathematical concepts, this paper is impressively concise. Especially the breakdown of the proposed method through Section 4.1 to 4.3 with Figure 3 makes this paper accessible to a general audience, I think.\n* The combination of Neural PDE Solver / ROM and Implicit Neural Networks is an interesting and promising direction, although implicit neural networks are usually nothing but MLPs. \n* The proposed method achieved significant performance gain as well as speed-up.  \n* The proposed method has the flexibility of INR that can be applied to any scale/resolution. \n\t\t \nWeaknesses:\n* The authors did not compare the training time. The key to high-quality solutions here is the quality of manifolds/basis functions. Also, the learned basis functions could be obtained by modulation using mu or features from the encoder e_{\\theta_e}. It usually requires a non-trivial amount of data and training time. The efficiency of training needs to be compared. Also, the model sizes need to be compared for a fair comparison with the baseline method. \n\t\t\nQuestions:\n1. How was g modulated by mu? In the literature, INR is often modulated by Hypernetworks that output the model parameters for another neural network.\n2. Training of the encoder is key to accurate solutions. How do you generate the training data? Is the method sensitive to the distribution of training data?\n", "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": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}]

Author rebuttal:

Decision: Accept: notable-top-25%