Physical Gradients for Deep Learning
Authors: Philipp Holl
Nils Thuerey
Vladlen Koltun
Venue: ICLR 2022 Submitted (2022)
Keywords: deep learning
simulation
optimization
inverse problems
physics

Abstract: Solving inverse problems, such as parameter estimation and optimal control, is a vital part of science. Many experiments repeatedly collect data and employ machine learning algorithms to quickly infer solutions to the associated inverse problems. We find that state-of-the-art training techniques are not well-suited to many problems that involve physical processes since the magnitude and direction of the gradients can vary strongly. We propose a novel hybrid training approach that combines higher-order optimization methods with machine learning techniques. We replace the gradient of the physical process by a new construct, referred to as the physical gradient. This also allows us to introduce domain knowledge into training by incorporating priors about the solution space into the gradients. We demonstrate the capabilities of our method on a variety of canonical physical systems, showing that physical gradients yield significant improvements on a wide range of optimization and learning problems.

Under review as a conference paper at ICLR 2022
PHYSICAL GRADIENTS FOR DEEP LEARNING
Anonymous authors
Paper under double-blind review
ABSTRACT
Solving inverse problems, such as parameter estimation and optimal control, is a vital part of science. Many experiments repeatedly collect data and rely on ma- chine learning algorithms to quickly infer solutions the associated inverse prob- lems. We find that state-of-the-ar training techniques are not well-suited to many problems that involve physical processes since the magnitude and direction of the gradients can vary strongly. We propose novel hybrid training approach that combines higher-order optimization methods with machine learning techniques. We take updates of higher-order domain-specific optimizer and embed them into the gradient-descent- based learning pipeline as physical gradient, replac- ing the regular gradient of the physical process. This also allows us introduce domain knowledge into training by incorporating priors about the solution space into the gradients. We demonstrate the capabilities of our method of canonical physical systems, showing that physical gradients yield significant improvements on a wide range of optimization and learning problems.
variety
INTRODUCTION
Inverse problems that involve physical systems play a central role in computational science. This class of problems includes parameter estimation (Tarantola] 2005) and optimal control (Zhou al., 1996). Solving inverse problems is integral in determining the age of the universe (Freedman et al.,
2001 measuring the amount of dark matter (Padmanabhan & White 2009), detecting gravitational
waves (George & Huertal 20188). controlling plasma flows (Maingi et al. 2019), searching for
neutrinoless double-beta decay (Agostini et al.) 2013) (Aalseth et al. 2018), reconstructing particle trajectories and energies (Belayneh et al. (020). and testing general relativity (Dyson et al. 1920.
Kraniotis Whitehouse 2003).
Decades of research in optimization have produced a wide range of iterative methods and algorithms for solving inverse problems (Press et al.) (2007). Higher-order methods such as limited- umited--emory memory BFGS (Liu& & Nocedal 1989) have been especially successful. These methods have two fundamental drawbacks. First, they require a good initial guess of the solution. When the initial guess far away from the true solution in parameter space, iterative methods may diverge converge to local optimum instead. Second, their computational cost is strongly problem-dependent. as the gradient and possibly even Hessian matrix need to be evaluated for each iteration of the optimization. This may be negligible for simple tasks but when the optimization involves a physical simulation, each iteration can take seconds longer. However, many experiments continuously repeatedly collect data, yielding data sets with millions or billions of individual events, each posing an inverse problem. Solving all of these with iterative solvers poses major challenge which is why, in such settings, data analysis has often turned to machine learning algorithms to quickly infer approximate solutions given the observations (Carleo et al. 2019 Delaquis et al., 2018 George & Huerta 2018b Agostini et al. 2013). The inferred solutions can then, for example, be used to filter data sets for interesting events or can serve as the starting point for a more accurate analysis (CMS Collaboration 2020
George & Huertal 20.
Interestingly, using machine learning models to solve the inverse problems also reduces the chances of divergence or convergence to undesired local optima. Unlike iterative solvers, learning models are jointly trained on large collections of data. The shared parameters allow updates from one example to avoid local minima from different example. Fig. shows generic case illustrating this effect: The task is to localize a wave packet A -in ffxx--exp {(x (x-zz)) //²) and determine its amplitude from noisy recorded time series (grey curve). Iterative solvers fail reconstruct the wave packet

---

Under review as a conference paper at ICLR 2022
when the initial guess is not already close to the true parameters, even diverging in 0.12% of our cases. In contrast, a neural network trained on the same data using the same objective quickly learns to fit all wave packets. Details given in appendix B.I In this paper. we derive hybrid optimization algorithm that integrates higher-order inverse solvers into the traditional deep learning pipeline without imposing these limitations and without discarding the progress made in first-order optimization. Our method can work with arbitrary second-order solvers when the solution space of the inverse problem allows for computation of the Hessian. Otherwise, it requires additional domain knowledge which is used to formulate the inverse solver in more efficient manner. Our method computes the weight updates using both backpropagation as well as the embedded inverse- problem solver. Specifically, the first-order gradient of the physical process used in backpropagation replaced by a physical grudient (PG) derived from the update of the anverse-problem solver. When used in conjunction with training deep learning models, physical gradients allow each update step to encode nonlinear physics, which yiclds substantially more accurate optimization directions. A key advantage of our approach is its compatibility with acceleration schemes (Duchi et al. 201T
Kingma & Ba 2015) and stabilization techniques (Toffe & Szegedy
2015: Huang et al. 2016 Ba
al. 2016) developed for training deep learning models.
The presented hybrid training approach differs fundamentally from pure first-order training. While PGs take the place of regular gradients in the optimization algorithm, they are based on real physical states instead of the adjoint ones calculated in backpropagation. This feature makes it possible to integrate solution priors into the training without restricting the model architecture (Kim et al. 2019), adding soft constraints to the objective (Raissi et al. 2018), or relying exclusively on first- order updates (Rackauckas et al. 2020).
Network L-BFGS-B
0
100
200
We test our method on a wide variety of inverse problems in- cluding the highly challenging Navier-Stokes equations. In all of our experiments, neural networks trained with PGs signif- icantly outperform state-of-the-arr training schemes, yielding accuracy gains of up to three orders magnitude.
Figure 1: BFGS fails to fit wave packet when the initial guess is too far from the true value. A neural network trained on the same objec- tive learns to fit all examples.
METHOD
We consider unconstrained inverse problems that involve differentiable physical process P X C Rd= + Y C Rd, which can be simulated. Here X denotes the physical parameter space and Y the space of possible observed outcomes. Given an observed or desired output y* E Y, the inverse problem consists of finding optimal parameters
1 x² arg min L(x) with L(x) -(x)= ((0)=11P() 1P(xx) (P(x)))) 1"
(1)
z
We are interested in approximating this problem using a neural network, x' = NN(y" [0), parame- terized by 0. Let y = (vi 1, ..., N} denote a set of N inverse problems, all sharing the same physical process P. Then training the network is equivalent to finding
N 1 I 0)) vl
arg min
Σ
This training scheme is unsupervised, i.c. no labels are required over the training data V.
2.1 ITERATIVE OPTIMIZATION
Inverse problems in Eq. I are classically solved by starting with an initial guess zo and iteratively applying updates Art. This yields an estimate of the inverse of the physical model P, ie., x" as P.--(p (y" (2)) where n² denotes the number of iterations required to reach chosen convergence threshold, A key decision in such an iterative optimization scheme is how to derive the update Art. Figure 2a visually compares gradient descent, Newton's method and PGs which will introduce
in section 2.2
2

---

Under review as a conference paper at ICLR 2022
(a)
(b)
Update
S HF NL
3
PG
dx -1 ( -y)
6
GD
GD
v
08
N
- (P(s) y y') ('(00)
2
0xx 00
0
-1
(x - Pv (v)))²
X1
y1
Figure 2: (a) Minimization of the function L 11 with y (z1, 22) in space. Contours of are shown in gray. Solid lines are the optimization trajectories gradient descent (GD), Newton's method (N), and physical gradients (PG), with infinitesimal step sizes. Circles represent the first 10 iterations with constant step size. (b) Comparison of the same optimization methods by their respective update steps and properties: whether they adapt to function sensitivity (S), are Hessian-free (HF) and take nonlinearities into account (NL.). H denotes the Hessian w.r.t. The properties of PG explained in section 2.2
First-order methods The scale of typical machine learning datasets necessitates the use of mini- batches, which has led popular optimization methods in machine learning to rely only first- order derivatives
These Hessian-free methods are popular deep learning due their low computational cost and stable convergence (Goodfellow et al. 2016). The simplest such method is
gradient descent (GD) (Curry 1944) and its stochastic version (SGD) where Az = -1- (OL) Note that 0L and therefore of Az, scales inversely with æ. When carries physical units, neither Ar 8x nor Arr lie in the same space z. This scaling behavior is problem when dealing with sensitive or insensitive functions, i.c functions whose gradient magnitudes de changes in cause large changes in L, this counter-intuitiv intuitively leads to large updates Az when using GD, resulting in even larger changes to L. This is often referred to as exploding gradients, and the opposite effect, vanishing gradients, occurs with insensitive functions. Machine learning models can be tuned to behave well given these updates (Toffe & Szegedy 2015 Loshchilov & Hutter] (019) but in physics-based optimization, where functions may be extremely sensitive in certain parts of the parameter space or on subsets of training examples, GD often prescribes suboptimal optimization directions.
Bz
2L are far from 1. When small
While there are a number of first-order optimizers that approximate higher-order information to improve convergence (Kingma & Ba, 2015 Hestenes et al. 1952] Duchi al.) 2011 Martens & Grosse 2015), they cannot take nonlinearities into account, which results in suboptimal optimization directions when optimizing nonlinear functions.
Second-order methods Newton's 's method (Atkinson 2008) uses the inverse Hessian determine
the optimization update Ax = -
8 81 0z
T
Direct computation of the Hessian adds
02
significant computational cost each iteration. Quasi-Newton methods (Broyden) 1970 Liu &
Nocedal 1989 Gill & Murray 1978 More 1978 Powell 1970 Berndt et al. 1974 Conn et al.
1991 Avriel 2003) alleviate this issue by approximating the inverse Hessian instead. Unlike GD, Newton-type methods prescribe small updates to sensitive functions and large updates to insensi-
tive ones. The resulting update steps typically progress towards optima much faster than first-order methods (Ye et al., 2019). However, the difficulty of approximating the Hessian for batched data sets (Schraudolph et al. (007) and the high cost of evaluating the Hessian directly have largely pre- vented second-order methods from being utilized in machine learning research (Goodfellow et al. 2016).
Domain-specific inverse solvers Domain knowledge can often used to formulate more effi- cient optimization updates directly approximate solution. A classical example would be preconditioners that, when chosen correctly, allow linear solvers to converge much more quickly and reliably. In nonlinear settings, part of the governing equations may be inverted analytically or
3

---

Under review as a conference paper at ICLR 2022
priors about the solution space may be incorporated into the solver, for example. A time-reversed simulation may be able to more accurately recover prior states than a generie inverse solver. These methods are commonly used in optimization but their use has cluded machine learning applications so far due to the fixed gradient-descent pipeline. When learning with PGs, domain-specific optimiz- can he integrated into training like any other optimizers, allowing for greatly training convergence. Appendix appendix A.2 lists necessary conditions that need to be fulfilled.
2.2 PHYSICAL GRADIENTS FOR DEEP LEARNING
We now consider the problem of learning solutions to the inverse problems as in Eq. 2 As discussed above, GD-based optimizers fail to take the nonlinearity and potentially strongly varying sensitivity of the joint problem into account while higher-order optimizers are hard to apply directly.
Instead, we aim to leverage the advantages of higher-order optimizers in a deep learning setting by embedding arbitrary inverse-pree prensse-roblemmmmme solvers physical gradients (PGs) into the first-order opti- mization pipeline. The network updates A can then be computed using any first-order optimizer such SGD Adam. Additionally, all state-of-the-art deep learning techniques, such as normal- ization loffe & Szegedy (2015); [a et al." (2016) or dropout (Srivastava et (2014)) can be employed, which often improves training speed and generalization performance Keskar et al. (2017).
Our method can be applied whenever one of the following available: (i)
8pp can be computed Ox""
can be derived analytically, or (iii) an inverse-physics solver can be derived using the available domain knowledge. If (i) is fulfilled, we can explicitly compute the Newton direction in in x space. Thus, with any of (i), (ii) or (iii) we have
numerically, (ii) a higher-order update such as 0xs
82p
aP
de
JP These Az are the basis for our method which we access to better updates Ar than ATeD x Oe describe in the following.
We derive our method starting from Eq. 2 which describes the combined problem of finding both solutions xn to the individual inverse problems y" E y and network weights that approximate them. One way to solve this problem is to first consider all inverse problems separately and precom- pute corresponding solutions x = {-,²² Pax)") = y } using an iterative optimizer. Then, x could be used as labels for supervised training of neural network:
N Σ
0. arg min
NN () - 22" -0-=>'1133
(3)
This precomputation makes possible use an efficient optimization method for x while training the neural network with a first-order optimizer, retaining the ability to employ state-of-the-art deep learning techniques, such as normalization (Toffe & Szegedy 2015 Ba et al., 2016) or dropout (Sri-] vastava et al.) 2014). However, this approach has severe drawbacks. For one, the individual opti- mizations can more easily get stuck in local optima than the combined training (Holl et al., 2020). Also, many inverse problems are inherently multi-modal, i.c. multiple solutions x' exist for one y" In such settings z heavily depends on the initial guess 2p used for precomputation, which may cause the network to interpolate between possible solutions, leading subpar convergence and generalization performance.
To avoid these problems, alter the training procedure from Eq. 3 in two ways. First, treat the inverse problem solver P-1 as part of a coupled optimization and consider it as part of the training loop, yielding
0, arg min
|NN (0) PM: (2:)
Next we can condition the embedded inverse problem solver on the neural network prediction, P.,.! (ui) - PT. (vf NN( [0)) by using it as an initial guess. The embedded optimizer can now find minimum close the prediction, which allows the neural network to choose any solu- tion of a multi-modal problem. Also, since all inverse problems from y" are optimized jointly, this reduces the likelihood that an individual solution gets stuck in local minimum.
is furthermore not necessary to run the iterative inverse solver to convergence. Any choice of N steps likely leads to convergence (proof given in appendix A.2), and choosing small n can
4

---

Under review as a conference paper at ICLR 2022
increase the training speed. We also observe this in practice: Optimizations with n converged reliably in all our experiments. With these modifications, we arrive at
0. = arg mii
- [NN (0) - P (2)' IN 10)
(5)
The central quantity here the difference between prediction and correction obtained from Pr= which we refer to the physical gradient. It encodes valid physical state and takes the place of (PP) T. which would otherwise be computed by backpropagation.
(ox (OP),,
Since the PG stems from a generic higher-order or domain-speeinee solver, the resulting updates A will also be non-linear without computing the Hessian w.r.t. 0. Assuming the physics optimizer can account for varying function sensitivity its updates Az, the resulting PG updates A will inherit that property since the network is generally tuned to retain gradient magnitude. The properties of PG as well as GD and Newton's 's method are summarized Fig. D
While the L22 measure in Eq. looks like supervised objective in space, the actual network optimization does not behave as such. To see this more clearly, consider p-1 performing a single
GD step with unit learning rate, (y- .) = ( (z)x) y²). Then the gradient for Eq
aP 0z is de 20
(P(z) y²), which is identical to the unsupervised training in Eq. (the proof is given
in appendix
Algorithm 1: Neural network training with physical gradients according to Eq.
for each training sample do Fo NN (0)
Solution inference with NN
for k = 1,
do
3& + P-1 (y² [K
Accumulate iterative updates for PG
end
AO
00
Gr'
(xo In)²
Backpropagation through Parameter update
0 - Step 0 A0,m)
end
Algorithm I details the corresponding neural network training procedure. There, the effective gra- dient &0 can also be written as 39 2g z,.12 with v. taken as constant as in Eq. which lends itself to straightforward integration into machine learning frameworks where a simple T loss can be used to optimize the network. This L2 loss acts as a proxy to embed the physical gradient into the network training pipeline and is not to be confused with traditional supervised loss in space (see appendix A.I for a more detailed discussion).
RESULTS
We first demonstrate the differences between various neural network training schemes on sim- ple one-parameter example before moving on to physical systems described by partial differen- tial equations. We specifically consider Poisson's equation, the heat equation, and the Navier- Stokes equations. This selection covers ubiquitous physical processes with diffusion. transport, and strongly non-local effects, featuring both explicit and implicit solvers. each experiment, we train identically-initiaiizee tive differences in convergence behavior. For the weight update, we primarily employ the Adam optimizer [Kingma & Ba 2015), standard method in deep learning that outperforms SGD in our experiments. All training data are randomly generated on the fly, resulting in effectively infinite data sets y All shown learning curves therefore representative of the performance unseen data. A detailed description of our experiments along with additional visualizations and performance measurements can be found in appendix B
networks using different gradient schemes for P to examine the rela-
5

---

Under review as a conference paper at ICLR 2022
10
Mean sin
10 -3
Adam supervised
10 5
Adam
Adam Newton
Adam Inv.phys
10
5000
10000
15000
20000
lteration
Figure 3: Network trained single-parameter optimization (section 3.1 using Adam in combina- tion with various gradient schemes. Running average over 1000 mini-batches.
3.1 SINGLE-PARAMETER OPTIMIZATION
First, we consider the task of finding solutions x² E to the inverse problem (Eq. Pa (ine)=si(((])))
1, given a E [0.1, 10.1). This problem may seem simple first glance but P shares some properties with chaotic physical systems, which makes this problem hard to optimize. We train neural network with three fully-comnected hidden layers to solve this task; the learning curves are shown in Fig. 3
Supervised training schemes (Eq. 3) perform poorly this task due multi-modality. Using Adam for the joint optimization with backpropagation through the network and physics (Eq.2) avoids the aforementioned problem as the gradients dynamically guide the optimization towards proximal minimum. However, 8P oscillates without bound when x - 0, which causes overflow errors in the dr. optimization. To avoid this, clip the gradients @P 0r to [-1, which keeps the optimization relatively stable, although in most cases it never fully converges. This demonstrates one of the key weaknesses of first-order learning: when the magnitude of the gradients varies strongly between examples, the optimization is either slow or unstable, depending on the learning rate.
Finally, we test two variants of physical gradients (Eq. 5)) this problem: Newton's method and an analytic solver. As Fig. shows, the PG methods vastly outperform both supervised and unsu- pervised first-order training, converging
accuracy of around 10-5 For this example, using PG with the inverted physics results in a 1000x improvement in accuracy over conventional network training, while being more stable. The gradients of all variants are visualized in appendix B.2
3.2 POISSON'S EQUATION
Poisson's 's equation, ( 2r, plays an important role in electrostatics, Newtonian gravity, and fluid dynamics (Ames 2014). It has the property that local changes in x can affect ((x) globally. Here we consider a -vo-dimnnsional system and train U-net (Ronneberger al. 2015) to solve inverse problems (Eq. D on pseudo-randomly generated y". Fig. 4p shows the learning curves.
We construct PGs based on the analytic inverse and couple them with Adam for training the neural network. When learning with Adam or SGD with momentum, learning drastically slows after around 200 and 300 iterations, respectively. Meanwhile, the learning curve involving PGs closely resembles an exponential curve, which indicates linear convergence, the ideal case for first-order methods optimizing an Le objective. During all of training, the PG variant converges exponentially faster than the first-order alternatives, its relative performance difference compared to Adam continually increasing from a factor of 3 at iteration 60 to full order of magnitude after 5k iterations.
3.3 HEAT EQUATION
Next, we consider system with fundamentally non-invertible dynamics. The heat equation, V/ 22u, models heat flow in solids but also plays part in many diffussve systems (Droniou
8u. 7/
2014
It gradually destroys information the temperature equilibrium is approached (Grayson 1987
causing VP to become sing-singular. Inspired by heat conduction in microprocessors, we generate examples øer by randomly scattering four to ten heat generating rectangular regions plane and

---

Under review as a conference paper at ICLR 2022
(a)
y
YAdam
+
XAdam
XA
x
(b)
(c)
10
Inferred MAE
Oplim.X MAE
10"
n
10
10
GD, 6)
SGD
Adam, 10-3
Adam PG phw.ph
10
L-BFGS-B
10²
phys.
-3
40
60
100
0
2500
5000
7500
10000
Optimization Iteration
Training Iteration
Figure 4: Optimization of Poisson's equation to reconstruct a desired output state (section 3.2 (a) Example from the data observed distribution (y²) and inferred solutions, ground truth solu- tion (z*), (b) Convergence curves of various optimization methods on a random example. (c) Neural network learning curves, running average over 64 mini-batches. The vertical axes are logarithmic.
(a)
XAdam
XA
XA
(b)
(c)
500
GD n= 1
Adam
Optim, X MAE
L-BFGS-B
Inferred x MAE
500
Adam PG (L-BFGS-B)
Inv. phys.
Adam PG (Inv.phys.)
200
100
100
20
40
60
80
100
2500
5000
7500
10000
Optimization Iteration
Training lteration
Figure 5: Optimization involving the heat equation (section 3.3 (a)
Example from the data set:
observed distribution (y"), inferred solutions, ground truth solution (x:"). (b) Optimization curves for one example. (c) Learning curves.
simulating the heat profile y" P(IGT) observed from outside a heat-conducting casing. We train a U-net (Ronneberger et al. 2015) to solve the corresponding inverse problem (Eq. I.. The
learning curves are shown in Fig. 5.
Since aP is stable, individual inverse problems can be optimized with GD (see Fig. 5p) For the dz unsupervised network training (Eq. 2). use Adam with 7/ 10-3 and observe that the distance to the solution starts rapidly decreasing before decelerating between iterations 30 and 40 to a slow but mostly stable convergence. The sudden deceleration is rooted in the adjoint problem, which is also a diffusion problem. Backpropagation through P removes detail from the gradients, which makes hard for first-order methods to recover the solution.

---

Under review as a conference paper at ICLR 2022
For training with PGs (Eq. 5), we consider L-BFGS-B B and the analytically derived inverse-physien solver L-BFGS-B B outperforms GD when optimizing single inverse problems but is limited by the ill-conditioning (see Fig. (5) 5b) of the problem. Instead, we use the available domain knowl- edge to formulate an inverse problem solver P-1 that finds approximate x= in a stable manner by introducing probabilistic reasoning into the algorithm (see appendix (B.4))
Using L-BFGS-B with n = 32 as a PG to train the network leads to solutions that are about 20% more accurate and noticeably sharper after the same number of training iterations. This advantage is achieved early during training and stays relatively constant throughout. From iteration 40 on, L-BFGS-B training converges about as fast as pure first-order training. However, the 17x longer computation time per iteration largely negates these advantages.
The PGs based on the inverse-physie solver, on the other hand, manage to improve the convergence speed substantially. They perform on par with L-BFGG-B B training for the first ~35 iterations. How- ever, afterwards the increase in accuracy does not slow down as much as with the other methods, resulting in an exponentially faster convergence. At iteration 100, the predictions are around 34% more accurate compared to Adam, and the difference increases to 130% after 10k iterations. To reach the accuracy of regular Adam training with 10k iterations, L-BFGS-B training requires 628 iterations and training with the inverse- phersics solver requires only 138 iterations.
Additional examples and learning curves are shown in section appendix B.4
3.4 NAVIER-STOKES EQUATIONS
Fluids and turbulence are among the most challenging and least understood areas of physics due to their highly nonlinear behavior and chaotic nature (Galdi, 2011). We consider a two-dimensional system governed by the incompressible Navier-Stokes equations: OH 2/2 V Ve Vp, v-v= 0, V x p = O, where p denotes pressure and v the viscosity. At = 0, a region of the fluid is randomly marked with a massless colorant mg that passively moves with the fluid, am - Vm, कं and after a time f. the marker is observed again, my. An example observation pair y² = {mo, m,) is shown in Fig. 6a. We target the inverse problem (Eq. ( of finding an initial fluid velocity x = vo such that the fluid simulation P matches my at time t. Since P is deterministic, x encodes the complete fluid flow from 0 to We define the objective in frequency space with lower frequencies being weighted more strongly. This definition considers the marker distribution match on all scales, from the coarse global match to fine details, and is compatible with the definition in Eq. [ We train a U-net (Ronneberger et al., 2015) to solve these inverse problems; the learning curves are shown in Fig. Cc.
de
For unsupervised training (Eq. 2 with Adam, the error decreases for the first 100 iterations while the network learns to infer velocities that lead to an approximate match The error then proceeds to decline at a much lower rate, nearly coming to standstill. This is caused by an overshoot in terms of vorticity, as visible Fig. right. While the resulting dynamics can roughly approximate the shape of the observed m,, they fail to match its detailed structure. Moving from this local optimum to the global optimum is very hard for the network as the distance in x space is large and the gradients become very noisy due the highly non-linear physics. A similar behavior can also be seen when optimizing single inverse problems with GD (F)..66)). Here the strongly varying magnitude of the gradients poses an additional challenge for the first-order optimization, causing its progress to stall. With 1, it takes more than 20K iterations for GD to converge for single problems.
For training with PGs (Bq. 5). we exploit our knowledge about the flow e.g. smoothness to formulate an inverse simulator P-1 that incorporates corresponding priors on the solution space. It runs the simulation in reverse, starting with my, then estimates z by comparing the marker densities from the forward and reverse pass (see appendix B.5 for details). This method converges for single inverse problems within around 100 iterations and its computational costs are on the same level as evaluating the gradient. When used for network training, we observe vastly improved convergence behavior compared to conventional first-order training. The error rapidly decreases during the first 200 iterations, at which point the inferred solutions are more accurate than pure Adam training by a factor of 2.3. The error then continues to improve at a slower but still exponentially faster rate, reaching a relative accuracy advantage of 5x after 20K iterations. To match the network trained with pure Adam for 20K iterations, the PG variant only requires 55 iterations. This improvement is possi- ble because the inverse-physics solver and associated PG does not suffer from the strongly-varying
8

---

Under review as a conference paper at ICLR 2022
(a)
mo
Yddam
YA+ PG
Xadam
XA+ PG
(b)
500
(c)
GD
Optim. y MSE
Inv.phys.
Inferred x MAE
10
100
5
Adam
Adam PG (Inv.phys)
2
20
10
10
10
10
10"
0
5000
10000
15000
20000
Optimization Iteration
Training Iteration
(a) One example from the data set: initial marker distribution (mo) simulated marker distribution after time t using ground-truth velocity (y") and network predictions (y.): predicted initial velocities (x.); ground truth velocity (x*). (b) Optimiza-
Figure 6: Incompressible fluid flow (section 3.4
tion curves with = 1, averaged over 4 examples. (c) Learning curves, running average over 64 mini-batches.
Table 1: Time to reach equal solution quality in the fluid experiment, measured MAE in
Method
Training time
Inference time per example
Neural network with PG
17.6 (15.6 iterations) n/a n/a
0.11 1 ms (immediate)
nomain-speciffi solver
2.2 (7 iterations)
Gradient descent optimizer
> 4h (20k iterations)
gradient strengths and directions, which drown the first-order signal in noise. Instead, physical gradients behave much more smoothly, both in magnitude and direction.
To reach the same solution quality the neural network prediction, the domain knowledge solver needs more than 10,000 times as long. This difference is caused by solver having to run the full forward and backward simulation several times. Table measured training and inference times. With both iterative solver and network, we used batch size of 64 and divided the total time by the batch size, More examples from the fluid data set including time evolution sequences are shown in
appendix B.5
CONCLUSION
We have shown how to integrate arbitrary optimization schemes for inverse problems into machine learning settings. Our results show that first-order optimization is suboptimal in many situations, especially when the magnitude or direction of gradients can vary during the optimization or from example example. While more advanced optimizers such as Adam alleviate this behavior, they fail to determine good optimization directions for many problems. By embedding higher-order optimizers and especially inverse-physics solvers into the first-order training pipeline, we were able to mitigate these shortcomings without discarding the progress made in the field of machine learning. The resulting training scheme converges exponentially faster than traditional first-order training wide variety of inverse problems involving partial differential equations while being on par with regular gradient evaluation terms of computational cost.
We anticipate method to be applicable to wide range of differentiable processes. Beyond physics, interesting avenues lie in applying PGs to differentiable rendering or training invertible neural networks.

---

Under review as a conference paper at ICLR 2022
REPRODUCIBILITY STATEMENT
We will publish the full source code and trained network models upon acceptance. They will be sufficient retrain the networks and plot the figures shown in this work. All data sets were generated randomly, so slight variations are expected. Repeated experiments are shown in the appendix.
REFERENCES
Craig Aalseth, N Abgrall, Estanislao Aguayo, SI Alvis, M Amman, Isaac J Arnquist, FT Avi- gnone III, Henning O Back, Alexander S Barabash, PS Barbeau, et al. Search for neutrinoless double-3 decay in ge 76 with the majorana demonstrator. Physical review letters, 120(13):132502, 2018.
Martín Abadi, Paul Barham, Jianmin Chen, Zhifeng Chen, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Geoffrey Irving, Michael Isard, et al. Tensorflow: system for large- scale machine learning. In 12th {USENIX) symposium on operating systems design and imple- mentation ({OSDI) 16), pp. 265-283, 2016.
M Agostini, M Allardt, E Andreotti, AM Bakalyarov, M Balata, Barabanov, M Barnabé Heider, N Barros, L Baudis, C Bauer, et al. Pulse shape discrimination for gerda phase data. The European Physical Journal C, 73(1))2558, 2013.
William Ames. Numerical methods for partial differential equations. Academic press, 2014.
Kendall E Atkinson. An introduction to numerical analysis. John wiley & sons, 2008
Mordecai Avriel. Nonlinear programming: analysis and methods. Courier Corporation, 2003.
Jimmy Lei Ba, Jamie Ryan Kiros, and Geoffrey Hinton. Layer normalization. In International Conference on Learning Representations (ICLR), 2016.
Dawit Belayneh, Federico Carminati, Amir Farbin, Benjamin Hooberman, Gulrukh Khattak, Miaoyuan Liu, Junze Liu, Dominick Olivito, Vitória Barin Pacela, Maurizio Pierini, et al. Calorimetry with deep learning: particle simulation and reconstruction for collider physics. The European Physical Journal C, 80(7):1031. 2020.
Emst R Berndt, Bronwyn H Hall, Robert E Hall, and Jerry A Hausman. Estimation and inference in nonlinear structural models. In Annals of Economic and Social Measurement, Volume 3, number 4, pp. 653-665, NBER, 1974.
Charles George Broyden. The convergence of class of double-rank minimization algorithms 1. general considerations. IMA Journal of Applied Mathematics, 6(1):76-90) 1970.
Giuseppe Carleo, Ignacio Cirac, Kyle Cranmer, Laurent Daudet, Maria Schuld, Naftali Tishby, Leslie Vogt-Maranto, and Lenka Zdeborová. Machine learning and the physical sciences. Re- views of Modern Physics, 91(4):045002, 2019.
CMS Collaboration. A deep neural network to search for new long-lived particles decaying to jets. Machine Learning: Science and Technology, 13))0035012, 2020.
Andrew R Conn, Nicholas IM Gould, and Ph Toint. Convergence of quasi-newton matrices generated by the symmetric rank update. Mathematical programming, 50(1-3):177-15.. 1991.
John J Craig. Introduction robotics: mechanics and control. Prentice Hall, 2005.
Haskell B Curry. The method of steepest descent for non-linear minimization problems. Quarterly of Applied Mathematics, 2(3):258-28 2211, 1944.
George Cybenko. Approximation by superpositions of sigmoidal function. Mathematics of control, signals and systems, 2(4):303-314, 1989.
S Delaquis, MJ Jewell, I Ostrovskiy, M Weber, T Ziegler, Dalmasson, LJ Kaufman, T Richards, JB Albert, G Anton, et al. Deep neural networks for energy and position reconstruction exo- 200. Journal Instrumentation, 13(08): 2018.
10

---

Under review as a conference paper at ICLR 2022
Jerome Droniou. Finite volume schemes for diffusion equations: introduction and review of modern methods. Mathematical Models and Methods in Applied Sciences, 24(08):1555-1619 2014.
Simon Du, Xiyu Zhai, Barnabas Poczos, and Aarti Singh. Gradient descent provably optimizes over-parameterized neural networks. arXiv preprint arXiv:1810.02054. 2018.
John Duchi, Elad Hazan, and Yoram Singer. Adaptive subgradient methods for online learning and stochastic optimization. Journal of machine learning research, 12(7), 2011.
F. W. Dyson, A. S. Eddington, and C. Davidson. A Determination of the Deflection of Light by the Sun's Gravitational Field, from Observations Made at the Total Eclipse of May 29, 1919, January 1920. htttss:::: //doi wtwp:/o.orr..... 098/ rsta. 1900. 0009
Wendy Freedman, Barry F Madore, Brad K Gibson, Laura Ferrarese, Daniel D Kelson, Shoko Sakai, Jeremy R Mould, Robert C Kennicutt Jr, Holland C Ford, John A Graham, et al. Final results from the hubble space telescope key project to measure the hubble constant. The Astro- physical Journal, 553(1):47, 2001.
Giovanni Galdi. An introduction to the mathematical theory of the Navier-Stokes equations: Steady- state problems. Springer Science & Business Media, 2011.
Daniel George and EA Huerta. Deep neural networks to enable real-time multimessenger astro- physics. Physical Review D, 97(4):044039, 2018a.
Daniel George and EA Huerta. Deep learning for real-time gravitational wave detection and param- eter estimation: Results with advanced ligo data. Physics Letters B, 71777664770 70, 2018b.
Philip E Gill and Walter Murray. Algorithms for the solution of the nonlinear least-squares problem.
SIAM Journal on Numerical Analysis, 15(5):977-99 1978.
Alexander Gluhovsky and Christopher Tong. The structure of energy conserving low-order models. Physics of Fluids, 1(())334-343, 1999.
Ian Goodfellow, Yoshua Bengio, Aaron Courville, and Yoshua Bengio. Deep learning. MIT press Cambridge, 2016.
Matthew A Grayson. The heat equation shrinks embedded plane curves to round points. Journal of Differential geometry, (6(2):285-314 1987.
FH Harlow. The marker-and-cell method. Fluid Dyn. Numerical Methods, 38, 1972.
Francis H Harlow and J Eddie Welch. Numerical calculation of time-dependent viscous incompress- ible flow of fluid with free surface. The physies fluids, 8(12)2182-2118,, 1965.
Magnus Hestenes, Eduard Stiefel, et al. Methods of conjugate gradients for solving linear systems. Journal of research of the National Bureau of Standards, 49(6):409 436, 1952.
Philipp Holl, Vladlen Koltun, and Nils Thuerey. Learning to control pdes with differentiable physics. In International Conference on Learning Representations (ICLR), 2020.
Gao Huang, Yu Sun, Zhuang Liu, Daniel Sedra, and Kilian Q Weinberger. Deep networks with stochastic depth. In European conference on computer vision, pp.
Springer, 2016.
Sergey Ioffe and Christian Szegedy. Batch normalization: Accelerating deep network training by reducing internal covariate shift. In arXiv:/502.03 167, February 2015.
Nitish Shirish Keskar, Dheevatsa Mudigere, Jorge Nocedal, Mikhail Smelyanskiy, and Ping Tak Pe- ter Tang. On large-batch training for deep learning: Generalization gap and sharp minima. In International Conference Learning Representations (ICLR), 2017.
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(2). pp- 59-70. Wiley Online Library, 2019.
Diederik Kingma and Jimmy Ba. Adam: A method for stochastic optimization. In International Conference on Learning Representations (ICLR), 2015.

---

Under review as a conference paper at ICLR 2022
GV Kraniotis and SB Whitehouse. Compact caleulation of the perihelion precession of mercury in general relativity, the cosmological constant and jacobi's inversion problem. Classical and Quantum Gravity, 20(22):4817, 2003.
Dong Liu and Jorge Nocedal. On the limited memory bfgs method for large scale optimization. Mathematical programming. 45(1-3):503-528, 1989.
Ilya Loshchilov and Frank Hutter. Decoupled weight decay regularization. International Confer- ence on Learning Representations (ICLR), 2019.
Rajesh Maingi, Amold Lumsdaine, Jean Paul Allain, Luis Chacon, SA Gourlay, CM Greenfield, JW Hughes, Humphreys, V Izzo, H McLean, et al. Summary of the fesac transformative enabling capabilities panel report. Fusion Science and Technology. 75(3):167-177, 2019.
James Martens and Roger Grosse. Optimizing neural networks with kronecker- approximate curvature. In International conference machine learning, pp. 2408-2411 PMLR, 2015.
Jorge J Moré. The levenberg-marquarde algorithm: implementation and theory. In Numerical anal- ysis. pp. 105-116. Springer, 1978.
Patrick Mullen, Keenan Crane, Dmitry Pavlov, Yiying Tong, and Mathieu Desbrun. Energy- preserving integrators for fluid animation. ACM Transactions on Graphics (TOG), 28(3):1-8, 2009.
Nikhil Padmanabhan and Martin White. Calibrating the baryon oscillation ruler for matter and halos. Physical Review D, 8066)006308,, 2009.
Adam Paszke, Sam Gross, Soumith Chintala, Gregory Chanan, Edward Yang, Zachary DeVito, Zeming Lin, Alban Desmaison, Luca Antiga, and Adam Lerer. Automatic differentiation in pytorch. In Advances in Neural Information Processing Systems, pp. 1-8, 2017.
Michael Powell. A hybrid method for nonlinear equations. In Numerical methods for nonlinear algebraic equations. Gordon and Breach, 1970.
William H. Press, Saul A. Teukolsky, William T. Vetterling, and Brian P. Flannery. Numerical Recipes. Cambridge University Press, 3 edition, 2007. ISBN 9780521880688.
Christopher Rackauckas, Yingbo Ma, Julius Martensen, Collin Warner, Kirill Zubov, Rohit Supekar, Dominic Skinner, and Ali Ramadhan. Universal differential equations for scientific machine learning. arXiv:2001.04385, 2020.
Maziar Raissi, Alireza Yazdani, and George Em Karniadakis. Hidden fluid mechanics: A navier-stokes informed deep learning framework for assimilating flow visualization data. arXXv:1808.04327 2018.
2018.
Olaf Ronneberger, Philipp Fischer, and Thomas Brox. U-net: Convolutional networks for biomedi- cal image segmentation. In International Conference on Medical image computing and computer- assisted intervention, pp. 234-241. Springer, 2015.
Nicol N Schraudolph, Jin Yu, and Simon Günter. A stochastic quasi-newton method for online convex optimization. In Artificial intelligence and staristics, pp. 436 436-443, 2007.
Andrew Selle, Ronald Fedkiw, ByungMoon Kim, Yingjic Liu, and Jarek Rossignac. An uncon- ditionally stable maccormack method. Journal of Scientific Computing, 35(2-3):350-371, June doi: 10. (0.10077/1001155077-9166
2008. ISSN 0885-7474, 1573-7691.
Nitish Srivastava, Geoffrey Hinton, Alex Krizhevsky, Ilya Sutskever, and Ruslan Salakhutdinov. Dropout: a simple way to prevent neural networks from overfitting. The journal of machine learning research, 19(1):1929-199 1958, 2014.
Albert Tarantola. Inverse problem theory and methods for model parameter estimation. SIAM, 2005.
Nan Ye, Farbod Roosta-Khorasani, and Tiangang Cui. Optimization methods for inverse problems. In 2017 MATRIX Annals, pp. 121-140. Springer, 2019.
Kemin Zhou, John Comstock Doyle, Keith Glover, al. Robust and optimal control, volume 40. Prentice hall New Jersey, 1996.
12

---

Under review as a conference paper at ICLR 2022
A
METHOD
A.
EMBEDDING THE PHYSICAL GRADIENT INTO THE GRADIENT DESCENT PIPELINE
The key for embedding tivers-physiee optimizers in the machine learning pipeline is the proxy L2 loss in Eq. [ which passes the gradients from the physics optimizer to the GD-based neural network optimizer.
This la loss may look like a supervised loss in x space but due to the dependency of the labels P- on the prediction, it behaves very differently. To demonstrate this, we consider the case that GD is being used as the physics optimizer. Then the total loss is purely defined in y space, reducing to regular first-order optimization. The proxy L2 loss simply connects the computational graphs used in backpropagation.
Physical gradients based on gradient descent are equal unsupervised training With the shorthands NN(x, (0) yı P(x;) and Ay; y/: vi, we can write the objective functions of unsupervised training and physical gradient training as
N
U(0)
= Σ
(unsupervised)
(6)
N
M(0) = Σ 2 P [=1
2;)1 (physical gradients)
(7)
Theorem 1. Minimizing (6) identical to minimizing (() using any gradient-descent
optimizer P is implemented as single gradient descent step with 17 = 1.
Proof. The gradient of the unsupervised objective function is
aU
Oy Ox
Ati
00
dx 00
For the mixed training, we recall that the objective function of single inverse problem is L(x) }1[p(z) y²1 The corresponding gradient is O/2 Ay-
By Inserting gradient descent step for
De
P-11 the gradient of the mixed objective function becomes
N
Ox 00
=
Σ
(I.
- P
-2 xx))
00
1=1
N
aL
8x 00
2i 17
dx
Oy Ox
n
A y;
Oz 00
which is equal to aU for 1) 1.
This generic La formulation is even applicable in settings where different parts of an end-to-end pipeline are computed by different software frameworks or hardware accelerators.
A.2 CONVERGENCE WITH PHYSICAL GRADIENTS
We consider the convergence of a coupled optimization of L,(x;) with .; = fo(y)) involving gradi- ent descent optimization of general function approximator fe and some other optimization scheme P for minimizing objective functions L3 (Zx4). Here, i denotes the index of a specific example from a data set of size N. In the context of learning to solve inverse problems, P is the inverse- inverse-problem problem solver which physical gradient P(x) æ is based, and fe represents the neural network. In the main text, we considered objectives of the form L; = (x;;) y l|2. but now we consider a more general form.
13

---

Under review as a conference paper at ICLR 2022
B
Xn
01.
0x Km
Figure 7: Convergence visualization in space for one example The point v, So(v) represents the current solution estimate. The grey line and area mark [x] < [(xn)] and x" minimum of L. H, (œv,) is the next point in the optimization trajectory of and the green and blue circles around represent open sets with radii and 1| - In ll. respectively. the area shaded in orange, the distance to f, decreases but L increases.
Assumptions
1. Let (L; RR R I...., N) be a finite set of functions and let x e Rd be any global minima of Li.
2. Let P = =(P Rd + B = 1,.., A,.......... be a set of update functions for Ls such that Vr Rd
Br 0 Vr Rd [((x) L.((P())) 2 r (L(z) L(x")) and 3K
[P((r)) -x ((L(x)) - L(x))).
3. Let So Rm + Rdd be differentiable w.r.t. 0 with the property that By > : Vi E
1. ..., N Va, E Rd ve 0 In E EN I15o3 (y.) xall2 S where Or is the sequence of
gradient descent steps with 0= 0 0% 17
0fu (So.. x).,))n1 0.
80
Here make minimal assumptions about L; not even continuity is required long have access to some form of a physical gradient that can optimize in x. We denote repeated applica- tions of P by pk (x). The assumption in (3) states that the function fo is flexible enough fit our problem. It must be able to converge to every point æ for all examples using gradient descent. has been shown that neural networks with sufficiently many parameters fulfill this condition under certain assumptions (Du et al. 2018) and the universal approximation theorem guarantees that such configuration exists even for shallow neural networks with enough parameters (Cybenko) 1989), We now show that P can be used in combination with gradient descent to optimize the coupled problem.
Theorem 2. For all E N, there exists an update strategy 0x+1 Ux(O,) based on single evaluation of pk for which L(fo, (y:)) converges to a minimum x; or minimum region of L; Vi 1.....N. N.
Proof. Dropping the example index i, we denote 2, ->*(x,) and AL* = L(xx) L(x*). Let ln denote the open of all for which L(x) I(xa) AAT Definition (2) provides that J. Tz-
14

---

Under review as a conference paper at ICLR 2022
: Vx E (..) [(x,) - L(x) > 2 AAL" where B, denotes the open set containing all s' Rd for which || æl12 < <, i.e. there exists small ball around ön which is fully contained in 12 (see sketch in Fig. 7
Since 1n is open, Be 0
Using the convergence assumption for fa. definition (3), we can find a finite N for which fo. E (Tn) and therefore L(x,) - L(fo,, JAL". We can thus use the following strategy Un for minimizing :: First, compute 8, p Then perform gradient descent steps in 0 with the effective objective function 1 fo En l2 until L(x,) L(fe) 2 SAL'. With this strategy, each application of Us reduces the loss to + S (1 A so any arbitrary value of <> L(x*) can
be reached within a finite number of steps. Definition (2) also ensures that [[(()) - =|| - as the
optimization progresses which guarantees that the optimization converges minimum region. O
While this guarantees convergence, it requires potentially many gradient descent steps in 0 for each physical gradient evaluation. This can be advantageous in special circumstances, e.g. when the physical gradient is more expensive to compute than an update in 8, and 0 is far away from solution. However, many cases, we want to update the physical gradient after each update to O. Without additional assumptions about P and Se, there is no guarantee in this case that L will decrease every iteration, even for infinitesimally small step sizes. Despite this, there is good reason to assume that the optimization decreases L over time.
When performing a gradient descent step in for the objective f - 1n1 the next value of fa must lie closer to ã than Tn Assuming that, average. the gradient descent steps do not prefer certain direction, Øn-1 is more likely to lie in I {x E Rd L(x) < L(x,)} than outside it. The
region of increasing loss is shaded orange in Fig. 7 Since this region always fills less than half of the sphere around Tn. L is more likely to decrease than increase.
While this shows that the loss should decrease average, it does guarantee convergence. We now look at two specific formulations of P for which convergence to the correct solution is guar- anteed. The first involves the physical gradient pointing directly towards a solution x², This is, for example, the case for unimodal problems when we choose k large enough. In the second case, we consider the physical gradient being aligned with the gradient descent vector in x space.
Theorem 3. Vr Rd YE (0,1) P(x) + X(x - x), then the sequence fo., with
On+1 = On
01 (x - P(x)) converges to x².
T
20
T
Proof. Rewriting x - P
- yields the update on+1 0n = =-
aL of
de 30
1--==1 11. This describes pure gradient descent towards x" with the gradients scaled by A. Since A E (0, 1]. the convergence proof of gradient descent applies.
where L = ~ ||ı -
Theorem 4. If vx
E Rd YEE E [0 T
then the sequence fo. with (x - P()) converges to minima of L.
P(x) = x - (
On+1 = =On - 80
Proof.
This is equivalent to gradient descent in L(0). Rewriting the update yields 0x41 - 8, 1 ))) - Vu
-17 by X.
0f 88
x
(L 0f Oz 89
which is the gradient descent update scaled
B EXPERIMENTS
Here, we give a more detailed description of our experiments including setup and analysis. The implementation of our experiments is based on the &Row (PhiFlow) framework (Holl et al. (2020) and uses TensorFlow (Abadi et al., 2016) and PyTorch (Paszke et al. (017)) for automatic differentiation. Our code is open source and will be made available upon acceptance.
15

---

Under review as a conference paper at ICLR 2022
Loss
Eo MAE
A MAE
10
10²
10
Guessing
Guessing
Guessing
L-BFGS-B B
L-BFGS-B
L-BFGS-B
Network
Network
Network
10
10"
0
500
1000
0
500
1000
0
500
1000
lterations
Iterations
Iterations
Learning curves of the network trained to fit wave packets. The performance of L-BFGS- B and random guessing are shown on the same data for reference. The left graph shows the objective lly - y" 113 and the right two graphs show the x-space deviation from the true solution position &o and amplitude A.
Figure 8:
B. WAVE PACKET FIT
This experiment is instance of generic curve fitting problem. The task is to find the parameters that result in least mean squared error between two curves. We use this problem to illustrate the advantages of using neural networks inverse solvers over classical optimization algorithms.
Data generation. We simulate an observed time series y" from a random ground truth position and amplitude x- {ED, A}. Each time series contains 256 entries and is made up of the wave packet and noise. For the wave packet, we sample Ço E [25.6, 230.4) and A E [1,4) from uniform distributions. The wave packet has the functional form
v(t) = A sin(s - exp (- 1 (5 - {o)
d²
where we f= and 20 constant for all data. For the noise, we superimpose random values sampled from the normal distribution (0, )) at each sample.
Network architecture. We construct the neural network from convolutional blocks, followed by fully-connected layers, all using the ReLU activation function. The input is first processed by five blocks, each containing max pooling operation and two ID convolutions with kernel size 3. Each convolution outputs 16 feature maps. The downsampled result is then passed to two fully connected layers with 64 and 32 and 2 neurons, respectively, before third fully-connected layer produces the predicted En and A. Eo is passed through Sigmoid activation function and normalized to values between 25.6 and 230.4.
Training and fitting. We fit the data using L-BBFGS-B with an initial guess of &o = 128 and and the network output is offset by the same amount. Both network and L-BFGS-B minimize the squared loss [|y(Eo, A) - --y" y1 and the resulting performance curves are shown Fig. g We observe that L-BFGS-B manages to fit the wave packet when it is reasonably close to the center where the initial guess predicts it. When the wave packet is located to either side, L-BFGS-B B does not find and instead fits the noise near the center.
The neural network is trained using Adam with learning rate 0.001 and batch size 100. Despite the simpler, Hessian-free, optimization updates, the network quickly learns to localize all wave packets,
16

---

Under review as a conference paper at ICLR 2022
outperforming L-BFGS-B after 30 40 training iterations. This improvement is possible because of the network's reparameterization of the problem, allowing for joint parameter optimization using all data. When the prediction for one example is close to a local optimum, updates from different examples can prevent it from converging to that sub-optimal solution.
B.2 SINGLE-PARAMETE OPTIMIZATION
The optimization consists of finding x= argmin, sin( ()). Fig. shows the gradients for the case a = 1. The PG as well as the supervised methods require finding a solution x² closest to an initial guess To- We solve this analytically by first determining the index of the closest minimum,
2T which we round to the closest integer. The furthest-out minima (n = 0 and n = ===-) 11)
require special handling. Then the position of the minimum corresponding to is x 3//2425-
Newton's method
The analytic form of a Newton update is
cos(2(//2)
P(x)
2x cos(a/2) sin(a/x)
which approaches both minima and maxima. We alter this update to walk only towards minimum points by flipping the optimization direction so that the update points in the same direction as Be. Since Newton's 's method relies on an inversion of the Hessian, the produced physical gradients di- aP Be to [- before they are
OP
8 O. To avoid overflow, we clip the gradients verge where ले backpropagated through the neural network.
Inverse gradients In addition to gradient inversion via the Hessian, we can invert the gradients directly. The principle of inverting the Jacobian is not new. Inverse kinematics applications such as controlling robots also employ first-order inversion (Craig 2005). This results in a similar behavior concerning dimensions and function sensitivity Newton-type methods. However, this first-order inversion cannot directly be applied to general optimization tasks because the Jacobian of a scalar objective function is row vector. For this specific inverse problem, however, the gradient can be inverted because only a single parameter is optimized. The physical gradient is computed as
h
sin(a/s3 +1
P(x)-x=
-
VP
a cos(a/a)
where h = sin (9) + measures the height above y". It approaches too when sin is maximal. This inversion requires domain knowledge since h is based on the fact that min P = 1. Like with Newton's 's method, we apply gradient clipping to bound diverging inverse gradients.
Neural network training We set up a multilayer perceptron with a single input and output well as three hidden layers containing 10, 20 and 10 neurons, respectively. This adds up to 461 total trainable parameters. A bias is applied at each layer and the ReLU activation function is applied after each hidden layer. The network is trained using Adam with learning rate of 0.002 and mini- batches containing 1000 uniformly sampled [0.1, 10.1]. We FHow together with Py Torch to train the neural network in this experiment.
Learning curves for different network initializations are shown in Fig. 10 the recorded compu- tation times are shown in Fig. [4 While supervised training shows very consistent behavior across training runs, the convergence curves of Adam show more erratic behavior. Training with PGs based on the analytic solver converges in all cases, but does not always reach the same level of accuracy. While yielding superior accuracy than pure Adam for the large majority of tested initializations, sometimes stops converging to accuracy larger than 10-4. This is most likely caused by nu- merical effects of the inverse solver. Training with physical gradients based on Newton's 's method inverse gradients yields convergence that strongly depends on the network initialization. This variant converges only for certain initializations. However, upon convergence, typically reaches an accuracy similar the inverse physics variant.
B.3 POISSON'S EQUATION
We consider Poisson's equation, Viy = x where x is the initial state and y is the output of the sim- ulator. We up two-dimensional simulation with 80 by 60 cubic cells. Our simulator computes
17

---

Under review as a conference paper at ICLR 2022
0.50
0.25
0.00
Supervised
-0.50
0.50
0.25
0.00
-0.50
0.50
0.25
0.00
-0.50
0.50
0.25
0.00
-0.50
0.50
0.25
0.00
PG
-0.25
-0.50
1.0
0.5
0.0
sin(1
10"
10"
10
Gradients for minimizing sin(a/x), shown for 1. From top bottom: Supervised learning with Tp = 1, Gradient descent, inverse gradient, Newton's 's method, physical gradient, objective function. Singularities for small not properly resolved in the IG / Newton plots.
Figure 9:
18

---

Under review as a conference paper at ICLR 2022
Seed
10
10"
10
Newton
I
Supervised
10 4
2500
5000
7500
10000
12500
15000
17500
20000
Seed
10
10
PG
GD
10
Newton
Supervised
2500
5000
7500
10000
12500
15000
17500
20000
Seed
10
10"
GD
10
Newton
IG
Supervised
10 1
2500
5000
7500
10000
12500
15000
17500
20000
Seed
10
10
PG
GD
10
Newton
Supervised
10
2500
5000
7500
10000
12500
15000
17500
20000
Networks trained on single-parameter optimization using Adam and various gradient schemes. Each figure shows the learning curves for network initialized with fixed seed between 5 and The x axis denotes the number of training iterations, Solid lines show the running average 1000 mini-batches.
Figure 10:
19

---

Under review as a conference paper at ICLR 2022
3/ P(x) = V implicitly via the conjugate gradient method. The inverse problem consists of finding an initial value x² for a given target y" such that V2y" y" = =x" We formulate this problem as minimizing L(x) p1) v²1 3 3xx
x²) ll.-
We now investigate the computed
updates Ar of various optimization methods for this problem.
Gradient descent Gradient descent prescribes the update Ar = -1j 8t.
= - Vy² y")
(LL
which requires additional implicit solve for each optimization step. This backward solve produces much larger values than the forward solve, causing GD-based methods to diverge from oscillations unless 17 very small. We found that GD requires 17 S 2 10-5, while the momentum in Adam allows for larger 7). For both GD and Adam, the optimization converges extremely slowly, making GD-based methods unfeasible for this problem.
Physical gradients via analytic inversion Poisson's equation can easily be inverted analytically, Oz
yielding x = V2y. Correspondingly, we formulate the update step as Az -1 By (v y²) =
-7 2 (y- 0 y*) which directly points to for n 1. Here the Laplace operator appears in the computation of the optimization direction. This is much easier to compute numerically than the Poisson operator used by gradient descent. Consequently, no additional implicit solve required for the optimization and the cost per iteration is less than with gradient descent. This computational advantages also carries over to neural network training where this method can be integrated into the backpropagation pipeline as a physical gradient.
Neural network training We first generate ground truth solutions x* by adding fluctuations of varying frequencies with random amplitudes. From these x", we compute y" = P(x") to form the set of target states V. Both generation of x- and y" is performed on the fly, resulting in data set of effectively infinite size, | = oc. This has the advantage that learning curves are representative of both test performance well as training performance. The top of Fig- 1T shows some examples generated this way. We train a U-net (Ronneberger et al.] (015)) with a total of 4 resolution levels and skip connections. The network receives the feature map y as input. Max pooling is used for downsampling and bilinear interpolation for upsampling, After each downsampling or upsampling operation, two blocks consisting of 2D convolution with kernel size of 3x3, batch normalization and ReLU activation are performed. All of these convolutions output 16 feature maps and a final 1x1 convolution brings the output down to one feature map. The network contains a total of 37,697 trainable parameters.
For SGD and Adam training, the composite gradient of NN P is computed with TensorFlow or PyTorch, enabling an end-to-end optimization. The learning rate is set to 1) = 10-3 with Adam and ? = =110- for SGD. The extremely small learning rate for SGD is required to balance out the large gradients and is consistent with the behavior of gradient-descent optimization on single examples where an 77 = 2 10 was required. We use typical value of 0.9 for the momentum of SGD and Adam. For the training using Adam with physical gradients, we compute Ar as described above and keep 1/ = 10-3 For each case, we set the learning rate to the maximum value that consistently converges. The learning curves for three additional random network initializations shown the bottom of Fig. [1 while Fig. 14|shows the computation time per iteration.
B.4 HEAT EQUATION
du V- V Given an We consider two-dimensional system governed by the heat equation 57 initial state x = Ug at to. the simulator computes the state at a later time t. via = u. = P(r). Exactly inverting this system is only possible for t v larger t. because initially distinct heat levels even out over time, drowning the original information in noise. Hence the Jacobian of the physics 8h a domain consisting of 64x64 cells of unit length. This level of diffusion challenging, and diffuses most details while leaving the large-scale structure intact.
and becomes increasingly unstable for
dy is singuai.. In our experiment we set t- = 8 on
We apply periodic boundary conditions and compute the result in frequency space where the physics can be computed analytically as 9 = â e ---((t- where f = F(y), denotes the k-th element of the trunfeeetransformmed vector y. Here, high frequencies are dampened exponentially. The inverse
20

---

Under review as a conference paper at ICLR 2022
Adam y
Adam+PG y
Adam x
Adam+PG x
x
Adam
Adam+PG
Adam x
Adam+ x
x
Adam y
Adam+PG
Adam x
Adam+ x
x*
10a
10a
10a
10a
Adam + PG
10a
Adam PG
10a
Adam PG
Adam
Adam
Adam
0
5000
10000
5000
10000
5000
10000
Figure 11: Inverse problems involving Poisson's 's equation. Top: Three examples from the data set, from left to right: observed target (y"), simulated observations resulting from network predictions (Adam y. Adam+PG y). predicted solutions (Adam æ, Adam+PG æ), ground truth solution (x*), Networks were trained for 12k iterations. Bottom: Neural network learning curves for three random network initializations, measured as || - x ll1-
21

---

Under review as a conference paper at ICLR 2022
problem can thus be written as minimizing L(x) = ||P(x) - v011 = ||= (F(x)) -e Ke -- e
-
y" 2
Gradient descent
Using the analytic formulation, we can compute the gradient descent update as Ar = -1 - ()
GD applies the forward physics to the gradient vector itself, which results in updates that are stable but lack high frequency spatial information. Consequently, GD-based optimization methods con- verge slowly on this task after fitting the coarse structure and have severe problems in recovering high-frequency details. This is not because the information fundamentally missing but because GD cannot adequately process high-frequency details.
Stable physical gradients ically, yielding t. Úx e² (t.--t)) This allows to define the update
The frequency formulation of the heat equation can be inverted analyt-
Ar --=-n-F F-1 (cH F(y- (
Here, high frequencies are multiplied by exponentially large factors, resulting in numerical instabil- ities. When applying this formula directly to the gradients, it can lead to large oscillations in Ar. This is the opposite behavior compared to Poisson's 's equation where the GD updates were unstable and the PG stable.
The numerical instabilities here can, however, be avoided by taking a probabilistic viewpoint. The observed values y contain certain amount of noise n, with the remainder constituting the signal s = y n. For the noise, we assume a normal distribution n ~ (N y) with c > and for the signal, we assume that it arises from reasonable values of a so that (N, 5 -e e-*²) with 5 > 0. With this, we can estimate the probability of an observed value arising from the signal using Bayes' theorem p(s/v) p(plss)ppss)))ppilaappn where we assume the priors p(s) = p(n)
1 Based on this probability, we dampen the amplification of the inverse physics which yields stable inverse. Gradients computed this way hold as much high-frequency information as can be extracted given the noise that is present. This leads to much faster convergence and more precise solution than any generic optimization method.
p(s(s)-p(s)
Neural network training For training, we generate x by randomly placing between and 10 hot rectangles of random size and shape in the domain and computing are executed on the fly so that Y = oo and the learning curves are representative both training and test performance. For the neural network, we use the same U-net architecture in the previous experiment. We train with a batch size of 128 and a constant learning rate of 1) 10 3 The network updates are computed with TensorFlow's or PyTorch's automatic differentiation. Fig. T2 shows two examples from the data set, along with the corresponding inferred solutions, as well as the network learning curves for two network initializations. The measured computation time per iteration is shown in Fig. 14
P(x*). Both operations
B.5 NAVIER-STOKES EQUATIONS
Here, we give additional details on the simulation, data generation, physical gradients and network training procedure for the fluid experiment.
Simulation details We simulate the fluid dynamics using direct numerical solver. We adopt the marker-in-cell (MAC) method (Harlow & Welch 1965) 1965) Harlow 1972) which guarantees stable sim- ulations even for large velocities or time increments. The velocity vectors are sampled in staggered form at the face centers of grid cells while the marker density is sampled at the cell centers. The initial velocity Up is specified at cell centers and resampled to a staggered grid for the simulation. Our simulation employs second-order advection scheme (Selle al. 2008) to transport both the marker and the velocity vectors. This step introduces sigmificant amount of numerical diffusion which clearly be seen in the final marker distributions. Hence, do not numerically solve
22

---

Under review as a conference paper at ICLR 2022
GD
BFGS
Inv.phys.
Adam
A.+ BFGS
A.+PG
GD
BFGS
Inv. Inv.phys.
Adam
A.+ BFGS
A.+PG
1
103
Adam PG (Inv.phys.)
10a
Adam + (Inv.phys.)
Adam
Adam
Adam PG (L-BFGS-B)
Adam + PG (L-BFGS-B)
10a
102
2000
4000
6000
8000
10000
2000
4000
6000
8000
10000
Inverse problems involving the heat equation. Top: Two examples from the data set. The top row shows observed target (y") and simulated observations resulting from inferred solutions. The bottom row shows the ground truth solution (x*) and inferred solutions. From left to right: ground truth; gradient descent (GD), L-BFGS-B B (BFGS) and inverse physics (Inv.phys.). running for 100 iterations each, starting with 2g O; Networks trained for 10k iterations. Bottom: Neural network learning curves for two random network initializations, measured in terms of || x*I|i-
Figure 12:
23

---

Under review as a conference paper at ICLR 2022
for adding additional viscosity. Incompressibility is achieved via Helmholz decomposition of the velocity field using a conjugate gradient solve.
Neither pressure projection nor advection energy-conserving operations. While specialized energy-conserving simulation schemes for fluids cxist (Gluhovsky & Tong 1999, Mullen et al.
2009). instead enforce energy conservation by normalizing the velocity field at each time step the total energy of the previous time step. Here, the energy is computed as E = -gg= dx u(z) since assume constant fluid density.
Data generation The data set consists of marker pairs (mo. me) which are randomly generated on-the-fly. For each example, center position for mg chosen on grid of 64x64 cells. mg is then generated from discretized noise fluctuations to fill half the domain size in each dimension. The number marked cells is random.
Next, a ground truth initial velocity is generated from three components. First, uniform velocity field moves the marker towards the center of the domain to avoid boundary collisions. Second, large vortex with random strength and direction is added. The velocity magnitude of the vortex falls off with a Gaussian function depending on the distance from the vortex center. Third, smaller-scale vortices of random strengths and sizes added additionally perturb the flow fields. These generated by assigning a random amplitude and phase to each frequency making up the velocity field. The range from which the amplitudes are sampled depends on the magnitude frequency.
Given mi and to, a ground truth simulation is run for = 2 with AL = 0.25. The resulting marker density is then used the target for the optimization. This ensures that there exists solution for each example.
Computation of physical gradients To compute the physical gradients for this example, we con- struct an explicit formulation in = P- pm- my (o) that produces estimate for given an initial guess zo by locally inverting the physics. From this information, fits the coarse velocity, ie. the uniform velocity and the vortex present in the data. This use of domain knowledge, i.e., enforeing the translation and rotation components of the velocity field as a prior, is what allows to produce a much better estimate of vp than the regular gradient. More formally, it assumes that the solution lies on a manifold that is much more low-dimensional than do- On the other hand, this estimator ignores the small-scale velocity fluctuations which limits the accuracy it can achieve. However, the difficulty of fitting the full velocity field without any assumptions outweighs this lim- itation. Nevertheless, GD could eventually lead to better results if trained for an extremely long time.
To estimate the vortex strength, the estimator runs a reverse Navier-Stokes simulation. The reverse simulation is initialized with the marker m. = m, and velocity vi uf from the forward simula- tion. The reverse simulation then computes mov and yor for all time steps by performing simulation steps with At = -0.25. Then, the update to the vortex strength is computed from the differences no m at each time step and an estimate of the vortex location these time steps.
Neural network training We train a U-net (Ronneberger et al. 2015) similar the previous periments but with 5 resolution levels. The network contains a total of 49,570 trainable parameters. The network is given the observed markers my and my resulting in input consisting of two fea- ture maps. It outputs two feature maps which are interpreted as a velocity field sampled at cell centers.
The objective function is defined as [F(P(x) y²)l w where F denotes the two-dimensional Fourier transform and w is weighting vector that factors high frequencies exponentially less than low frequencies.
We train the network using Adam with learning rate 0.005 and mini-batches containing 64 ex- amples each, using Porch''' automatic differentiation to compute the weight updates. We found that second-order optimizers like L-BFGS-B yield no significant advantage over gradient descent, and typically overshoot in terms of fregueryy motions. Example trajectories and reconstruc- tions are shown Fig. 13 and performance measurements are shown Fig. 14
24

---

Under review as a conference paper at ICLR 2022
GT
A+PG
Adam
32
GT
N
A+PG
N
Adam
GT
A+PG
Adam
Three example inverse problems involving the Navier-Stokes equations. For each exam- ple, the ground truth (GT) and neural network reconstructions using Adam with physical gradient (A+PG) and pure Adam training (Adam) are displayed as rows. Each row shows the initial veloc- ity Do z well five frames from the resulting marker density sequence m(t), at time steps t E [0, 0.5, 1, 5, 2}. The differences the Adam version are especially clear in terms of up-
Figure 13:
25

---

Under review as a conference paper at ICLR 2022
Sirgle-parameter Optimization (Sin)
Sec. / /iteration
0.006
Adam + PG (Inv.phys.)
Adam
Adam PG (Newton)
0.005
Adam + PG (Inv.grad.)
Adam (supervised)
2500
5000
7500
10000 Iteration
12500
15000
17500
20000
Poisson Equation
Sec. /iteration
1.0
Adam PG (Inv.phys.) Adam
0.5
2000
4000
6000 Iteration
8000
10000
12000
Heat Equation
Sec. /iteration
0.14
Adam + PG (Inv.phys.)
Adam
0.12
0.10
2000
4000
6000
8000
10000
Iteration
Incompressible Fluid
4.50
Sec. Itteaatoon
Adam + PG (Inv.phys.)
4.25
Adam
Msaa Aamal
4.00
2000
4000
6000
8000
10000
Iteration
Figure 14: Measured time per neural network training iteration for all experiments, averaged over 64 mini-batches. Step times were measured using Python's 's perf.counter () function and in- clude data generation, gradient evaluation and network update. In all experiments, the computa- tional cost difference between the various gradients marginal, affecting the overall training time by less than 10%.
26


Reviews:
[{"review_id": "qniXCfzVveo", "reviewer": "Reviewer_ML7L", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The omnipresence of inverse problems asks for stable and efficient deep learning-based solvers and new methods are highly desired.\nWhile the proposed usage of so-called physical gradients is based on a rather minor modification of the supervised approach, it could overcome several drawbacks of existing methods and thus seems to be a research direction well worth pursuing. \nUnfortunately, the present paper does not sufficiently motivate and justify the usage of physical gradients and has several other issues.\n\n\n- **Motivation:**\n\n    The proposed method employs an iterative solver in each optimization step in order to create the current targets for the neural network. Thus, it seems that the neural network cannot outperform the iterative solver except for perhaps averaging out noise or adapting the initialization. This raises the question, why one should prefer the neural network over the iterative solver?  Similar to the supervised approach, it might be beneficial to replace the iterative solver by a neural network, as this could lead to faster inference times in case of a computationally expensive solver and allow for sensitivity analysis / uncertainty estimation via the trained neural network. The former, however, would drastically slow down the neural network training when using physical gradients as the iterative solver needs to be applied in every optimization step and the latter is not explored in the paper. \n\n\n- **Empirical evidence:**\n\n    Unfortunately, the numerical experiments do not showcase the potential strengths of the proposed method. While a crucial part of this method is the initialization of the iterative solver using the current network prediction, this seems to only be used in the first toy example. In the other examples, the difference between the use of physical gradients and the supervised approach is not clear. Moreover, the proposed method might also work with only a few iterations of the iterative solver in each step. However, for all but the last example, an explicit solution is used instead of an iterative solver and there is no numerical evidence of how the method depends on the accuracy of the iterative solver.  Finally, all experiments use an infinite stream of synthetically created (mostly noiseless) data, which is usually not given in practice. \n\n\n- **Further comments:**\n\n    Further issues, comments, and suggestions for improvements can be found in the following:\n     \n    (i) The statement that physical gradients behave more smoothly, both in magnitude and direction, appears in Subsection 3.4 and in the abstract. However, apart from the toy example in Figure 7, it seems that there is not enough theoretical or empirical evidence for that statement. \n    \n    (ii) In Equations (1) and (2), $x^*$ and $\\theta_*$ belong to the $\\operatorname{argmin}$ (which is in general a set of minimizers, given that the minimum is attained), rather than being equal to the $\\min$.\n    \n    (iii) I guess that $\\mathcal{P}^{-1}_{n^*}(y^* | x_0)$ should, in general, not be equal to $x^*$ as stated in the introduction of Subsection 2.1. It should probably only be an approximation of the latter.\n    \n    (iv) In my opinion the sentences \"Note that the dimensions of ...\" and the subsequent one, are rather confusing. What is meant by \"dimensions\" in this context, why is this phenomenon counter-intuitive, and how does this relate to the vanishing/exploding gradient problem in *deep* networks?\n    \n    (v) It seems that the derivative of the forward operator is missing in the update rules in Figure 1 and also the expression for the Hessian given in the caption should be derived in more detail.\n    \n    (vi) I am not aware of a generic definition of a \"sensitive\" function (as mentioned several times in the paper) and think that this should be properly defined.\n    \n    (vii) The statement that jointly optimizing all inverse problems reduces the likelihood that an individual solution gets stuck in a local minimum, should, if it holds true in this generality, be referenced or sufficiently motivated.\n    \n    (viii) In Algorithm 1, it is unclear what it means to condition the inverse solver on two values $x_{k-1}$ and $y_k$.\n    \n    (ix) The first example does not fit into the framework described in Equations (1) and (2) as now $\\mathcal{P}$ is parametrized by $a$ and $y^*$ is fixed.\n    \n    (x) In the last example, $y^*$ seems to correspond to the tuple $(m_0,m_t)$ rather than only $m_t$.\n    \n    (xi) It is not clear, how Assumption 3 on page 15 relates to the universal approximation theorem (as stated below), which deals with uniform approximation of a continuous function on a compactum rather than convergence of gradient descent.\n    \n    (xii) The factor $1/2$ seems to be missing in the loss $L$ in Subsection B.2.\n    \n    (xiii) In Subsection B.4 there appears $x_0$ which seems not to be defined.\n    \n    (xiv) Typos: *an* generic inverse solver; *we* found that second-order optimizers", "summary_of_the_paper": "The present paper proposes a novel ansatz for solving inverse problems by means of deep learning. Specifically, the loss attached to a specific observation is defined as the difference between the neural network prediction and the output of an approximate iterative solver initialized with the current neural network prediction. The gradient of such loss (for a mini-batch of observations) w.r.t. to the neural network parameters, coined physical gradient, can then be used to optimize the neural network by standard variants of gradient descent. Potential advantages of this approach over the classical (un-)supervised approaches include the following: (1) One circumvents the need for (possibly ill-conditioned) automatic differentiation of the forward operator, (2) the varying initialization of the iterative solver enables adaptive convergence to different solutions (in case of a non-injective forward operator) throughout the neural network training, and (3) a few iterations of the solver in each step might suffice to provide a useful loss for the neural network. Moreover, the paper also explores the use of (Quasi-)Newton methods.  \nThe functionality of the proposed method is backed by theoretical results and four numerical experiments.", "summary_of_the_review": "The paper proposes an interesting adaptation to the supervised approach of solving inverse problems by means of deep learning. Unfortunately, the included experiments are not representative enough to show whether the promised benefits of this approach hold true. In addition, there are several issues regarding the presentation of the material.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "2: Several of the paper’s claims are incorrect or not well-supported.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "A815YhxJ0fK", "reviewer": "Reviewer_mue1", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 2, "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": "Frankly speaking, I did not understand the main goal of this paper. I am happy to adjust my scores when the concerns are addressed and there is a high chance that I did not understand vital parts of this paper.\n\nThe main method or the goal of the paper is illustrated page 2 Section 2 with equations (1) and (2). The inverse problems are usually framed into two different frameworks:\n1. Optimization framework where x* = min_x L(x). The goal is to minimize the loss function which is defined as the difference between the target y and F(x) (F is the forward physical process). This is consistent with the paper's equation (1). However, this is not what this physics gradient is meant for.\n2. Mapping learning framework where we learn (usually using NN) either a probabilistic mapping or deterministic mapping from target y space to the control parameter space x. That is we learn Inv: Y -> X. This seems to be what the authors are trying to illustrate from Equation (2). However I struggle to understand the notion of \"This training scheme is unsupervised, i.e. no labels are required over the training data Y\". How can the mapping paradigm be learned without any training data pairs (x, y) which comes from the actual physical process? Do the authors mean the P(NN(y*|theta)) can somehow skip being evaluated by the actual physical process? If not, which means the underlying physical process is called during training, then we have the ground truth pair of (x,y) during training, which means this counts as an supervised training scheme.\n\nUnless the authors are referring to some more advanced techniques like Physics informed neural network where the PDE equation prior is hardcoded in the loss, I really can't see why this is a unsupervised task.\n\nSecondly, I didn't quite understand the assumption of having an inverse problem solver P^{-1} in the first place. If we have a numerical inverse problem solver, why do we need to use a neural network to solve the inverse problem anyway?  The assumption of having access to the forward numerical process P is ok, but the access to an oracle P^{-1} really made me struggle to understand. \n\nThis is the biggest concern I have, this method proposes to solve inverse problems using neural networks, but it needs access to a numerical inverse solver from domain knowledge, which seems to beats the whole purpose of using neural networks.\n\nFinally, the empirical results. Given the above questions raised a mere misunderstanding on my part, all the \"convergence speed comparison\" is done on the unit of \"iteration\". However, I would assume that the numerical inverse solver P^{-1} would take a much longer time in wall clock sense than a neural network update. Therefore without this critical information about the actual running time of the inverse solver, the mere comparison between performance upon iterations is meaningless as the total time of Inv. Phys gradient can well exceed the Adam and even Adam + Newton.\n\nMy apology if I missed crucial parts of this paper, but I did spend some time on it and still struggle to understand things I mentioned above. If the above concerns are ultimately addressed by the rebuttal, the clarity of the delivery would need a lot of improvement.", "summary_of_the_paper": "The paper proposed a new update method for neural network solving inverse problem by incorporating the physical information as prior into the gradient computation step. It leads to faster convergence and better convergence performance.", "summary_of_the_review": "The basis for problem set up is unclear to me in general, therefore the contributions and results are completely compromised by that.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "2: Several of the paper’s claims are incorrect or not well-supported.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "8spuWh6jocM", "reviewer": "Reviewer_FRkr", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Strengths:\n* The idea the paper proposes is simple but powerful: That writing down a particular (intuitive) proxy L2 loss for the original inverse problem yields stable gradient updates with direct physical intuition in the physical output space.\n* The empirical demonstrations are extensive, and demonstrate the superior performance of the proposed method in physical settings with different salient properties (e.g., multi-modality, large sensitivity, singularities, etc.) \n\nWeaknesses:\n* The paper is written so as to significantly oversell the generality of the proposed framework, which in turn significantly impacts clarity. There are several major examples of this:\n\t* The title \"physical gradients\" and several phrases throughout (e.g. in the abstract: \"We replace the gradient of the physical process by a new construct\") strongly imply that a new method of taking gradients will be introduced. This is not the case. Instead, the original inverse problem is rewitten such at the gradients that \"fall out\" have physical intuition.\n    * The formulations in Equations 4 and 5 depend on $\\mathcal{P}_n^{-1}$, which in various places is described as being derived using domain knowledge. However, $\\mathcal{P}_n^{-1}$ is just a single step of gradient descent, and while this is mentioned, it is not emphasized. Notably, all experiments use this \"single step of gradient descent\" interpretation, and the physical interpretation of the resultant gradient depends on this as well. This point should be made much clearer.\n    * Relatedly, the abstract indicates that the training approach \"combines higher-order optimization methods with machine learning techniques.\" To the best of my knowledge, this does not seem to be true.\n* Certain claims made in the paper are not fully substantiated. For instance:\n    * The claims in Figure 1(b) about the strengths of the present method vs. GD and Newton's method are not fully substantiated/described.\n    * On page 4, it is noted that the method can be applied under one of three conditions. However, only condition (iii) is ever explicitly used (as far as I understand), and the rationale as to why either of the other two conditions is sufficient is never described.\n* The main body of the paper is not self-contained, causing issues with clarity that are then exacerbated by typos. In particular, Appendix A.1 should definitely be in the main body of the paper in order to make the overall method/intuition understandable. (Also, the last paragraph of page 4 mentions that x_n can be taken as a constant - this is, importantly, not the case, but requires digging into Appendix A.1 in order to understand.)\n\nMinor comments (not affecting my score):\n* The paragraph on \"inverse physics via domain knowledge\" is hard to understand until after the entire paper has been read; its clarity could be improved.\n* Above equation 3: The authors might consider using a different subscript than \"sup.\" Presumably this is meant to stand for \"supervised,\" but could be confused for \"supremum\" given the optimization setting.\n* In the experiments, the descriptions of what methods are being compared could be made much crisper and clearer. In addition, interpretation of the legend in e.g. Figure 2 should be made more self-contained.\n* In all graphs, legends should be moved such that the lines are visible.\n* In Figure 3, it is not clear what the x-axis is. The caption of 3(a) presumably has a typo, as convergence curves are not shown.", "summary_of_the_paper": "This paper proposes an approach for learning neural network-based approximators for inverse problems in physical domains. The crux of this approach involves rewriting the training loss function to capture the difference between (a) the NN's estimate of the problem input and (b) the solution obtained when using this estimate as an initial point for an approximate inverse problem solver. This loss function is then optimized via gradient descent. This approach is demonstrated on several physical systems, and shown to avoid issues due to multi-modality (which supervised training schemes fall prey to), instability in the gradients through the physical process (faced in some settings by first-order training methods applied to the \"standard\" inverse problem), and issues such as expense of computation or ill-conditioning (faced in some settings by training method using e.g. inverse Hessian information).", "summary_of_the_review": "The method provided by this paper is simple, powerful, and physically intuitive. However, there are significant issues with overselling, substantiation of claims, and clarity that mean I cannot recommend acceptance of this paper in its current form.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "2: Several of the paper’s claims are incorrect or not well-supported.", "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": "33hHRHsc0S", "reviewer": "Reviewer_pjrC", "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": "This work is clearly presented and proposes an interesting idea: using inverse solvers to provide targets that iteratively improve the neural network predictions of the initial state.\n\nThe case n=1 used is probably the most interesting, because the inverse step is indeed simply providing a slightly improved target to supervise the training of the network, without the need of labeled data and without having to run a full inverse solve.\n\nThe experiments performed, though on fairly simple domains, demonstrate that the method works and is able to generalize at least in the given domains (ie interpolate), given that data is sampled randomly.\n\nAs with many of these types of applications of machine learning to physical inference, it is questionable if the benefit of the fully trained model achieved at the end compensates the cost of training it for a long time on a limited domain, while also having to provide an already existing inverse solver (or something capable of approximating this inverse). It would be good for the authors could discuss further what use cases they foresee for this type of model in practice, where the costs of training and formulating a known inverse could outweight simply using a traditional methods. For the cases where the inverse is actually known, it would be important to provide a comparison to simply providing the inverse optimization.\n\n-------\n\nAfter reading the authors response, I maintain my previous evaluation.", "summary_of_the_paper": "This paper looks at the (inverse) problem of finding the some initial state given a final state of a physical system. It proposes training neural networks to predict the initial state for a given final state. These networks are trained iteratively by using inverse solvers (based on domain knowledge) that provide physical \"targets\" that are closer to the desired initial state than the current prediction. ", "summary_of_the_review": "Despite some of the comments above, I believe the conceptual contribution and the experiments presented in this paper are of enough interest and robustness to warrant its acceptance. I am rating it as marginally above the acceptance threshold for now, expecting the authors to improve and respond to the comments above.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "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:

Meta-review: The paper introduces a method to solve inverse problems: given y, find x such that P(x)=y, for a given physical simulator P. A standard approach is to learn a neural net such that the inverse x=NN(y;\theta). The authors state that this is problematic because it is difficult to take "higher order" gradient information into consideration when using this standard approach. The method assumes that there is an approximate inverse solver inv(P) and discusses an alternative "Physical Gradient" objective that can incorporate knowledge of an approximate inv(P) and a neural network. The experiments are good though comparing performance on an iteration basis is not always fair since an iteration of the PG method can be much more expensive than standard approaches.

The biggest issue that reviewers had was the clarity of the presentation. The authors have made a reasonable attempt to correct this, but I'm inclined to agree with the general reviewer sentiment that the presentation is still not at the required level. I agree that there are many things that are not clear, including the confusing discussion in section 2.1 about how the method takes higher order information into consideration. It only becomes partially transparent later in the experiments what is meant by higher order information.

Overall, I feel this is the basis of a potentially valuable contribution but that the current presentation is quite confusing. As mentioned by others, I would also suggest to find a different name since Physical Gradient is also rather misleading. 

The following points were not part of the review process and I do not base the final decision on them, but the authors may want to consider the following: 

I believe there is also an error in the basic approach, or at least an approximation is made which is not explained. The error is that the approximate inv(P) depends also on the parameter \theta (since this is used to initialise inv(P)). This dependency is not taken into consideration in the paper. For example, in theorem 1 in appendix A, the calculation of the gradient dM/d\theta is incorrect since the authors assume that inv(P) is independent of \theta, which it is not (since the preconditioner value depends on \theta). If we do take this into account, we would need to know the derivative of inv(P|x) with respect to the preconditioner x. This dependency would alter the gradient, potentially considerably. The gradient in figure 2 for the PG is also incorrect. One may of course simply say that the paper discounts this correction term in order to retain tractability; however, this would need to be stated as an approximation.
Decision: Reject