Year: 2021

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 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 take updates of a higher-order or domain-specific optimizer and embed them into the gradient-descent-based learning pipeline as a physical gradient, replacing the regular gradient of the physical process. 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.

heading: 1 INTRODUCTION
text: 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 et 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 & Huerta, 2018b), 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., 2020), 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-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 is far away from the true solution in parameter space, iterative methods may diverge or converge to a 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 or longer. However, many experiments continuously or 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 a 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 & Huerta, 2018a;b).
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 a different example. Fig. 1 shows a generic case illustrating this effect: The task is to localize a wave packetA·sin(f ·x)·exp(− 12 (x−x0)
2/σ2) and determine its amplitude from a noisy recorded time series (grey curve). Iterative solvers fail to reconstruct the wave packet
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 are given in appendix B.1. In this paper, we derive a 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 a more efficient manner. Our method computes the weight updates using both backpropagation as well as the embedded inverseproblem solver. Specifically, the first-order gradient of the physical process used in backpropagation is replaced by a physical gradient (PG) derived from the update of the inverse-problem solver. When used in conjunction with training deep learning models, physical gradients allow each update step to encode nonlinear physics, which yields substantially more accurate optimization directions. A key advantage of our approach is its compatibility with acceleration schemes (Duchi et al., 2011; Kingma & Ba, 2015) and stabilization techniques (Ioffe & Szegedy, 2015; Huang et al., 2016; Ba et 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 firstorder updates (Rackauckas et al., 2020).
We test our method on a wide variety of inverse problems including the highly challenging Navier-Stokes equations. In all of our experiments, neural networks trained with PGs significantly outperform state-of-the-art training schemes, yielding accuracy gains of up to three orders of magnitude.

heading: 2 METHOD
text: We consider unconstrained inverse problems that involve a differentiable physical process P : X ⊂ Rdx → Y ⊂ Rdy 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∗ ∈ Y , the inverse problem consists of finding optimal parameters
x∗ = arg min x
L(x) with L(x) = 1
2 ‖P(x)− y∗‖22. (1)
We are interested in approximating this problem using a neural network, x∗ = NN(y∗ | θ), parameterized by θ. Let Y∗ = {y∗i | i = 1, ..., N} denote a set of N inverse problems, all sharing the same physical process P . Then training the network is equivalent to finding
θ∗ = arg min θ N∑ i=1 1 2 ‖P ( NN(y∗i | θ) ) − y∗i ‖22 . (2)
This training scheme is unsupervised, i.e. no labels are required over the training data Y .

heading: 2.1 ITERATIVE OPTIMIZATION
text: Inverse problems as in Eq. 1 are classically solved by starting with an initial guess x0 and iteratively applying updates ∆xk. This yields an estimate of the inverse of the physical model P , i.e., x∗ ≈ P−1n∗ (y∗ |x0) where n∗ denotes the number of iterations required to reach a chosen convergence threshold. A key decision in such an iterative optimization scheme is how to derive the update ∆xk. Figure 2a visually compares gradient descent, Newton’s method and PGs which we will introduce in section 2.2.
First-order methods The scale of typical machine learning datasets necessitates the use of minibatches, which has led popular optimization methods in machine learning to rely only on firstorder derivatives of L. These Hessian-free methods are popular in deep learning due to 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 ∆x = −η· ( ∂L ∂x )T . Note that ∂L∂x , and therefore of ∆x, scales inversely with x. When x carries physical units, neither ∆x nor ∆xT lie in the same space as x. This scaling behavior is a problem when dealing with sensitive or insensitive functions, i.e. functions whose gradient magnitudes ∣∣∂L ∂x
∣∣ are far from 1. When small changes in x cause large changes in L, this counter-intuitively leads to large updates ∆x 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 (Ioffe & Szegedy, 2015; Loshchilov & Hutter, 2019) 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.
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 et 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 method (Atkinson, 2008) uses the inverse Hessian to determine the optimization update ∆x = −η · ( ∂2L ∂x2 )−1 ( ∂L ∂x )T . Direct computation of the Hessian adds a significant computational cost to 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 insensitive 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., 2007) and the high cost of evaluating the Hessian directly have largely prevented second-order methods from being utilized in machine learning research (Goodfellow et al., 2016).
Domain-specific inverse solvers Domain knowledge can often be used to formulate more efficient optimization updates or to directly approximate a 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
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 generic inverse solver. These methods are commonly used in optimization but their use has eluded machine learning applications so far due to the fixed gradient-descent pipeline. When learning with PGs, domain-specific optimizers can be integrated into training like any other optimizers, allowing for greatly-increased training convergence. Appendix appendix A.2 lists necessary conditions that need to be fulfilled.

heading: 2.2 PHYSICAL GRADIENTS FOR DEEP LEARNING
text: 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-problem solvers as physical gradients (PGs) into the first-order optimization pipeline. The network updates ∆θ can then be computed using any first-order optimizer such as SGD or Adam. Additionally, all state-of-the-art deep learning techniques, such as normalization Ioffe & Szegedy (2015); Ba et al. (2016) or dropout Srivastava et al. (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 is available: (i) ∂ 2P ∂x2 can be computed numerically, (ii) a higher-order update such as ( ∂2P ∂x2 )−1 ∂P ∂x 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 access to better updates ∆x than ∆xGD ∝ ∂P∂x . These ∆x are the basis for our method which we describe in the following.
We derive our method starting from Eq. 2 which describes the combined problem of finding both solutions x∗ to the individual inverse problems y∗ ∈ Y∗ and network weights θ that approximate them. One way to solve this problem is to first consider all inverse problems separately and precompute corresponding solutions X sv = {xsvi : P(xsvi ) = y∗i } using an iterative optimizer. Then, X sv could be used as labels for supervised training of a neural network:
θ∗ = arg min θ N∑ i=1 1 2 ‖NN(y∗i | θ)− xsvi ‖22. (3)
This precomputation makes it possible to use an efficient optimization method forX sv 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 (Ioffe & Szegedy, 2015; Ba et al., 2016) or dropout (Srivastava et al., 2014). However, this approach has severe drawbacks. For one, the individual optimizations 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.e. multiple solutions x∗ exist for one y∗. In such settings xsv heavily depends on the initial guess x0 used for precomputation, which may cause the network to interpolate between possible solutions, leading to subpar convergence and generalization performance.
To avoid these problems, we alter the training procedure from Eq. 3 in two ways. First, we treat the inverse problem solver P−1 as part of a coupled optimization and consider it as part of the training loop, yielding
θ∗ = arg min θ N∑ i=1 1 2 ‖NN(y∗i | θ)− P−1n∗i (y ∗ i )‖22 . (4)
Next we can condition the embedded inverse problem solver on the neural network prediction, P−1n∗i (y ∗ i ) → P −1 n∗i
(y∗i |NN(y∗i | θ)), by using it as an initial guess. The embedded optimizer can now find a minimum close to the prediction, which allows the neural network to choose any solution 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 a local minimum.
It is furthermore not necessary to run the iterative inverse solver to convergence. Any choice of n ∈ N steps likely leads to convergence (proof given in appendix A.2), and choosing small n can
increase the training speed. We also observe this in practice: Optimizations with n = 1 converged reliably in all our experiments. With these modifications, we arrive at
θ∗ = arg min θ N∑ i=1 1 2 ‖NN(y∗i | θ)− P−1n (y∗i |NN(y∗i | θ))‖22. (5)
The central quantity here is the difference between prediction and correction obtained from P−1n , which we refer to as the physical gradient. It encodes a valid physical state and takes the place of( ∂P ∂x )T , which would otherwise be computed by backpropagation.
Since the PG stems from a generic higher-order or a domain-specific solver, the resulting updates ∆θ will also be non-linear without computing the Hessian w.r.t. θ. Assuming the physics optimizer can account for varying function sensitivity in its updates ∆x, the resulting PG updates ∆θ 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 method are summarized in Fig. 2b.
While the L2 measure in Eq. 5 looks like a supervised objective in x 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, P−1(y∗ |x) = x− ( ∂P ∂x )T (P(x)− y∗). Then the gradient for Eq. 5
is ( ∂P ∂x ∂x ∂θ )T (P(x)− y∗), which is identical to the unsupervised training in Eq. 2 (the proof is given in appendix A.1).
Algorithm 1: Neural network training with physical gradients according to Eq. 5 for each training sample do
x0 ← NN(y∗ | θ) Solution inference with NN for k = 1, ..., n do
xk ← P−1(y∗ |xk−1) Accumulate iterative updates for PG end ∆θ ← ( ∂x ∂θ
)T · (x0 − xn)T Backpropagation through NN θ ← Step(θ,∆θ, η) Parameter update
end
Algorithm 1 details the corresponding neural network training procedure. There, the effective gradient ∆θ can also be written as ∂∂θ 1 2‖x0 − xn‖ 2 2 with xn taken as constant as in Eq. 5, which lends itself to straightforward integration into machine learning frameworks where a simple L2 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 a traditional supervised loss in x space (see appendix A.1 for a more detailed discussion).

heading: 3 RESULTS
text: We first demonstrate the differences between various neural network training schemes on a simple one-parameter example before moving on to physical systems described by partial differential equations. We specifically consider Poisson’s equation, the heat equation, and the NavierStokes equations. This selection covers ubiquitous physical processes with diffusion, transport, and strongly non-local effects, featuring both explicit and implicit solvers. In each experiment, we train identically-initialized networks using different gradient schemes for P to examine the relative differences in convergence behavior. For the weight update, we primarily employ the Adam optimizer (Kingma & Ba, 2015), a 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 are therefore representative of the performance on unseen data. A detailed description of our experiments along with additional visualizations and performance measurements can be found in appendix B.

heading: 3.1 SINGLE-PARAMETER OPTIMIZATION
text: First, we consider the task of finding solutions x∗ ∈ R to the inverse problem (Eq. 1) Pa(x) = sin(ax ), y
∗ = −1, given a ∈ [0.1, 10.1). This problem may seem simple at first glance but P shares some properties with chaotic physical systems, which makes this problem hard to optimize. We train a neural network with three fully-connected hidden layers to solve this task; the learning curves are shown in Fig. 3.
Supervised training schemes (Eq. 3) perform poorly on this task due to 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 a proximal minimum. However, ∂P∂x oscillates without bound when x→ 0, which causes overflow errors in the optimization. To avoid this, we clip the gradients ∂P∂x to [−1, 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) on this problem: Newton’s method and an analytic solver. As Fig. 3 shows, the PG methods vastly outperform both supervised and unsupervised first-order training, converging to an 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.

heading: 3.2 POISSON’S EQUATION
text: Poisson’s equation, P(x) = ∇−2x, plays an important role in electrostatics, Newtonian gravity, and fluid dynamics (Ames, 2014). It has the property that local changes in x can affect P(x) globally. Here we consider a two-dimensional system and train a U-net (Ronneberger et al., 2015) to solve inverse problems (Eq. 1) on pseudo-randomly generated y∗. Fig. 4c 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 L2 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 a full order of magnitude after 5k iterations.

heading: 3.3 HEAT EQUATION
text: Next, we consider a system with fundamentally non-invertible dynamics. The heat equation, ∂u∂t = ν · ∇2u, models heat flow in solids but also plays a part in many diffusive systems (Droniou, 2014). It gradually destroys information as the temperature equilibrium is approached (Grayson, 1987), causing ∇P to become near-singular. Inspired by heat conduction in microprocessors, we generate examples xGT by randomly scattering four to ten heat generating rectangular regions on a plane and
simulating the heat profile y∗ = P(xGT) as observed from outside a heat-conducting casing. We train a U-net (Ronneberger et al., 2015) to solve the corresponding inverse problem (Eq. 1). The learning curves are shown in Fig. 5c.
Since ∂P∂x is stable, individual inverse problems can be optimized with GD (see Fig. 5b). For the unsupervised network training (Eq. 2), we use Adam with η = 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 it hard for first-order methods to recover the solution.
For training with PGs (Eq. 5), we consider L-BFGS-B and the analytically derived inverse-physics solver P−1. L-BFGS-B outperforms GD when optimizing single inverse problems but is limited by the ill-conditioning (see Fig. 5b) of the problem. Instead, we use the available domain knowledge 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-physics solver, on the other hand, manage to improve the convergence speed substantially. They perform on par with L-BFGS-B training for the first ∼35 iterations. However, 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-physics solver requires only 138 iterations.
Additional examples and learning curves are shown in section appendix B.4.

heading: 3.4 NAVIER-STOKES EQUATIONS
text: 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: ∂v∂t = ν∇
2v − ∇ · ∇v − ∇p, ∇·v = 0,∇×p = 0, where p denotes pressure and ν the viscosity. At t = 0, a region of the fluid is randomly marked with a massless colorant m0 that passively moves with the fluid, ∂m∂t = −v · ∇m, and after a time t, the marker is observed again, mt. An example observation pair y∗ = {m0,mt} is shown in Fig. 6a. We target the inverse problem (Eq. 1) of finding an initial fluid velocity x ≡ v0 such that the fluid simulation P matches mt at time t. Since P is deterministic, x encodes the complete fluid flow from 0 to t. 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. 1. We train a U-net (Ronneberger et al., 2015) to solve these inverse problems; the learning curves are shown in Fig. 6c.
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 a standstill. This is caused by an overshoot in terms of vorticity, as visible in Fig. 6a right. While the resulting dynamics can roughly approximate the shape of the observed mt, 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 to the highly non-linear physics. A similar behavior can also be seen when optimizing single inverse problems with GD (Fig. 6b). 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 (Eq. 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 mt, then estimates x 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 possible because the inverse-physics solver and associated PG does not suffer from the strongly-varying
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 as 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 1 lists measured training and inference times. With both iterative solver and network, we used a 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.

heading: 4 CONCLUSION
text: 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 to 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 on a wide variety of inverse problems involving partial differential equations while being on par with regular gradient evaluation in terms of computational cost.
We anticipate our method to be applicable to a wide range of differentiable processes. Beyond physics, interesting avenues lie in applying PGs to differentiable rendering or training invertible neural networks.

heading: REPRODUCIBILITY STATEMENT
text: We will publish the full source code and trained network models upon acceptance. They will be sufficient to 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.

heading: A METHOD

heading: A.1 EMBEDDING THE PHYSICAL GRADIENT INTO THE GRADIENT DESCENT PIPELINE
text: The key for embedding inverse-physics optimizers in the machine learning pipeline is the proxy L2 loss in Eq. 5 which passes the gradients from the physics optimizer to the GD-based neural network optimizer.
This l2 loss may look like a supervised loss in x space but due to the dependency of the labels P−1 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 a regular first-order optimization. The proxy L2 loss simply connects the computational graphs used in backpropagation.
Physical gradients based on gradient descent are equal to unsupervised training With the shorthands xi = NN(y∗i | θ), yi = P(xi) and ∆yi = yi − y∗i , we can write the objective functions of unsupervised training and physical gradient training as
U(θ) = N∑ i=1 1 2 ||∆yi||22 (unsupervised) (6)
M(θ) = N∑ i=1 1 2 ||xi − P−1(y∗i |xi)||22 (physical gradients) (7)
Theorem 1. Minimizing U(θ) is identical to minimizing M(θ) using any gradient-descent-based optimizer if P−1 is implemented as a single gradient descent step with η = 1.
Proof. The gradient of the unsupervised objective function is
∂U ∂θ = N∑ i=1 ∆yi · ∂y ∂x ∂x ∂θ
For the mixed training, we recall that the objective function of a single inverse problem is L(x) = 1 2 ||P(x)− y ∗||22. The corresponding gradient is ∂L∂x = ∆y · ∂y ∂x . Inserting a gradient descent step for P−1, the gradient of the mixed objective function becomes
∂M
∂θ = N∑ i=1 ( xi − P−1(y∗i |xi) ) ∂x ∂θ
= N∑ i=1 ( xi − ( xi − η · ∂L ∂x )) ∂x ∂θ
= N∑ i=1 η ·∆yi · ∂y ∂x ∂x ∂θ
which is equal to ∂U∂θ for η = 1.
This generic L2 formulation is even applicable in settings where different parts of an end-to-end pipeline are computed by different software frameworks or hardware accelerators.

heading: A.2 CONVERGENCE WITH PHYSICAL GRADIENTS
text: We consider the convergence of a coupled optimization of Li(xi) with xi = fθ(yi) involving gradient descent optimization of a general function approximator fθ and some other optimization scheme P for minimizing objective functions Li(xi). Here, i denotes the index of a specific example y from a data set of size N . In the context of learning to solve inverse problems, P is the inverse-problem solver on which physical gradient P (x) − x is based, and fθ represents the neural network. In the main text, we considered objectives of the form Li = 12 ||P(xi)−y ∗ i ||22, but now we consider a more general form.

heading: Assumptions
text: 1. Let L = {Li : Rd → R | i = 1, ..., N} be a finite set of functions and let x∗i ∈ Rd be any global minima of Li.
2. Let P = {Pi : Rd → Rd | i = 1, ..., N} be a set of update functions for Li such that ∃τ > 0 : ∀x ∈ Rd : Li(x)−Li(Pi(x)) ≥ τ (L(x)− L(x∗)) and ∃K > 0 : ∀x ∈ Rd : ||Pi(x)− x|| ≤ K(L(x)− L(x∗)).
3. Let fθ : Rm → Rd be differentiable w.r.t. θ with the property that ∃η > 0 : ∀i ∈ 1, ..., N ∀xi ∈ Rd ∀ > 0 ∃n ∈ N : ||fθxin (yi) − xi||2 ≤ where θ x n is the sequence of
gradient descent steps with θxn+1 = θ x n − η ( ∂fθ ∂θ )T (fθn − x), η > 0.
Here we make minimal assumptions about L; not even continuity is required as long as we have access to some form of a physical gradient that can optimize it in x. We denote k repeated applications of P by P k(x). The assumption in (3) states that the function fθ is flexible enough to fit our problem. It must be able to converge to every point x for all examples using gradient descent. It 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 a 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 k ∈ N, there exists an update strategy θn+1 = Uk(θn) based on a single evaluation of P k for which L(fθn(yi)) converges to a minimum x ∗ i or minimum region of Li ∀i =

heading: 1, ..., N .
text: Proof. Dropping the example index i, we denote x̃n ≡ P k(xn) and ∆L∗ = L(xn)−L(x∗). Let I2 denote the open set of all x for which L(x)−L(xn) > τ2∆L ∗ Definition (2) provides that x̃n ∈ I2.
Since I2 is open, ∃ > 0 : ∀x ∈ B (x̃n) : L(xn) − L(x) > τ2∆L ∗ where B (x) denotes the open set containing all x′ ∈ Rd for which ||x′ − x||2 < , i.e. there exists a small ball around x̃n which is fully contained in I2 (see sketch in Fig. 7).
Using the convergence assumption for fθ, definition (3), we can find a finite n ∈ N for which fθn ∈ B (x̃n) and therefore L(xn)− L(fθn) > τ2∆L
∗. We can thus use the following strategy Uk for minimizing L(fθ): First, compute x̃n = P k(xn). Then perform gradient descent steps in θ with the effective objective function 12 ||fθ− x̃n|| 2 2 until L(xn)−L(fθ) ≥ τ2∆L
∗. With this strategy, each application ofUk reduces the loss to ∆L∗n+1 ≤ (1− τ2 )∆L ∗ n so any arbitrary value ofL > L(x
∗) can be reached within a finite number of steps. Definition (2) also ensures that ||P (x)− x|| → 0 as the optimization progresses which guarantees that the optimization converges to a minimum region.
While this guarantees convergence, it requires potentially many gradient descent steps in θ 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 θ, and θ is far away from a solution. However, in many cases, we want to update the physical gradient after each update to θ. Without additional assumptions about P and fθ, 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 12 ||fθ − x̃n|| 2 2, the next value of fθ must lie closer to x̃n than xn. Assuming that, on average, the gradient descent steps do not prefer a certain direction, xn+1 is more likely to lie in I = {x ∈ Rd : L(x) < L(xn)} 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 x̃n, L is more likely to decrease than increase.
While this shows that the loss should decrease on average, it does not guarantee convergence. We now look at two specific formulations of P for which convergence to the correct solution is guaranteed. 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. If ∀x ∈ Rd ∃λ ∈ (0, 1] : P (x) = x + λ(x∗ − x), then the sequence fθn with
θn+1 = θn − η ( ∂f ∂θ )T (x− P (x)) converges to x∗. Proof. Rewriting x − P (x) = ∂∂x ( λ 2 ||x− x ∗||22 ) yields the update θn+1 − θn = −η ( ∂L̃ ∂x ∂f ∂θ
)T where L̃ = λ2 ||x− x
∗||22. This describes pure gradient descent towards x∗ with the gradients scaled by λ. Since λ ∈ (0, 1], the convergence proof of gradient descent applies.
Theorem 4. If ∀x ∈ Rd ∃λ ∈ (0, 1] : P (x) = x − λ ( ∂L ∂x )T , then the sequence fθn with
θn+1 = θn − η ( ∂f ∂θ )T (x− P (x)) converges to minima of L.
Proof. This is equivalent to gradient descent in L(θ). Rewriting the update yields θn+1 − θn = −η ( ∂f ∂θ )T ( x− ( x− ( λ∂L∂x )T)) = −ηλ ( ∂L ∂x ∂f ∂θ )T which is the gradient descent update scaled by λ.

heading: B EXPERIMENTS
text: Here, we give a more detailed description of our experiments including setup and analysis. The implementation of our experiments is based on the ΦFlow (PhiFlow) framework (Holl et al., 2020) and uses TensorFlow (Abadi et al., 2016) and PyTorch (Paszke et al., 2017) for automatic differentiation. Our code is open source and will be made available upon acceptance.

heading: B.1 WAVE PACKET FIT
text: This experiment is an instance of a 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 as inverse solvers over classical optimization algorithms.
Data generation. We simulate an observed time series y∗ from a random ground truth position and amplitude x∗ = {ξ0, A}. Each time series contains 256 entries and is made up of the wave packet and noise. For the wave packet, we sample ξ0 ∈ [25.6, 230.4) and A ∈ [1, 4) from uniform distributions. The wave packet has the functional form
y(ξ) = A · sin(f · ξ) · exp ( −1
2
(ξ − ξ0)2
σ2 ) where we set f = 45 and σ = 20 constant for all data. For the noise, we superimpose random values sampled from the normal distribution N (0, 25 ) 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 a max pooling operation and two 1D 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 a third fully-connected layer produces the predicted ξ0 and A. ξ0 is passed through a Sigmoid activation function and normalized to values between 25.6 and 230.4.
Training and fitting. We fit the data using L-BFGS-B with an initial guess of ξ0 = 128 andA = 2 and the network output is offset by the same amount. Both network and L-BFGS-B minimize the squared loss ||y(ξ0, A)−y∗||22 and the resulting performance curves are shown in Fig. 8. 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 does not find it and instead fits the noise near the center.
The neural network is trained using Adam with learning rate 0.001 and batch size of 100. Despite the simpler, Hessian-free, optimization updates, the network quickly learns to localize all wave packets,
outperforming L-BFGS-B after 30 to 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.

heading: B.2 SINGLE-PARAMETER OPTIMIZATION
text: The optimization consists of finding x∗ = argminx sin( a x ). Fig. 9 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 x0. We solve this analytically by first determining the index of the closest minimum, n = a2π·x0 − 3 4 which we round to the closest integer. The furthest-out minima (n = 0 and n = −1) require special handling. Then the position of the minimum corresponding to n is x∗ = a3π/2+2π·n .
Newton’s method The analytic form of a Newton update is
P(x)− x = − x 2 cos(a/x)
2x cos(a/x)− a 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 ∂P∂x . Since Newton’s method relies on an inversion of the Hessian, the produced physical gradients diverge where ∂
2P ∂x2 = 0. To avoid overflow, we clip the gradients ∂P ∂x to [−1, 1] before they are
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 as 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 a 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
P(x)− x = − h ∇P = −x2 sin(a/x) + 1 a cos(a/x)
where h = sin ( a x ) + 1 measures the height above y∗. It approaches ±∞ when sin(ax ) is maximal. This inversion requires domain knowledge since h is based on the fact that minP = −1. Like with Newton’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 as 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 a learning rate of 0.002 and minibatches containing 1000 uniformly sampled a ∈ [0.1, 10.1]. We use ΦFlow together with PyTorch to train the neural network in this experiment.
Learning curves for different network initializations are shown in Fig. 10 and the recorded computation times are shown in Fig. 14. 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, it sometimes stops converging to an accuracy larger than 10−4. This is most likely caused by numerical effects of the inverse solver. Training with physical gradients based on Newton’s method or inverse gradients yields a convergence that strongly depends on the network initialization. This variant converges only for certain initializations. However, upon convergence, it typically reaches an accuracy similar to the inverse physics variant.

heading: B.3 POISSON’S EQUATION
text: We consider Poisson’s equation, ∇2y = x where x is the initial state and y is the output of the simulator. We set up a two-dimensional simulation with 80 by 60 cubic cells. Our simulator computes
y = P(x) = ∇−2x implicitly via the conjugate gradient method. The inverse problem consists of finding an initial value x∗ for a given target y∗ such that ∇2y∗ = x∗. We formulate this problem as minimizing L(x) = 12 ||P(x) − y ∗||22 = 12 ||∇ −2(x − x∗)||22. We now investigate the computed updates ∆x of various optimization methods for this problem. Gradient descent Gradient descent prescribes the update ∆x = −η · ( ∂L ∂x )T = −η ·∇−2 (y − y∗) which requires an 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 η is very small. We found that GD requires η ≤ 2 · 10−5, while the momentum in Adam allows for larger η. 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, yielding x = ∇2y. Correspondingly, we formulate the update step as ∆x = −η · ∂x∂y · (y − y ∗) = −η · ∇2 (y − y∗) which directly points to x∗ for η = 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 is 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 Y . Both generation of x∗ and y∗ is performed on the fly, resulting in a data set of effectively infinite size, |Y| = ∞. This has the advantage that learning curves are representative of both test performance as well as training performance. The top of Fig. 11 shows some examples generated this way. We train a U-net (Ronneberger et al., 2015) 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 η = 10−3 with Adam and η = 10−9 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 η = 2 · 10−5 was required. We use a typical value of 0.9 for the momentum of SGD and Adam. For the training using Adam with physical gradients, we compute ∆x as described above and keep η = 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 are shown at the bottom of Fig. 11, while Fig. 14 shows the computation time per iteration.

heading: B.4 HEAT EQUATION
text: We consider a two-dimensional system governed by the heat equation ∂u∂t = ν · ∇ 2u. Given an initial state x = u0 at t0, the simulator computes the state at a later time t∗ via y = u∗ = P(x). Exactly inverting this system is only possible for t · ν = 0 and becomes increasingly unstable for larger t ·ν because initially distinct heat levels even out over time, drowning the original information in noise. Hence the Jacobian of the physics ∂y∂x is near-singular. In our experiment we set t ·ν = 8 on a domain consisting of 64x64 cells of unit length. This level of diffusion is challenging, and diffuses most details while leaving the large-scale structure intact.
We apply periodic boundary conditions and compute the result in frequency space where the physics can be computed analytically as ŷ = x̂ · e−k2(t∗−t0) where ŷk ≡ F(y)k denotes the k-th element of the Fourier-transformed vector y. Here, high frequencies are dampened exponentially. The inverse
y * Adam y Adam+PG y Adam x Adam+PG x x *
problem can thus be written as minimizing L(x) = ||P(x)−y∗||22 = ||F−1 ( F(x) · e−k2(t∗−t0) ) − y∗||22
Gradient descent Using the analytic formulation, we can compute the gradient descent update as ∆x = −η · F−1 ( e−k 2(t∗−t0)F(y − y∗) ) .
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 converge slowly on this task after fitting the coarse structure and have severe problems in recovering high-frequency details. This is not because the information is fundamentally missing but because GD cannot adequately process high-frequency details.
Stable physical gradients The frequency formulation of the heat equation can be inverted analytically, yielding x̂k = ŷk · ek 2(t∗−t0). This allows us to define the update
∆x = −η · F−1 ( ek 2(t∗−t0)F(y − y∗) ) .
Here, high frequencies are multiplied by exponentially large factors, resulting in numerical instabilities. When applying this formula directly to the gradients, it can lead to large oscillations in ∆x. This is the opposite behavior compared to Poisson’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 a certain amount of noise n, with the remainder constituting the signal s = y − n. For the noise, we assume a normal distribution n ∼ N (0, · y) with > 0 and for the signal, we assume that it arises from reasonable values of x so that y ∼ N (0, δ · e−k2) with δ > 0. With this, we can estimate the probability of an observed value arising from the signal using Bayes’ theorem p(s|v) = p(v|s)·p(s)p(v|s)·p(s)+p(v|n)·p(n) where we assume the priors p(s) = p(n) = 1 2 . Based on this probability, we dampen the amplification of the inverse physics which yields a stable inverse. Gradients computed in this way hold as much high-frequency information as can be extracted given the noise that is present. This leads to a much faster convergence and more precise solution than any generic optimization method.
Neural network training For training, we generate x∗ by randomly placing between 4 and 10 hot rectangles of random size and shape in the domain and computing y = P(x∗). Both operations are executed on the fly so that |Y| = ∞ and the learning curves are representative of both training and test performance. For the neural network, we use the same U-net architecture as in the previous experiment. We train with a batch size of 128 and a constant learning rate of η = 10−3. The network updates are computed with TensorFlow’s or PyTorch’s automatic differentiation. Fig. 12 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.

heading: B.5 NAVIER-STOKES EQUATIONS
text: 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 a direct numerical solver. We adopt the marker-in-cell (MAC) method (Harlow & Welch, 1965; Harlow, 1972) which guarantees stable simulations 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 v0 is specified at cell centers and resampled to a staggered grid for the simulation. Our simulation employs a second-order advection scheme (Selle et al., 2008) to transport both the marker and the velocity vectors. This step introduces significant amount of numerical diffusion which can clearly be seen in the final marker distributions. Hence, we do not numerically solve
for adding additional viscosity. Incompressibility is achieved via Helmholz decomposition of the velocity field using a conjugate gradient solve.
Neither pressure projection nor advection are energy-conserving operations. While specialized energy-conserving simulation schemes for fluids exist (Gluhovsky & Tong, 1999; Mullen et al., 2009), we instead enforce energy conservation by normalizing the velocity field at each time step to the total energy of the previous time step. Here, the energy is computed as E = ∫ R2 dx v(x)
2 since we assume constant fluid density.
Data generation The data set consists of marker pairs {m0,mt} which are randomly generated on-the-fly. For each example, a center position for m0 is chosen on a grid of 64x64 cells. m0 is then generated from discretized noise fluctuations to fill half the domain size in each dimension. The number of marked cells is random.
Next, a ground truth initial velocity v0 is generated from three components. First, a uniform velocity field moves the marker towards the center of the domain to avoid boundary collisions. Second, a 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 are added additionally perturb the flow fields. These are 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 m0 and v0, a ground truth simulation is run for t = 2 with ∆t = 0.25. The resulting marker density is then used as the target for the optimization. This ensures that there exists a solution for each example.
Computation of physical gradients To compute the physical gradients for this example, we construct an explicit formulation v̂0 = P−1(m0,mt |x0) that produces an estimate for v0 given an initial guess x0 by locally inverting the physics. From this information, it fits the coarse velocity, i.e. the uniform velocity and the vortex present in the data. This use of domain knowledge, i.e., enforcing the translation and rotation components of the velocity field as a prior, is what allows it to produce a much better estimate of v0 than the regular gradient. More formally, it assumes that the solution lies on a manifold that is much more low-dimensional than v0. 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 limitation. 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 mrevt = mt and velocity v rev t = v
t from the forward simulation. The reverse simulation then computesmrev and vrev for all time steps by performing simulation steps with ∆t = −0.25. Then, the update to the vortex strength is computed from the differences mrev −m at each time step and an estimate of the vortex location at these time steps.
Neural network training We train a U-net (Ronneberger et al., 2015) similar to the previous experiments but with 5 resolution levels. The network contains a total of 49,570 trainable parameters. The network is given the observed markers m0 and mt, resulting in an input consisting of two feature 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∗)| · w where F denotes the two-dimensional Fourier transform and w is a weighting vector that factors high frequencies exponentially less than low frequencies.
We train the network using Adam with a learning rate of 0.005 and mini-batches containing 64 examples each, using PyTorch’s 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 high-frequency motions. Example trajectories and reconstructions are shown in Fig. 13 and performance measurements are shown in Fig. 14.

Review summaries:
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. 

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).

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.   The functionality of the proposed method is backed by theoretical results and four numerical experiments.