# Nature’s Cost Function: Simulating Physics by Minimizing the Action

Tim Strang
ML Collective
Holliston, MA, 01746
<EMAIL_ADDRESS>
&Isabella Caruso
MIT Materials Science
Cambridge, MA, 02142
<EMAIL_ADDRESS>
&Sam Greydanus
ML Collective
Corvallis, Oregon, 97330
<EMAIL_ADDRESS>

###### Abstract

In physics, there is a scalar function called the action which behaves like a
cost function. When minimized, it yields the “path of least action” which
represents the path a physical system will take through space and time. This
function is crucial in theoretical physics and is usually minimized
analytically to obtain equations of motion for various problems. In this
paper, we propose a different approach: instead of minimizing the action
analytically, we discretize it and then minimize it directly with gradient
descent. We use this approach to obtain dynamics for six different physical
systems and show that they are nearly identical to ground-truth dynamics. We
discuss failure modes such as the unconstrained energy effect and show how to
address them. Finally, we use the discretized action to construct a simple but
novel quantum simulation. Code: github.com/greydanus/ncf

Figure 1: Finding a path of least action with gradient descent. Left: We
compare the normal approach of ODE integration to our approach of action
minimization. The action is the sum, over every point in the path, of kinetic
energy $T$ minus potential energy $V$. We compute the gradients of this
quantity with respect to the path coordinates and then deform the initial path
(yellow) into the path of least action (green). This path resolves to a
parabola, matching the path obtained via ODE integration. Right: We plot the
path’s action $S$, kinetic energy $T$, and potential energy $V$ over the
course of optimization. All three quantities asymptote at the $S$, $T$, and
$V$ values of the ODE trajectory.

##  1 Introduction

Over the past two centuries, many works have studied the action in a variety
of physical and mathematical settings (Goldstine, 2012; Ferguson, 2004). In
all these studies, the authors treated the action analytically. They only
solved for paths of least action indirectly by first invoking the Euler-
Lagrange equation, then solving for a system of differential equations, and
finally integrating those equations over time with an ordinary differential
equation (ODE) solver.

In this paper, we take a different approach:

> Instead of solving for paths of stationary action analytically, we
> discretize the action function and then minimize it directly with gradient
> descent.

We use this approach to obtain dynamics for six different physical systems: a
body in free fall, a pendulum, a double pendulum, the three body problem, a
gas with a Lennard-Jones potential, and the planetary ephemerides of the inner
solar system. We document a failure mode called the unconstrained energy
effect and discuss how to mitigate it. Finally, we note that the action
controls the resonance modes of paths in quantum mechanics and use our
discretized version to construct a simple quantum simulation.

##  2 Background

The Lagrangian method. The Lagrangian method begins by considering all the
paths a physical system could take from an initial state $\bf x(t_{0})$ to a
final state $\bf x(t_{1})$. Then it provides a simple rule for selecting the
path $\hat{\bf x}$ that nature will actually take: the action $S$, defined in
Eqn. 1, must have a stationary value over this path. Here $T$ and $V$ are the
kinetic and potential energy functions for the system at any given time $t$ in
$[t_{0},t_{1}]$.

$\displaystyle S$ $\displaystyle:=\int_{t_{0}}^{t_{1}}L({\bf x},~{}\dot{\bf
x},~{}t)~{}dt\quad\textrm{where}\quad L=T-V$ (1) $\displaystyle\quad\hat{\bf
x}$ $\displaystyle~{}~{}\textrm{has the
property}\quad\frac{d}{dt}\left(\frac{\partial
L}{\partial\dot{\hat{x}}(t)}\right)=\frac{\partial
L}{\partial\hat{x}(t)}\quad\textrm{for}\quad t\in[t_{0},t_{1}]$ (2)

Finding $\hat{\bf x}$ with Euler-Lagrange. When $S$ is stationary, we can show
that the Euler-Lagrange equation (Eqn. 2) holds true over the interval
$[t_{0},t_{1}]$ (Morin, 2008). This observation is valuable because it allows
us to solve for $\hat{\bf x}$: first we apply the Euler-Lagrange equation to
the Lagrangian $L$ and derive a system of partial differential equations. Then
we integrate those equations to obtain $\hat{\bf x}$. Importantly, this
approach works for all problems spanning classical mechanics, electrodynamics,
thermodynamics, and relativity. It provides a coherent theoretical framework
for studying classical physics as a whole.

Finding $\hat{\bf x}$ with action minimization (this work). A more direct
approach to finding $\hat{\bf x}$ begins with the insight that paths of
stationary action are almost always also paths of least action Morin (2008).
Thus, without much loss of generality, we can exchange the Euler-Lagrange
equation for the simple minimization objective shown in the third part of Eqn.
3. Meanwhile, as shown in the first part of Eqn. 3, we can redefine $S$ as a
discrete sum over $N$ evenly-spaced time slices:

$\displaystyle S:=\sum_{i=0}^{N}L({\bf x},~{}\dot{{\bf x}},~{}t_{i})\Delta
t\quad\textrm{where}\quad\dot{{\bf x}}(t_{i}):=\frac{{\bf x}(t_{i+1})-{\bf
x}(t_{i})}{\Delta t}\quad\textrm{and}\quad\hat{\bf
x}:=\operatorname*{arg\,min}_{\bf x}S(\bf x)$ (3)

One problem remains: having discretized $\hat{\bf x}$ we can no longer take
its derivative to obtain an exact value for $\dot{\bf x}(t_{i})$. Instead, we
must use the finite-differences approximation shown in the second part of
Equation 3. Of course, this approximation will not be possible for the very
last $\dot{\bf x}$ in the sum because $\dot{\bf x}_{N+1}$ does not exist. For
this value we will assume that, for large $N$, the change in velocity over the
interval $\Delta t$ is small and thus let $\dot{\bf x}_{N}=\dot{\bf x}_{N-1}$.
Having made this last approximation, we can now compute the gradient
$\frac{\partial S}{\partial{\bf x}}$ numerically and use it to minimize $S$.
This can be done with PyTorch (Paszke et al., 2019) or any other package that
supports automatic differentiation.

Related work: minimizing the Onsager-Machlup action. One line of related work
involves minimizing the action of the Onsager-Machlup (OM) function (Onsager &
Machlup, 1953; Faccioli et al., 2006). It is used to describe the time
dynamics of the probability density of a stochastic process. The action of the
OM function has a different physical interpretation from the action of a
Lagrangian but shares certain mathematical properties. Adib (2008) discretize
diffusion trajectories and calculate their gradients (Eqn. 22) but do not
minimize the OM action. Conversely, Lee et al. (2017); Zuckerman & Woolf
(2000) minimize the OM action – to analyze protein folding dynamics – but use
gradient-free optimization methods: Metropolis MCMC and Conformational Space
Annealing respectively (Metropolis et al., 1953; Lee et al., 1997). Like our
work, these works fix the initial and final states and (sometimes) minimize
the action over the interval. Unlike our work, they use the OM function and do
not perform minimization with automatic differentiation or gradient descent.
We discuss additional related work in the context of gradient-based
optimization of physics problems in Appendix A.1.

Figure 2: Visualizing our technique applied to a variety of different physics
simulations. First row: We show diagrams of the setups and variables of the
six systems under investigation. Second row: we visualize the ODE paths
(purple), initial semi-random paths (yellow), and final paths of minimized
action (green). These plots give an intuition for how the paths deform during
optimization. Third row: We visualize the optimization dynamics in terms of
action, kinetic energy, and potential energy.  Figure 3: The unconstrained
energy effect in the pendulum task. Optimization converges to the baseline
path initially. After additional optimizer steps, it converges to a different
but also valid solution where the pendulum swings up, hangs vertically, and
then falls into the final state. This occurs because our current methods do
not constrain the total energy $T+V$ in the way the ODE initial state $({\bf
x}(t_{0}),\dot{{\bf x}}(t_{0}))$ does.

##  3 Six Experiments

Minimizing the action in the way we have described has not been studied in
detail. We have discussed a few related works which do this in the context of
the OM function but their scope and details diverge from this work and each
other. Thus in our experiments we prioritized simplicity. Unless otherwise
specified, we set all constants such as mass and gravity to one. When
selecting physical systems, we began with two toy problems (for debugging): a
free body and a pendulum. Then we investigated four more complex systems: a
double pendulum, the three body problem, a simple gas, and a real ephemeris
dataset of planetary motion. These systems presented an interesting challenge
because they were all nonlinear, chaotic, and high-dimensional111The state of
the simple gas, for example, has a hundred degrees of freedom.. In each case,
we compared our results to a baseline path obtained with a simple ODE solver
using Euler integration. Appendix A.2 gives more details, including the
Lagrangians and equations of motion.

The unconstrained energy effect. Early in our experiments we encountered the
unconstrained energy effect. This happens when the optimizer converges on a
valid physical path with a different total energy from the baseline. Figure 3
shows an example. The reason this happens is that, although we fix the initial
and final states, we do not constrain the path’s total energy $T+V$. Even
though paths like the one in Figure 3 are not necessarily invalid, they make
it difficult for us to recover baseline ODE paths. For this reason, we use the
baseline ODE paths to initialize our paths, perturb them with Gaussian noise,
and then use early stopping to select for paths which are similar (often,
identical) to the ODE baselines. This approach matches the mathematical ansatz
of the “calculus of variations” where one studies perturbed paths in the
vicinity of the true path of least action. We note that there are other ways
to mitigate this effect which don’t require an ODE-generated initial path. We
discuss them in Appendix A.2, as they are beyond the main scope of this work.

Results. On all six physical systems we obtained paths of least action which
were nearly identical to the baseline paths of the ODE solver. Figure 2 shows
optimization dynamics and qualitative results while Table 2 shows quantitative
results. These results suggest that action minimization can generate
physically-valid dynamics even for chaotic and strongly-coupled systems like
the double pendulum and three body problem. One interesting pattern we noticed
was that optimization dynamics were dominated by the kinetic energy term $T$
(third row of Figure 2). This occurs because $S$ tends to be more sensitive to
$T$ (which grows as $\dot{{\bf x}}^{2}$) than $V$. Future methods should focus
on stabilizing $T$.

##  4 The Quantum Limit

Figure 4: Visualizing wave packet dynamics at different scales ($\hbar$
$\propto$ Scale-1). As $\hbar\rightarrow 0$, the waves retain coherence longer
and interfere with the right wall at higher frequencies, approaching the
behavior of a classical particle at Scale = 2. Coherence grows as oscillations
along the off-diagonals of the phase-action matrix, $K$, increase in
frequency.

The equation for the action, which Lagrange introduced in the context of
classical mechanics, is one of the few equations which does not need to be
redefined in the quantum context. Here, systems are described with a wave
function and allowed to take on a superposition (linear combination) of paths
at once. The best-known example of this phenomena is the famous double-slit
experiment, where a single electron fired at two separate slits travels
through both and interferes with itself on the other side. This interference
occurs because the phase evolves differently with each path and thus is out of
sync when the paths recombine. The exponential portion of Eqn. 4 describes how
a system evolves as it travels along a given path: it oscillates at a rate
proportional to the path’s action. Thus the phase of the wave function, after
traveling along a given path for some interval of time, is also proportional
to the action. The rest of Eqn. 4 describes how to recombine the contributions
of all paths after some interval of time. First, it rotates the wave function
by the phase change associated with each path and then it sums (integrates)
these states to obtain the full state at the end of the time interval. In the
double slit example, this summation would recombine the waves from the two
slits into the full wave function $\Psi(\bf x,t)$ sampled at the sensor.

$\displaystyle\Psi(x,t)=\int\limits_{\bf x(0)=x}\Psi(\bf
x(t),0)*exp\\{\frac{i}{\hbar}\int L[\bf
x(\tau),\;\bf\dot{x}(\tau),\tau]d\tau\\}\mathcal{\;}\mathcal{D}\bf x$ (4)

Discretizing the path integral. Feynman et al. (1965) and many others have
treated Eqn. 4 analytically, deriving from it both Schrödinger’s equation and
the classical principle of least action. However, as was the case with the
classical action, few works have discretized the paths and used the action to
construct a numerical simulation. Once again, that is our goal. It can be done
by dividing the spatial and time dimensions of Eqn. 4 into small, uniform
intervals. Then, for each time interval, we can restrict our group of paths to
only the $N^{2}$ linear paths that connect the $N$ discrete spatial locations
to one another. We can calculate their actions and, with the help of Eqn. 4,
their phase changes. Then we can organize these phase changes in an $N\times
N$ phase-action matrix, $K$, where rows correspond to locations of path start
points, columns correspond to locations of path endpoints, and the matrix
elements themselves are the phase changes of the paths. This matrix, plotted
in Figure 4, provides an intuitive way to visualize how nearby paths interfere
with one another at the quantum-classical limit.

A simple quantum simulation. Now, to simulate dynamics over one discrete time
slice, we can simply multiply the wave function (a discretized vector of
complex amplitudes at the $N$ spatial locations) by $K$. Then we can evolve it
arbitrarily far into the future by multiplying by $K$ repeatedly (see Figure
4). The quantum simulation we recover is fundamentally different from the more
popular Hamiltonian-based simulations (Hartree & Hartree, 1935; Richings et
al., 2015) and is remarkable for its simplicity. It mirrors the philosophical
approach we took in the six classical experiments in that, in both cases, we
express computation in terms of the action as directly as possible.

##  5 Conclusion

Basic physical principles and the elegant reasoning behind them are often
obscured in the midst of numerical approximations and domain-specific
notation. In this work we have sought to prevent this from happening by
expressing computation in terms of the action as directly as possible.

#### Acknowledgments

All authors thank Jason Yosinski, Rosanne Liu, and other members of the ML
Collective research group for helpful discussion, ideas, and feedback on
experiments.

TS thanks the Central Rock Gym for an ambient environment encouraging of
intellectual reflection. IC thanks the MIT Materials Science Department for
unwittingly allowing her to indulge in this side project while waiting for
feedback on her thesis application. SG thanks TS and IC for joining him on
this research journey which was neither the shortest path nor the one of least
action.

## References

  * Adib (2008) Artur B Adib.  Stochastic actions for diffusive dynamics: Reweighting, sampling, and minimization.  _The Journal of Physical Chemistry B_ , 112(19):5910–5916, 2008.
  * Alder & Wainwright (1959) Berni J Alder and Thomas Everett Wainwright.  Studies in molecular dynamics. i. general method.  _The Journal of Chemical Physics_ , 31(2):459–466, 1959.
  * Bar-Sinai et al. (2019) Yohai Bar-Sinai, Stephan Hoyer, Jason Hickey, and Michael P Brenner.  Learning data-driven discretizations for partial differential equations.  _Proceedings of the National Academy of Sciences_ , 116(31):15344–15349, 2019.
  * Chen et al. (2018) Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud.  Neural ordinary differential equations.  _Advances in neural information processing systems_ , 31, 2018.
  * Cranmer et al. (2020) Miles Cranmer, Sam Greydanus, Stephan Hoyer, Peter Battaglia, David Spergel, and Shirley Ho.  Lagrangian neural networks.  _arXiv preprint arXiv:2003.04630_ , 2020.
  * Faccioli et al. (2006) P Faccioli, M Sega, F Pederiva, and H Orland.  Dominant pathways in protein folding.  _Physical review letters_ , 97(10):108101, 2006\.
  * Ferguson (2004) James Ferguson.  A brief survey of the history of the calculus of variations and its applications.  _arXiv preprint math/0402357_ , 2004.
  * Feynman et al. (1965) Richard P Feynman, Robert B Leighton, and Matthew Sands.  The feynman lectures on physics; vol. ii.  _American Journal of Physics_ , 33(9):750–752, 1965.
  * Freeman et al. (2021) C Daniel Freeman, Erik Frey, Anton Raichuk, Sertan Girgin, Igor Mordatch, and Olivier Bachem.  Brax–a differentiable physics engine for large scale rigid body simulation.  _arXiv preprint arXiv:2106.13281_ , 2021.
  * Goldstine (2012) Herman Heine Goldstine.  _A History of the Calculus of Variations from the 17th through the 19th Century_ , volume 5.  Springer Science & Business Media, 2012.
  * Hartree & Hartree (1935) Douglas Rayner Hartree and William Hartree.  Self-consistent field, with exchange, for beryllium.  _Proceedings of the Royal Society of London. Series A-Mathematical and Physical Sciences_ , 150(869):9–33, 1935\.
  * Hollingsworth & Dror (2018) Scott A Hollingsworth and Ron O Dror.  Molecular dynamics simulation for all.  _Neuron_ , 99(6):1129–1143, 2018.
  * Kingma & Ba (2014) Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Lee et al. (1997) Jooyoung Lee, Harold A Scheraga, and Shelly Rackovsky.  New optimization method for conformational energy calculations on polypeptides: conformational space annealing.  _Journal of computational chemistry_ , 18(9):1222–1232, 1997.
  * Lee et al. (2017) Juyong Lee, In-Ho Lee, InSuk Joung, Jooyoung Lee, and Bernard R Brooks.  Finding multiple reaction pathways via global optimization of action.  _Nature Communications_ , 8(1):15443, 2017.
  * Lutter et al. (2019) Michael Lutter, Christian Ritter, and Jan Peters.  Deep lagrangian networks: Using physics as model prior for deep learning.  _arXiv preprint arXiv:1907.04490_ , 2019.
  * Metropolis et al. (1953) Nicholas Metropolis, Arianna W Rosenbluth, Marshall N Rosenbluth, Augusta H Teller, and Edward Teller.  Equation of state calculations by fast computing machines.  _The journal of chemical physics_ , 21(6):1087–1092, 1953.
  * Morin (2008) David Morin.  _Introduction to classical mechanics: with problems and solutions_.  Cambridge University Press, 2008.
  * Onsager & Machlup (1953) Lars Onsager and Stefan Machlup.  Fluctuations and irreversible processes.  _Physical Review_ , 91(6):1505, 1953.
  * Paszke et al. (2019) Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al.  Pytorch: An imperative style, high-performance deep learning library.  _Advances in neural information processing systems_ , 32, 2019.
  * Richings et al. (2015) Gareth W Richings, Iakov Polyak, Kaite E Spinlove, Graham A Worth, Irene Burghardt, and Benjamin Lasorne.  Quantum dynamics simulations using gaussian wavepackets: the vmcg method.  _International Reviews in Physical Chemistry_ , 34(2):269–308, 2015.
  * Schoenholz & Cubuk (2020) Samuel Schoenholz and Ekin Dogus Cubuk.  Jax md: a framework for differentiable physics.  _Advances in Neural Information Processing Systems_ , 33:11428–11441, 2020.
  * Schroeder (2015) Daniel V Schroeder.  Interactive molecular dynamics.  _American Journal of Physics_ , 83(3):210–218, 2015.
  * Sweet et al. (2018) Chelsea Sweet, Oyewumi Akinfenwa, and Jonathan J Foley IV.  Facilitating students’ interaction with real gas properties using a discovery-based approach and molecular dynamics simulations.  _Journal of Chemical Education_ , 95(3):384–392, 2018.
  * Zhang et al. (2013) BL Zhang, J Wang, M Li, and Q Hou.  A molecular dynamics study of helium bubble formation and gas release near titanium surfaces.  _Journal of nuclear materials_ , 438(1-3):178–182, 2013.
  * Zuckerman & Woolf (2000) Daniel M Zuckerman and Thomas B Woolf.  Efficient dynamic importance sampling of rare events in one dimension.  _Physical Review E_ , 63(1):016702, 2000.

##  Appendix A Appendix

###  A.1 Additional background

History and interpretation of the action. Physicists discovered the action in
the $19^{\textrm{th}}$ century while attempting to derive the mechanics of
solid and fluid bodies on a purely analytic basis Ferguson (2004). There was a
desire to build a coherent, mathematically-consistent framework for describing
the physical world. Notions of physics at the time were highly intuitive and
lacked rigorous analytic treatment. Thus at its conception the utility of the
action was the main focus, not its physical interpretation. Even so, the
interpretation of the action is of fundamental interest in physics. For the
most part (there are exceptions in relativity) the action is written as the
difference between the summed kinetic energy and the summed potential energy
along a (closed) system’s path through time (see Eqn. 1. Our intuition is that
these two terms roughly correspond to the total flow of energy due to entropic
forces and inertial forces respectively. Paths of least action, then,
represent the notion that the dynamics of the universe are the result of a
delicate balance between forces associated with entropy and forces associated
with inertia.

The Newtonian vs. Lagrangian view of physics. When studying physics, one
begins with concepts that are most accessible to human intuition and then
proceeds towards more abstract concepts that describe nature in generality.
This, for example, is why introductory physics courses emphasize a Newtonian
view of physics. Here a group of objects are studied at a particular snapshot
in time, their forces tabulated, and their positions updated to yield the next
state of the system. This view is intuitive because it matches our natural
experience of time and space. The Lagrangian method, by contrast, treats the
entire path of a physical system as its own mathematical entity. It is the
more abstract and unintuitive approach, but it is also far more general and
mathematically consistent across domains.

Related work: gradient-based optimization for physics. We have already
addressed how other works studied minimization of the OM action without using
gradient-based methods. There are also a number of works that have explored
the use of automatic differentiation and gradient-based optimization in the
physical sciences, but have not specifically looked at the action. Schoenholz
& Cubuk (2020) use these tools to study phenomena such as phase transitions
and crystal packing in the context of molecular dynamics. Freeman et al.
(2021) introduce a general-purpose differentiable physics engine for research
in policy optimization. Another recent line of work combines gradient-based
optimization with differential equations (DE). For example, Chen et al. (2018)
proposed using neural networks to parameterize ODEs and Lutter et al. (2019);
Cranmer et al. (2020) showed how to do the same for Lagrangian DEs. Meanwhile,
work by Bar-Sinai et al. (2019) and others has focused on improving adaptive
integration methods.

###  A.2 Additional methods

Optimization. We experimented with standard gradient descent, gradient descent
with momentum, and the Adam optimizer Kingma & Ba (2014), which adaptively
scales momentum. Adding momentum – either directly or via Adam – improved
convergence time by at least a factor of two. We chose the Adam optimizer so
as to avoid setting separate momentum hyperparameters. We used 500
optimization steps for all of our experiments except for the three body
problem, which recovered, with near-exact precision, the ODE path after 1000
steps.

Constructing the baseline paths. To simulate a classical system’s dynamics
using traditional methods, we must first derive a set of second-order
differential equations of motion from physical principles. This is generally
done via application of Newton’s Second Law or the Euler-Lagrange equation
(Eqn. 2). In all systems discussed in this paper, these equations of motion
are directly derivable and expressed as sets of coupled ordinary differential
equations (ODEs) in Table 1. In order to set up a traditional baseline for
contrasting with the new action optimization method, each of these was solved
numerically using Euler’s method prior to optimizer trials, using boundary
conditions known to produce complex behaviors. Euler’s method is the most
basic and broadly understood numerical ODE-solving technique, and was chosen
so as to maximize the baseline’s analogy to our equally basic but as-yet
untested optimization-based method.

Other ways to address the unconstrained energy effect. Several alternative
methods were implemented in attempts to address the unconstrained energy
effect before the path perturbation method discussed in the main body of the
paper was finally selected. Initially, upon realizing the problem lay in a
lack of sufficient constraints, we attempted to feed the optimizer the same
${\bf x}(t_{0})$, $\dot{{\bf x}}(t_{0})$ information that the ODE solver used
to resolve total system energy. The simplest way to do this in our discretized
path model was to give the optimizer more reference points from the baseline
ODE simulation. Practically, this meant freezing not just the initial and
final points of the path during optimization, but several of the adjacent
points as well. While this did constrain the initial and final velocities of
the optimized paths, it did not significantly affect the behavior of the
entire path unless large swathes of points were frozen. As this would greatly
reduce the space of potential applications, we shifted our attention to other
possible solutions.

Next, two distinct regularization schemes were implemented, both based on the
addition of energy conservation terms in the optimizer loss function. The
first added a global energy loss term, which penalized the optimizer in
proportion to the difference between total system energy at each point along a
path and the expected value, calculated from ODE ${\bf x}(t_{0})$ and
$\dot{{\bf x}}(t_{0})$ input conditions. The second regularizer did not rely
on the ODE’s initial conditions, but added a loss term proportional to the
derivative of total system energy over time, i.e. a local energy conservation
term. The first of these worked to a degree, but the path perturbation method
gave far cleaner results on the majority of systems tested. Additionally, path
perturbation showcases the potential ability of the optimizer method to refine
energy predictions based on coarser models. Especially in chemistry, precise
energy predictions are a common target for computational models. The local
energy conservation method, on the other hand, provided no noticeable
advantage. Indeed, it sometimes had an overtly adverse effect on optimizer
performance, likely due to the additional restrictions it enforced on the
space of potential intermediate paths.

For the purposes of our exploratory work, the path perturbation method was the
best fit. We suspect that there is a better method for properly constraining
the total energy and pointing action minimization towards one non-degenerate
path. One of the main focuses of future work in this area should be the
characterization of such a method.

System | Lagrangian $\mathcal{L}=T-V$ | Equations of motion | Noise $\sigma$ | $\Delta t$ | LR
---|---|---|---|---|---
Free Body | $\frac{1}{2}m\dot{x}^{2}-mgx$ | $\ddot{x}=g$ | $1.5$ | $0.25$ | $1$
Pendulum | $\frac{1}{2}ml^{2}\dot{\theta}^{2}-mgl(1-cos(\theta))$ | $\ddot{\theta}=-\frac{g}{l}sin(\theta)$ |  $2e$-$1$ | $1$ |  $5e$-$2$
Dbl. Pend. | Appendix | Appendix |  $6e$-$1$ |  $6e$-$2$ |  $1e$-$2$
Three Body | $\sum\limits_{i,j}\frac{1}{2}m_{i}\bf\dot{x}_{i}^{2}+\frac{Gm_{i}m_{j}}{|r_{ij}|}$ | $\bf\ddot{x}_{i}=-G\sum\limits_{j}\frac{m_{j}}{|r_{ij}|^{2}}\hat{r}_{ij}$ |  $3e$-$2$ | $0.5$ |  $2e$-$4$
Gas | $\sum\limits_{i,j}\frac{1}{2}m_{i}\bf\dot{x}_{i}^{2}+V_{LJ}(\bf r_{ij})$ | $\bf\ddot{x}_{i}=\frac{1}{m_{i}}\sum\limits_{j}F_{LJ}(\bf r_{ij})$ |  $1e$-$2$ | $0.5$ |  $1e$-$4$
Ephemeris | $\sum\limits_{i,j}\frac{1}{2}m_{i}\bf\dot{x}_{i}^{2}+\frac{Gm_{i}m_{j}}{r_{ij}}$ | $\bf\ddot{x}_{i}=-G\sum\limits_{j}\frac{m_{j}}{r_{ij}^{2}}\hat{r}_{ij}$ |  $2e$+$10$ |  $1$ day |  $1e$+$9$
Table 1: All optimizations were run for 500 steps. LR = learning rate.

Double pendulum Lagrangian:

$\bf
L=\frac{1}{2}(m_{1}+m_{2})l_{1}^{2}\dot{\theta}_{1}^{2}\;+\;\frac{1}{2}m_{2}l_{2}^{2}\dot{\theta}_{2}^{2}\;+\;m_{2}l_{1}l_{2}\dot{\theta}_{1}\dot{\theta}_{2}cos(\theta_{1}-\theta_{2})\;+\;(m_{1}+m_{2})gl_{1}cos(\theta_{1})\;+\;m_{2}gl_{2}cos(\theta_{2})$

Double pendulum equations of motion:

$\ddot{\theta}_{i}=\frac{f_{i}-\alpha_{i}f_{j}}{1-\alpha_{i}\alpha_{j}}\quad\textrm{where}\quad
f_{1}\equiv-\frac{l_{2}}{l_{1}}[\frac{m_{2}}{m_{1}+m_{2}}\dot{\theta}_{2}^{2}sin(\theta_{1}-\theta_{2})-\frac{g}{l_{1}}sin(\theta_{1})]$
$f_{2}\equiv\frac{l_{1}}{l_{2}}\dot{\theta}_{1}^{2}sin(\theta_{1}-\theta_{2})-\frac{g}{l_{2}}sin(\theta_{2})$

Lennard-Jones Potential:

$V_{LJ}(\bf
r)=4\epsilon\left[\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^{6}\right]$

where r is the interparticle distance, $\epsilon$ is the depth of the
potential energy well in the interaction between the particles, and $\sigma$
is related to the value of r where the potential energy is at the minimum.
Using the formula $F=-\nabla U$, we can also derive a Lennard-Jones force:

$F_{LJ}(\bf
r)=4\epsilon\left[\left(\frac{12\sigma^{12}}{r^{13}}\right)-\left(\frac{6\sigma^{6}}{r^{7}}\right)\right]$

| Sum Action S (Js) | Sum Kinetic T (Js) | Sum Potential V (Js) | Mean Square Error
---|---|---|---|---
System | ODE | Sim | ODE | Sim | ODE | Sim | Initial | Final
Free Body | -13.8 | -14.0 | 9.54 | 8.84 | 23.4 | 22.8 | 630 | 0.352
Pendulum | -16.3 | -16.3 | 81.9 | 81.7 | 98.2 | 98.0 | 0.118 | 1.20e-2
Dbl Pend | 2.59 | 2.59 | 1.58 | 1.55 | -1.01 | -1.04 | 0.170 | 3.22e-3
3 Body | 6.94e-4 | 6.96e-4 | 2.48e-4 | 2.49e-4 | -4.46e-4 | -4.47e-4 | 9.74e-3 | 5.74e-6
Gas | 7.81e-5 | 7.71e-5 | 1.11e-4 | 1.18e-4 | 3.30e-5 | 4.13e-5 | 3.09e-5 | 2.80e-6
Ephemeris | 3.76e33 | 3.75e33 | 1.28e33 | 1.26e33 | -2.48e33 | -2.49e33 | 1.50e20 | 7.12e17
Table 2: Quantitative results for the six classical systems described in
Section 3.

Figure 5: Visualizing the capped Lennard-Jones potential. The capped portion
of the potential has a linearized slope to prevent forces associated with
small radii from growing extremely large and creating numerical instabilities.

Figure 6: Visualizing solid and gas states generated by using two different
heat dissipation coefficients.

Gas simulation. Molecular dynamics (MD) simulations enable investigation of
dynamic systems that would be difficult to observe experimentally
(Hollingsworth & Dror, 2018). These types of simulations date back to the
1950s, with the earliest work studying hard sphere gasses (Alder & Wainwright,
1959). Today, MD is used to understand far more complicated systems, from
helium diffusion in titanium to protein folding (Zhang et al., 2013; Lee et
al., 2017). These lines of work have advanced in step with our computing
abilities. At their most basic and traditional, MD simulations operate by
taking the locations of all atoms in the system, using some potential to model
interactions between the atoms, calculating each pairwise interaction, and
updating the velocities and positions of the atoms accordingly (Schroeder,
2015).

Our simulation is similar to early MD simulations and simple MD simulations
such as Schroeder (2015) that are often used to help students understand the
properties of gasses in introductory chemistry and physics classes. This
simulation is not optimized for a particular material or problem, so there are
limits to the quantitative information that we can get from it. We use a
Lennard-Jones potential (Figure 6) in determining interaction forces for
simplicity, but this only describes pairwise interactions. More sophisticated
force fields are possible. Even with this simple simulation, we are able to
generate a gas-like state and a solid-like state (Figure 6), akin to those
possible with other simple MD simulations that use a Lennard-Jones potential
(Schroeder, 2015; Sweet et al., 2018).

The Ephemeris experiment. We downloaded raw ephemeris data for the inner
planets of the solar system for the calendar year 2022 (1 day resolution). To
do this we used the online interface provided by NASA’s Horizons project. We
used the Solar System Barycenter (SSB) for a coordinate center. In
constructing our simulation, we used the simple gravitational well potential
of $\frac{Gm_{i}m_{j}}{r_{ij}}$ and used SI units for the gravitational
constant $G$, planet masses $m$, and durations of time. In exploring whether
action minimization could reconstruct the inner planetary orbits, we perturbed
only the paths of the inner planets and not that of the Sun. We considered a
time duration of two months for this experiment because the orbits of Venus
and Mercury cycle more rapidly than that of the Earth (qualitative
visualizations of initial and final paths grew difficult as their orbits begin
to extend over more than one cycle and overlap their tails).

###  A.3 Additional quantum

The principle of least action as a limiting case of the quantum path integral.
The principle of least action is just the limiting case of a path integral
where the actions of paths considered are allowed to approach classical
scales. In this regime, paths whose neighbors have variation in their actions
contribute phases which are out of sync with one another. Summing
contributions from these regions is akin to integrating a rapidly-oscillating
sine function: every section with a positive (constructive) contribution is
inevitably averaged out by an adjacent negative (destructive) contribution. As
spatial scales increase, oscillations along any given path become more
frequent and the net contribution of any arbitrary group of paths approaches
zero. The sole exception occurs where groups of nearby paths have vanishing
variations in action. These paths contribute similar phases and interfere
constructively even at large spatial scales. From this point of view, the
classical principle of least action is actually a description of the resonance
modes of the quantum path integral.

Additional discussion related to Figure 4. The quantum-classical resonance
limit is depicted graphically in Figure 4, where an identical quantum
simulation is run at three different spatial scales. Notice that as we trend
up in scale towards classicality, even small movements along off-diagonals
result in a series of rapid oscillations. These are high-action paths, linear
traversals of a particle between distant points in our space grid over a
vanishing dt. Contributions from these paths will interfere destructively with
one another in almost every trial. We see the direct result of this in the
timed snapshots of each simulation; at a larger scale, the same initial wave
packet will decohere slowly and exhibit only fine high-frequency interference
structures compared to its lower-frequency and thus more-quantum counterparts.
As we approach the classical limit, the entire wave packet will follow the
classical path of least action for a particle bouncing off a wall with
increasing fidelity.


# Probabilistic Monocular 3D Human Pose Estimation with Normalizing Flows

Tom Wehrbein1 Marco Rudolph1 Bodo Rosenhahn1 Bastian Wandt2 1Leibniz
University Hannover, 2University of British Columbia
<EMAIL_ADDRESS>
###### Abstract

3D human pose estimation from monocular images is a highly ill-posed problem
due to depth ambiguities and occlusions. Nonetheless, most existing works
ignore these ambiguities and only estimate a single solution. In contrast, we
generate a diverse set of hypotheses that represents the full posterior
distribution of feasible 3D poses. To this end, we propose a normalizing flow
based method that exploits the deterministic 3D-to-2D mapping to solve the
ambiguous inverse 2D-to-3D problem. Additionally, uncertain detections and
occlusions are effectively modeled by incorporating uncertainty information of
the 2D detector as condition. Further keys to success are a learned 3D pose
prior and a generalization of the best-of-M loss. We evaluate our approach on
the two benchmark datasets Human3.6M and MPI-INF-3DHP, outperforming all
comparable methods in most metrics. The implementation is available on
GitHub111https://github.com/twehrbein/Probabilistic-Monocular-3D-Human-Pose-
Estimation-with-Normalizing-Flows.

##  1 Introduction

Estimating the 3D pose of a human from a single monocular image is an active
research field in computer vision. It has many applications _e.g_. in human
computer interaction, animation, medicine and surveillance. A common approach
is to decouple the problem into two stages. In the first stage, a 2D pose
detector is used to estimate 2D keypoints which are then lifted to 3D joint
locations in the second stage. By utilizing a 2D pose detector pretrained on
diverse and richly annotated data, the 3D pose estimator becomes invariant to
different scenes varying in lighting, background and clothing. However,
reconstructing the correct 3D pose from 2D joint detections is a highly ill-
posed problem because of depth ambiguities and occluded body parts. While some
ambiguities can be resolved by utilizing information from the image (_e.g_.
difference in shading due to depth disparity) or by exploiting known
proportions of the human body, such as joint angle and bone length
constraints, there still remain scenarios where multiple plausible 3D poses
are consistent with the same image. Fig. 1 shows such a situation where the
left arm is occluded by the upper body and therefore its position cannot be
determined unambiguously. Nevertheless, most existing works ignore the
ambiguities by assuming that only a single solution exists. In contrast, we
model monocular 3D human pose estimation as an ambiguous inverse problem with
multiple feasible solutions. Thus, in this work, we propose to estimate the
full posterior distribution of plausible 3D poses conditioned on a monocular
image.

Figure 1: Our model generates diverse 3D pose hypotheses that are consistent
with the input image. Compared to [28, 40] we achieve a higher diversity
mainly where 2D detections are uncertain, in this case for the occluded left
arm. For visualization purposes, more than three hypotheses are shown only for
the highly ambiguous left arm.

Recently, few methods [20, 28, 29, 35, 40] have been proposed that follow the
line of research to explicitly generate multiple 3D pose hypotheses from the
2D input. However, they only consider 2D joint coordinates and ignore the
uncertainty of the 2D detector. While it is reasonable to infer depth
ambiguities based on 2D coordinates only, directly modeling occlusions and
uncertain detections is not meaningful. Fortunately, most 2D human joint
detectors encode valuable information about uncertainties of the location of
human joints in the predicted heatmaps. Instead of discarding this
information, we propose to explicitly extract and utilize the uncertainties of
the 2D detector from the estimated heatmaps. As shown in Fig. 1, this enables
us to effectively model the uncertainties of the 2D detector together with the
inherent depth ambiguities.

In this work, we propose a normalizing flow based method inspired by the
framework for solving ambiguous inverse problems from Ardizzone _et al_. [3].
A normalizing flow [38, 46, 47] is a sequence of bijective transformations
which allows evaluation in both directions. We propose to view 3D human pose
estimation from a single image as an ambiguous inverse problem, since it is a
deterministic forward process (_i.e_. projection of the 3D pose to 2D) with
multiple different inverse mappings. Constructing a bijection between a 3D
pose and the combination of a 2D pose with a latent vector allows to utilize
the 3D-to-2D mapping (forward process) during training. Intuitively, depth
information that otherwise gets lost in the forward process is encoded in the
latent vector. Repeatedly sampling the latent vector and computing the inverse
path of the normalizing flow generates arbitrary many 3D pose hypotheses that
approximate the true posterior distribution. To incorporate the uncertainty
information from the heatmaps, we employ a conditional variant of normalizing
flows [4, 52]. We extract the uncertainty information by fitting 2D Gaussians
to the heatmaps which are then used to form a conditioning vector. We optimize
the model in both directions. The forward path learns the 3D-to-2D mapping and
to produce latent vectors following a predefined distribution. For the inverse
path, we utilize the 3D pose discriminator of [49] to penalize
anthropometrically unfeasible poses. Additionally, we apply a loss enforcing
the 3D pose hypotheses to reflect the uncertainties of the 2D detector.
Motivated by common practice in particle filters, we further propose a
generalization of the best-of-M loss [16] that minimizes the distance between
the mean of the $k$ best hypotheses and the corresponding ground truth.

We evaluate our approach on the two benchmark datasets Human3.6M [19] and MPI-
INF-3DHP [33] and outperform all comparable methods in most metrics. Given the
focus on ambiguous examples, we further evaluate on a subset of Human3.6M
containing only samples with a high degree of 2D detector uncertainty. On this
subset, our method outperforms the competitors by a large margin. To
summarize, our contributions are:

  * •

To the best of our knowledge, we are the first to employ a normalizing flow
based method for modeling the posterior distribution of 3D poses given a
single image.

  * •

Uncertainty information from the predicted heatmaps of the 2D detector is
incorporated into our method, enabling to effectively model occlusions and
uncertain detections.

  * •

We propose a generalization of the best-of-M loss that noticeably improves
prediction performance.

##  2 Related Work

In this section, we first give an overview of recent work in 3D human pose
estimation focused on two-stage approaches. Afterwards, existing methods for
multi-hypotheses 3D pose generation are discussed, followed by an overview of
relevant work on normalizing flows. While there has been recent interest in
estimating the 3D human body shape from monocular images [5, 21, 25, 26, 30,
37, 53, 55], this work focuses on predicting the 3D locations of a set of
predefined joints.

Lifting 2D to 3D: Our approach belongs to the vast body of work that estimate
3D poses from the output of a 2D pose detector [8, 9, 12, 17, 18, 31, 41, 49,
50, 51, 54]. These two-stage approaches decouple the difficult problem of 3D
depth estimation from the easier 2D pose localization. Furthermore, it allows
to use both indoor and in-the-wild data for training the 2D detector, which
effectively reduces the bias towards sterile indoor scenes. Akhter and Black
[1] learn a pose-conditioned joint angle limit prior to restrict invalid 3D
pose reconstructions. They perform 3D pose estimation using an over-complete
dictionary of poses. Moreno-Noguer [34] casts the problem as a regression
between 2D and 3D poses represented as distance matrices. Lifting 2D to 3D
joints was further sparked by Martinez _et al_. [32], who employ a simple
fully-connected network to lift 2D detections to 3D poses, surprisingly
outperforming past approaches. Due to its simplicity and strong performance,
it serves as a popular baseline for many following works.

Unlike the above mentioned approaches that assume a unimodal posterior
distribution and only predict a single 3D pose for each input, we are able to
generate a diverse set of plausible 3D poses. Additionally, anatomical
constraints are learned implicitly by utilizing a strong 3D pose
discriminator. In contrast to previous works that integrate uncertainty
information of the 2D detector (_e.g_. [7, 50, 54]), we fit a 2D Gaussian to
each heatmap instead of using only the maximum value of each heatmap as
confidence score, thus better capturing the uncertainty distribution.

Figure 2: An overview of our proposed method. We employ a normalizing flow
consisting of affine coupling blocks [11] for generating multiple 3D pose
hypotheses. By constructing a bijection between a 3D pose and the
concatenation of a 2D pose with a latent vector, we can exploit the 3D-to-2D
mapping (forward path) during training. The model is optimized in both
directions, whereas at inference, only the path from 2D to 3D (inverse path)
is computed. Arbitrary many 3D pose hypotheses can be generated by repeatedly
sampling the latent vector from a known distribution and computing the inverse
path. Uncertainty information of the 2D detector in form of fitted Gaussians
is incorporated by conditioning the coupling blocks. The architecture of a
single coupling block is visualized in the gray box. For visualization
purposes, only the forward computation of the coupling block is shown.

Multi-Hypotheses 3D Human Pose Estimation: There are early works [27, 42, 43,
44] that extensively analyze and discuss the ambiguities of monocular 3D human
pose estimation and sample multiple 3D poses via heuristics. More recently,
Jahangiri and Yuille [20] propose to generate multiple hypotheses from a
predicted seed 3D pose by uniformly sampling from learned occupancy matrices
[1]. Furthermore, they impose bone length constraints and reject hypotheses
with 2D reprojection error larger than some threshold. Li and Lee [28] employ
a mixture density network (MDN) [6] to learn the multimodal posterior
distribution. The conditional mean of each Gaussian kernel then denotes one 3D
pose hypothesis. Oikarinen _et al_. [35] utilize the semantic graph neural
network of [56] to improve upon the MDN approach of [28]. Contrary to our
normalizing flow based approach, the number of generated hypotheses needs to
be specified a priori and is fixed for every input. Furthermore, when
increasing the number of generated hypotheses, significantly more
computational resources are required. Sharma _et al_. [40] employ a
conditional variational autoencoder to synthesize diverse 3D pose hypotheses
conditioned on a 2D pose detection. They also propose to derive joint-ordinal
depth relations from the image to rank the 3D pose samples. In contrast to
[20], our normalizing flow based approach does not need to incorporate
computationally heavy rejection sampling or requires to define the number of
generated 3D pose hypotheses a priori. Our method is more flexible and is able
to model any posterior distribution without requiring explicit hard
constraints. Moreover, we are the only ones to incorporate the uncertainty
information of the 2D detector, enabling us to significantly improve on highly
ambiguous cases and to inherently handle an arbitrary number of occluded
joints.

Normalizing Flows: A normalizing flow [38, 46, 47] is a sequence of bijective
transformations that transforms a simple tractable distribution into a complex
target data distribution. Because of the bijectivity, evaluation in both
directions is possible. Namely, sampling data from the modeled distribution as
well as exact density estimation (_i.e_. assigning a likelihood to each data
point). Most common state-of-the-art flow architectures are based on auto-
regressive models that utilize the Bayesian chain rule to decompose the
density [10, 11, 13, 23, 36, 39]. For a more comprehensive introduction, we
refer the reader to [24].

Ardizonne _et al_. [3] extend the real-valued non-volume preserving (Real-NVP)
transformations from Dinh _et al_. [11] to the task of computing posteriors
for ambiguous inverse problems. Given such an ambiguous inverse problem, they
propose to learn the well-understood forward process in a supervised manner
and encode otherwise lost information in additional latent variables. Thus,
they learn a bijective mapping between the target data distribution and the
joint distribution of latent variables and forward process solutions. Due to
invertibility, the inverse is implicitly learned. By repeatedly sampling the
latent variables from a simple tractable distribution, they can approximate
the full posterior. Inspired by their work, we adopt and extent their
framework for modeling the full posterior distribution of plausible 3D poses
conditioned on a monocular image. We introduce a conditioning vector, a
learnable prior and two additional loss functions.

To the best of our knowledge, the only previous works in human pose estimation
utilizing normalizing flows are [5, 53, 55]. However, they employ normalizing
flows as 3D pose prior and not for directly modeling the posterior
distribution of 3D poses conditioned on an image.

##  3 Method

Our aim is to learn the full posterior distribution of plausible 3D poses
conditioned on a monocular image. We follow the popular two-stage approach by
first applying a state-of-the-art 2D joint detector [45] and subsequently
using its output to estimate corresponding 3D pose hypotheses. The core idea
is that instead of conditioning the posterior distribution only on the 2D
detections, we additionally utilize uncertainty information extracted from the
predicted heatmaps in a novel way. This enables to effectively model the
uncertainties of the 2D detector together with the inherent depth ambiguities.

An overview of the proposed method is shown in Fig. 2. To learn the posterior
distribution, we employ a normalizing flow to construct a bijective mapping
between a 3D pose $\bm{x}\in\mathbb{R}^{3J}$ and the concatenation of a 2D
pose $\bm{y}\in\mathbb{R}^{2J}$ with a latent vector
$\bm{z}\in\mathbb{R}^{J}$, where $J$ is the number of joints in one pose. The
introduction of the latent vector $\bm{z}$ allows to utilize the well-defined
forward process of projecting a 3D pose to its 2D observation during training.
Intuitively, $\bm{z}$ captures depth information that is otherwise lost in the
mapping from 3D to 2D. Instead of simply using the argmax of the heatmaps, we
incorporate the uncertainty information of the 2D detector by conditioning the
normalizing flow on Gaussians fitted to the heatmaps. At inference, the full
posterior is approximated by repeatedly sampling $\bm{z}$ from the
distribution of latent variables and computing the inverse path. If the
forward process is simulated successfully, all generated hypotheses reproject
to the corresponding 2D pose observation.

###  3.1 Conditional Normalizing Flow

As normalizing flow we adopt the Real-NVP [11] affine coupling block
architecture. This architecture can straight-forwardly be extended to
incorporate a conditional input [4, 52]. A single coupling block is shown in
the gray box in Fig. 2. The input $\bm{u}_{\mathrm{in}}$ is split into two
parts $\bm{u}_{\mathrm{in},1}$ and $\bm{u}_{\mathrm{in},2}$. Subsequently,
$\bm{u}_{\mathrm{in},1}$ and $\bm{u}_{\mathrm{in},2}$ undergo a scale and
translation transformation parameterized by the functions $s_{i}$ and $t_{i}$
($i\in\\{1,2\\}$) on two separate paths. The outputs
$\bm{u}_{\mathrm{\mathrm{out}},1}$ and $\bm{u}_{\mathrm{\mathrm{out}},2}$ are
concatenated to form the overall output of the coupling block. Given the
heatmap condition $\hat{\bm{c}}$, further encoded into the conditioning vector
$\bm{c}=h_{\bm{\theta}}(\hat{\bm{c}})$, the forward path of a coupling block
is defined as

$\displaystyle\bm{u}_{\mathrm{\mathrm{out}},2}=\bm{u}_{\mathrm{in},2}\odot
e^{s_{1}(\bm{u}_{\mathrm{in},1},\bm{c})}+t_{1}(\bm{u}_{\mathrm{in},1},\bm{c})$
(1) $\displaystyle\bm{u}_{\mathrm{\mathrm{out}},1}=\bm{u}_{\mathrm{in},1}\odot
e^{s_{2}(\bm{u}_{\mathrm{\mathrm{out}},2},\bm{c})}+t_{2}(\bm{u}_{\mathrm{\mathrm{out}},2},\bm{c}),$

where $\odot$ denotes the element-wise multiplication. The exponential
function is used to prevent multiplication by zero, which ensures the
invertibility of the block. Note that $s_{i}$ and $t_{i}$ represent functions
that do not need to be invertible. The only restriction is that their produced
output matches the dimensions of the data on the corresponding path in the
coupling block. Instead of regressing the scale and translation coefficients
separately, we employ a fully-connected network that jointly predicts them by
splitting its output. By construction, the coupling block can be trivially
inverted without any computational overhead. The overall network consists of
multiple chained blocks, each followed by a predefined random permutation
which shuffles the path assignment of the variables. The output of the last
block is split to form the 2D pose $\bm{y}$ and the latent vector $\bm{z}$.

Following Ardizzone _et al_. [4], we adopt a parameterized soft clamping
mechanism to prevent instabilities caused by the exponential function in the
coupling block. The soft clamping is defined as

$\sigma_{\alpha}(r)=\frac{2\alpha}{\pi}\arctan{\frac{r}{\alpha}},$ (2)

and is applied as the last layer of $s_{1}$ and $s_{2}$. It prevents scaling
components of exploding magnitude by restricting the output to the interval
$(-\alpha,\alpha)$.

###  3.2 Heatmap Condition

Recent 2D detectors are optimized by applying a supervised loss between the
predicted heatmap and a ground-truth heatmap consisting of a 2D Gaussian
centered at the joint location. This leads to the predicted heatmaps being a
valuable source of uncertainty of the 2D detector. Instead of estimating 3D
poses solely based on 2D joint coordinates, we incorporate the uncertainties
of the 2D detector encoded in the estimated heatmaps. Specifically, we fit a
2D Gaussian to each predicted heatmap to best capture the uncertainty
distribution. The fitting process is done using non-linear least squares. As
initial parameters, we set the amplitude to $1$, the mean of each Gaussian to
the corresponding regressed 2D joint location and the covariance matrix to a
diagonal matrix with $\sigma^{2}_{x}=\sigma^{2}_{y}=\sigma^{2}_{\text{gt}}$,
where $\sigma^{2}_{\text{gt}}$ is the ground-truth variance used for training
the 2D detector. For each image, the fitted coefficients are stacked to form a
single vector. We discard the Gaussian coefficients for the hip joints, since
the typical alignment of the root joint of the 3D poses heavily reduces the
possible variances in these joints. Thus, the heatmap conditioning vector is
denoted as $\hat{\bm{c}}\in\mathbb{R}^{6(J-3)}$. We employ a fully-connected
network as encoding network $h_{\bm{\theta}}$ that further encodes
$\hat{\bm{c}}$ into $\bm{c}=h_{\bm{\theta}}(\hat{\bm{c}})$. For the 3D pose
hypotheses to best reflect the uncertainties of the 2D detector, we explicitly
optimize the network to match the 2D Gaussian distributions in the $x$\- and
$y$-direction of the 3D hypotheses for each joint. Let
$\hat{\Sigma}\in\mathbb{R}^{2\times 2}$ be the covariance matrix of a single
joint estimated from the positions of that joint in the $L$ produced
hypotheses. Defining $\Sigma\in\mathbb{R}^{2\times 2}$ as the covariance
matrix of the fitted 2D Gaussian of the corresponding heatmap, we minimize a
masked lower bound Root Mean Square Error (RMSE) between both covariance
matrices:

$\displaystyle\mathcal{L}_{\text{HM}}\big{(}{\Sigma},$
$\displaystyle\hat{\Sigma}\big{)}=m\cdot\Big{(}\max\big{(}0,\Sigma_{1,1}-\hat{\Sigma}_{1,1}\big{)}^{2}$
(3)
$\displaystyle+\max\big{(}0,\Sigma_{2,2}-\hat{\Sigma}_{2,2}\big{)}^{2}+\big{(}\Sigma_{1,2}-\hat{\Sigma}_{1,2}\big{)}^{2}\Big{)}^{\frac{1}{2}},$

where the masking scalar $m$ is defined as

$m=\left\\{\begin{array}[]{ll}1&\sqrt{\Sigma_{1,1}}>\sigma_{t}\lor\sqrt{\Sigma_{2,2}}>\sigma_{t}\\\
0&\text{otherwise}\end{array}\right..$ (4)

Thus, the loss has no influence if the 2D detector is certain about the
location of the specific joint, indicated by a fitted Gaussian with standard
deviations smaller than the threshold $\sigma_{t}$. To not unnecessarily
restrict the network, we only penalize the diagonal entries of the covariance
matrices if they are smaller than the corresponding ground-truth values.

###  3.3 Optimization

The core idea of the optimization procedure is to train the 3D-to-2D mapping
(forward path) in a supervised manner, while the highly ambiguous 2D-to-3D
mapping (inverse path) is learned implicitly due to the invertibility of the
normalizing flow and supported by additional supervision with the inverse
process. Each training iteration consists of first calculating the forward
path, followed by $L$ computations of the inverse path and two additional, one
for the discriminator and one for the deterministic 3D reconstruction. The
gradients from both directions are accumulated before performing a parameter
update. Note that due to the Real-NVP coupling block architecture, both
directions can be computed efficiently.

Forward Path: In the forward process, the network predicts the corresponding
2D joint detections given a 3D pose. This is optimized using the $L_{1}$
distance:

$\mathcal{L}_{\text{2D}}=\left\lVert\bm{y}-\bm{\hat{y}}\right\rVert_{1},$ (5)

where $\bm{y}$ is the ground-truth and $\bm{\hat{y}}$ the estimated 2D
observation. The estimated latent variables are optimized to follow a zero-
mean isotropic Gaussian $p_{\bm{Z}}=\mathcal{N}(0,I)$ and to be independent
from the distribution of 2D observations $p_{\bm{Y}}$. Both properties are
enforced by minimizing the Maximum Mean Discrepancy (MMD) [14] between the
joint distribution of network outputs $q(\bm{y},\bm{z})$ and the product of
marginal distributions $p_{\bm{Y}}$ and $p_{\bm{Z}}$. Given samples
$\hat{\bm{V}}={\\{\hat{\bm{v}}_{i}}\\}^{n}_{i=1}$ drawn i.i.d. from
$q(\bm{y},\bm{z})$ and $\bm{V}={\\{\bm{v}_{i}}\\}^{n}_{i=1}$ with
$\bm{v}_{i}=[\bm{y},\bm{z}]$ and $\bm{y}\sim p_{\bm{Y}}$, $\bm{z}\sim
p_{\bm{Z}}$, the unbiased estimator of the squared MMD with kernel $\varphi$
is

$\displaystyle\mathcal{L}$
$\displaystyle{}_{\text{MMD}}=\text{MMD}^{2}_{u}(\bm{V},\hat{\bm{V}})=\frac{1}{n(n-1)}\sum_{i\neq
j}^{n}\varphi\left(\bm{v}_{i},\bm{v}_{j}\right)$ (6)
$\displaystyle+\frac{1}{n(n-1)}\sum_{i\neq
j}^{n}\varphi\left(\hat{\bm{v}}_{i},\hat{\bm{v}}_{j}\right)-\frac{2}{n^{2}}\sum_{i,j=1}^{n}\varphi\left(\bm{v}_{i},\hat{\bm{v}}_{j}\right).$

Following [3], we block the gradients of $\mathcal{L}_{\text{MMD}}$ with
respect to $\bm{y}$ to prevent the predictions of $\bm{y}$ from deteriorating.

Inverse Path: Given a 2D pose $\bm{y}$, a latent vector $\bm{z}$ is drawn from
the base distribution $p_{\bm{Z}}$ and concatenated to form the input
$[\bm{y},\bm{z}]$ of the inverse path. By repeatedly sampling $\bm{z}\sim
p_{\bm{Z}}$, arbitrary many 3D pose hypotheses can be created. Although
$\mathcal{L}_{\text{2D}}$ and $\mathcal{L}_{\text{MMD}}$ are in theory
sufficient to best approximate the true posterior distribution [3], we apply
additional losses to the inverse path to improve convergence. To penalize
geometrically unfeasible 3D pose hypotheses, we introduce a discriminator
network and adopt the Improved Wasserstein GAN training procedure of [15]. The
inverse path acts as the generator by producing 3D poses that minimize the
negated output of the discriminator. This loss is denoted as
$\mathcal{L}_{\text{gen}}$. The architecture of the discriminator is taken
from [49], including a Kinematic Chain Space layer [48] encoding bone lengths
and angular representations. Additionally, we generate a 3D pose for each 2D
input by using the corresponding latent vector $\bm{z}_{\text{det}}$ produced
in the forward path. Note that contrary to sampling a latent vector
$\bm{z}\sim p_{\bm{Z}}$, when using the estimated latent vector
$\bm{z}_{\text{det}}$, it is reasonable to apply a supervised loss
$\mathcal{L}_{\text{det}}$ linking a 2D input to a single 3D pose, since the
combination of a predicted latent vector and matching 2D pose detection should
correspond to the single exact solution of the ambiguous inverse problem. We
minimize the $L_{1}$ distance between the ground-truth 3D pose $\bm{x}$ and
the estimated 3D pose $\hat{\bm{x}}_{\text{det}}$

$\mathcal{L}_{\text{det}}=\left\lVert\bm{x}-\hat{\bm{x}}_{\text{det}}\right\rVert_{1}.$
(7)

To further guide the optimization process, we propose a generalization of the
best-of-M loss [16]. Given a set of 3D pose hypotheses
$\bm{H}={\\{\hat{\bm{x}}_{i}}\\}^{L}_{i=1}$ generated from the same 2D input,
we select the subset $\bm{H}_{\text{topk}}\subseteq\bm{H}$ consisting of the
$k$ pose hypotheses with the lowest Mean Per Joint Position Error (MPJPE) to
the corresponding ground-truth pose $\bm{x}$. We then minimize the $L_{1}$
distance between the ground-truth pose $\bm{x}$ and the mean of the $k$ best
hypotheses:

$\mathcal{L}_{\text{MB}}=\left\lVert\bm{x}-\frac{\sum_{\hat{\bm{x}}\in\bm{H}_{\text{topk}}}\hat{\bm{x}}}{k}\right\rVert_{1}.$
(8)

Overall: In total, the objective function of our normalizing flow is

$\displaystyle\mathcal{L}_{\text{NF}}=$
$\displaystyle\mathcal{L}_{{\text{2D}}}+\mathcal{L}_{\text{gen}}+\lambda_{\text{MMD}}\mathcal{L}_{\text{MMD}}$
(9)
$\displaystyle+\lambda_{\text{det}}\mathcal{L}_{\text{det}}+\lambda_{\text{MB}}\mathcal{L}_{\text{MB}}+\lambda_{\text{HM}}\mathcal{L}_{\text{HM}},$

where $\lambda_{\text{MMD}}$, $\lambda_{\text{det}}$, $\lambda_{\text{MB}}$
and $\lambda_{\text{HM}}$ represent the weights of the corresponding losses.
The discriminator network is optimized to distinguish between the 3D poses
produced by the normalizing flow and 3D poses from the training set by
minimizing the WGAN-GP objective function [15]. The encoding network
$h_{\bm{\theta}}$ is jointly optimized with the normalizing flow by
propagating gradients from $\mathcal{L}_{\text{NF}}$ through
$h_{\bm{\theta}}$.

##  4 Experiments

###  4.1 Datasets and Evaluation Metrics

Human3.6M [19] is the largest video pose dataset for 3D human pose estimation.
It features $7$ professional actors performing 15 different activities, such
as Sitting, Walking and Smoking. For each frame, accurate 2D and 3D joint
locations and camera parameters are provided. We follow the standard protocols
and evaluate on every $64^{\text{th}}$ frame of subjects 9 and 11. Protocol 1
computes the Mean Per Joint Position Error (MPJPE) between the reconstructed
and ground-truth 3D joint coordinates directly, whereas Protocol 2 first
applies a rigid alignment between the poses (PMPJPE). We additionally show
results for the Correct Poses Score (CPS) metric proposed by [50]. It
considers a pose as correct if and only if all joints have a Euclidean
distance to the ground-truth below a threshold value $\vartheta$. CPS is then
defined as the area under curve for
$\vartheta\in[0\text{\,mm},300\text{\,mm}]$. Instead of evaluating the
reconstruction joint by joint, the CPS considers the whole pose. Compared to
other common metrics, it is better suited to detect wrongly estimated poses
that could negatively influence downstream tasks.

Human3.6M Ambiguous (H36MA): To focus evaluation on highly ambiguous examples,
we select a subset222Information about the exact composition of the subset can
be found in the official GitHub repository. of the Human3.6M test split
according to the uncertainties of the 2D detector. This subset only contains
samples for which at least one fitted Gaussian has a standard deviation larger
than 5 px, which holds true for $6.4\%$ of all samples in the test split.
These samples are extremely challenging since the joint detector gives
inaccurate or wrong results.

MPI-INF-3DHP (3DHP) [33] is a 3D human pose dataset containing annotated
images recorded in three different settings: studio with green screen, studio
without green screen and outdoors. We evaluate on the test split without
utilizing the training data to assess the generalization capability of our
network. Following previous works, the Percentage of Correct Keypoints (PCK)
under 150 mm is adopted as the metric for 3DHP.

###  4.2 Implementation Details

2D Detector: We use the publicly available HRNet [45] pretrained on MPII [2]
as our 2D joint detector and finetune it on Human3.6M. Target ground-truth
heatmaps are created with $\sigma_{\text{gt}}=2\text{\,px}$.

Data Preprocessing: We center each 2D pose to its mean and divide it by its
standard deviation. The 3D poses are processed in metres and also mean
centered individually. Before evaluation, 3D poses are zero-centered around
the hip joint to follow the standard protocols.

Network Details: The normalizing flow consists of 8 coupling blocks with
fully-connected networks, denoted as subnetworks, as scale and translation
functions. Each subnetwork upscales its input to 1024 dimensions with a fully-
connected layer. This is followed by a ReLU and a second fully-connected layer
with dimension 48. The condition encoding network $h_{\bm{\theta}}$ follows
the same design with 256 and 56 as output dimensions of the fully-connected
layers. We set the clamping parameter inside the coupling blocks to
$\alpha=2.0$. For $\mathcal{L}_{\text{MMD}}$, we follow [3] and employ a
mixture of inverse multiquadratics kernels

$\varphi^{im}_{\mathcal{S}}(\bm{v},\hat{\bm{v}})=\sum_{b\in\mathcal{S}}\frac{b}{b+\left\|\bm{v}-\hat{\bm{v}}\right\|^{2}}$
(10)

with bandwidth parameters $\mathcal{S}=\\{0.0025,0.04,0.81\\}$.

Training: The overall network is trained for 155 epochs using Adam [22] with
an initial learning rate of $1\cdot 10^{-4}$ and momentum values
$\beta_{1}=0.5$ and $\beta_{2}=0.9$. The learning rate is halved after 150
epochs, and a batch size of 64 is used. To improve optimization stability, we
clip the gradients in the range $[-15,15]$. During training, the covariance
matrices are computed from $L=200$ 3D pose hypotheses and the standard
deviation threshold for the masking of $\mathcal{L}_{\text{HM}}$ (Eq. 4) is
set to $\sigma_{t}=1.05\cdot\sigma_{\text{gt}}=2.1$. The weights of the
different losses are set to $\lambda_{\text{MMD}}=10$,
$\lambda_{\text{det}}=\lambda_{\text{MB}}=4$ and $\lambda_{\text{HM}}=750$,
and the number of best hypotheses selected in $\mathcal{L}_{\text{MB}}$ to
$k=5$. Since we estimate 3D poses in metric scale, there needs to be a
conversion factor defined to relate between covariance matrices from pose
hypotheses and from heatmaps. We empirically found a good conversion factor to
be $1$ px $\mathrel{\widehat{=}}10$ mm.

Protocol 1 (MPJPE) | Direct. | Disc. | Eat | Greet | Phone | Photo | Pose | Purch. | Sit | SitD | Smoke | Wait | WalkD | Walk | WalkT | Avg.
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
Martinez _et al_. [32] ($M=1$) | 51.8 | 56.2 | 58.1 | 59.0 | 69.5 | 78.4 | 55.2 | 58.1 | 74.0 | 94.6 | 62.3 | 59.1 | 65.1 | 49.5 | 52.4 | 62.9
Li _et al_. [29] ($M=10$) | 62.0 | 69.7 | 64.3 | 73.6 | 75.1 | 84.8 | 68.7 | 75.0 | 81.2 | 104.3 | 70.2 | 72.0 | 75.0 | 67.0 | 69.0 | 73.9
Li _et al_. [28] ($M=5$) | 43.8 | 48.6 | 49.1 | 49.8 | 57.6 | 61.5 | 45.9 | 48.3 | 62.0 | 73.4 | 54.8 | 50.6 | 56.0 | 43.4 | 45.5 | 52.7
Oikarinen _et al_. [35] ($M=200$) | 40.0 | 43.2 | 41.0 | 43.4 | 50.0 | 53.6 | 40.1 | 41.4 | 52.6 | 67.3 | 48.1 | $\bm{44.2}$ | $\bm{44.9}$ | 39.5 | 40.2 | 46.2
Sharma _et al_. [40] ($M=200$) | $\bm{37.8}$ | 43.2 | 43.0 | 44.3 | 51.1 | 57.0 | $\bm{39.7}$ | 43.0 | 56.3 | 64.0 | 48.1 | 45.4 | 50.4 | 37.9 | $\bm{39.9}$ | 46.8
Ours ($\bm{z}_{0}$) ($M=1$) | 52.4 | 60.2 | 57.8 | 57.4 | 65.7 | 74.1 | 56.2 | 59.1 | 69.3 | 78.0 | 61.2 | 63.7 | 67.0 | 50.0 | 54.9 | 61.8
Ours ($M=200$) | 38.5 | $\bm{42.5}$ | $\bm{39.9}$ | $\bm{41.7}$ | $\bm{46.5}$ | $\bm{51.6}$ | 39.9 | $\bm{40.8}$ | $\bm{49.5}$ | $\bm{56.8}$ | $\bm{45.3}$ | 46.4 | 46.8 | $\bm{37.8}$ | 40.4 | $\bm{44.3}$
Protocol 2 (PMPJPE) | Direct. | Disc. | Eat | Greet | Phone | Photo | Pose | Purch. | Sit | SitD | Smoke | Wait | WalkD | Walk | WalkT | Avg.
Martinez _et al_. [32] ($M=1$) | 39.5 | 43.2 | 46.4 | 47.0 | 51.0 | 56.0 | 41.4 | 40.6 | 56.5 | 69.4 | 49.2 | 45.0 | 49.5 | 38.0 | 43.1 | 47.7
Li _et al_. [29] ($M=10$) | 38.5 | 41.7 | 39.6 | 45.2 | 45.8 | 46.5 | 37.8 | 42.7 | 52.4 | 62.9 | 45.3 | 40.9 | 45.3 | 38.6 | 38.4 | 44.3
Li _et al_. [28] ($M=5$) | 35.5 | 39.8 | 41.3 | 42.3 | 46.0 | 48.9 | 36.9 | 37.3 | 51.0 | 60.6 | 44.9 | 40.2 | 44.1 | 33.1 | 36.9 | 42.6
Oikarinen _et al_. [35] ($M=200$) | 30.8 | 34.7 | 33.6 | 34.2 | 39.6 | 42.2 | 31.0 | 31.9 | 42.9 | 53.5 | 38.1 | 34.1 | 38.0 | 29.6 | 31.1 | 36.3
Sharma _et al_. [40] ($M=200$) | 30.6 | 34.6 | 35.7 | 36.4 | 41.2 | 43.6 | 31.8 | 31.5 | 46.2 | 49.7 | 39.7 | 35.8 | 39.6 | 29.7 | 32.8 | 37.3
Ours ($\bm{z}_{0}$) ($M=1$) | 37.8 | 41.7 | 42.1 | 41.8 | 46.5 | 50.2 | 38.0 | 39.2 | 51.7 | 61.8 | 45.4 | 42.6 | 45.7 | 33.7 | 38.5 | 43.8
Ours ($M=200$) | $\bm{27.9}$ | $\bm{31.4}$ | $\bm{29.7}$ | $\bm{30.2}$ | $\bm{34.9}$ | $\bm{37.1}$ | $\bm{27.3}$ | $\bm{28.2}$ | $\bm{39.0}$ | $\bm{46.1}$ | $\bm{34.2}$ | $\bm{32.3}$ | $\bm{33.6}$ | $\bm{26.1}$ | $\bm{27.5}$ | $\bm{32.4}$

Table 1: Detailed results of MPJPE in millimetres on Human3.6M under Protocol
1 (no rigid alignment) and Protocol 2 (rigid alignment). Our model achieves
state-of-the-art results, outperforming all other methods in nearly every
activity. All scores are taken from the referenced papers, except the row
marked with * which is computed using the publicly available official code and
model from [40]. The number of samples estimated by the respective approaches
is denoted as $M$.

###  4.3 Evaluation on Human3.6M

We follow previous works and report metrics for the best 3D pose hypothesis
generated by our network. This is especially reasonable for ambiguous
examples, where multiple diverse 3D poses form a correct solution for the 3D
pose reconstruction. Therefore, instead of validating whether predictions are
equal to a specific solution, we evaluate if that specific solution is
contained in the set of predictions. Additionally, we show results for 3D
poses generated with an all-zero latent vector $\bm{z}_{0}$. Since we sample
$\bm{z}$ from $\mathcal{N}(0,I)$ during training, such poses are approximately
the highest likelihood solutions. Following [40], we produce $M=200$
hypotheses for each 2D input. The results of our approach and other state-of-
the-art methods are shown in Table 1. We outperform every competitor and
achieve a clear improvement of $4.1\%$ and $10.7\%$ over the previous best
scores under Protocol 1 and Protocol 2. Note that Li _et al_. [28] only show
detailed results for $M=5$, but state that their model performance does not
significantly improve when increasing $M$. We generated the numbers for [40]
under Protocol 2 (row marked with *) using their publicly available model,
code and data, because they only report scores for the PMPJPE on subject 11.
Outperforming the single prediction baseline of [32] with $\bm{z}_{0}$
generated poses (_i.e_. $M=1$) shows that our model is additionally able to
give strong single predictions.

To evaluate the performance on highly ambiguous examples, we compute results
for the challenging subset H36MA. We use the publicly available code from [40]
and [28] to compare with their approaches. As is shown in Table 2, we
outperform both competitors significantly and by a larger margin than on the
whole test set. This emphasizes the ability of our model to generate diverse
hypotheses for highly ambiguous examples. We argue that the CPS is especially
meaningful in this setting, since high individual joint errors that often
occur for challenging poses cannot be averaged out as in _e.g_. MPJPE or PCK.

Method | MPJPE$\downarrow$ | PMPJPE$\downarrow$ | PCK$\uparrow$ | CPS$\uparrow$
---|---|---|---|---
Li _et al_. [28] | 81.1 | 66.0 | 85.7 | 119.9
Sharma _et al_. [40] | 78.3 | 61.1 | 88.5 | 136.4
Ours | $\bm{71.0}$ | $\bm{54.2}$ | $\bm{93.4}$ | $\bm{171.0}$

Table 2: Evaluation results on the subset H36MA containing highly ambiguous
examples. For each metric, the best hypothesis score is reported.

###  4.4 Transfer to MPI-INF-3DHP

To assess the generalization capability of our model, we evaluate on MPI-
INF-3DHP. Note that neither the 2D detector nor the normalizing flow is
trained on this dataset. The results are shown in Table 3. Even though [28]
use the ground-truth 2D joints provided by the dataset, we clearly outperform
them in all three settings. We also achieve competitive performance compared
to the weakly supervised approach from [29] that focuses on transfer learning.
Our strong results for outdoor scenes further emphasize the generalization
capability to different settings.

Method | Studio GS | Studio no GS | Outdoor | All PCK
---|---|---|---|---
Li _et al_. [29] | $\bm{86.9}$ | $\bm{86.6}$ | 79.3 | $\bm{85.0}$
Li _et al_. [28] | 70.1 | 68.2 | 66.6 | 67.9
Ours | 86.6 | 82.8 | $\bm{82.5}$ | 84.3

Table 3: Quantitative results on MPI-INF-3DHP. We outperform the approach from
[28] by a large margin which even uses ground truth 2D joint positions. Note
that [29] is trained weakly supervised and therefore specifically built for
transfer learning. However, we still achieve on par results and even
outperform them in the challenging outdoor sequences.

###  4.5 Sample Diversity

Heatmap Variance: We visually inspect the distribution of generated joint
locations and compare them with the corresponding fitted Gaussians in Fig. 3.
For visualization purposes, only three hypotheses are shown for all joints
except the one with the highest uncertainty. As can be seen, the uncertainties
of the 2D detector are reflected in the 3D hypotheses.

Figure 3: The uncertainties of the 2D detector are successfully reflected in
the 3D pose hypotheses. For visualization purposes, we only show the fitted
Gaussian and a high number of hypotheses for the joint with highest
uncertainty.

Depth Ambiguities: Even though variance in the depth direction is not
explicitly optimized, our model learns to generate feasible hypotheses with
varying depth. In fact, the standard deviation of the hypotheses averaged over
all joints in the test set of Human3.6M is highest in the depth direction with
$42.4$ mm, compared to $18.3$ mm and $17.3$ mm in the $x$\- and
$y$-directions. Ankle, elbow and wrist joints account for the highest amount
of variance. A visual example of meaningful depth diversity is given in Fig.
4.

Figure 4: Depth ambiguities can be modeled together with the uncertainties of
the 2D detector.

Sample Set Size and Noise Baseline: In Fig. 5, we plot the MPJPE on the subset
H36MA with increasing number of samples. The best hypothesis performance of
our model continues to improve significantly, further enlarging the gap to
[40]. To validate that our approach is superior to directly sampling from the
fitted Gaussians, we also plot results for a sampling baseline. The baseline
is constructed by adding noise sampled from the fitted Gaussians to each joint
of the $\bm{z}_{0}$ predictions. A constant Gaussian for the depth dimension
is assumed. Evidently, the performance of this baseline saturates earlier and
at a higher error.

$0$$200$$400$$600$$800$$1{,}000$$65$$70$$75$$80$$85$$90$$95$# hypothesesMPJPE
(mm)$\bm{z}_{0}$ \+ noiseSharma _et al_. [40]Ours Figure 5:  Evaluation
results on the subset H36MA for an increasing number of hypotheses. Our model
further improves and enlarges the gap to [40] and to a noise baseline.

###  4.6 Ablation Studies

To quantify the influence of our proposed components and loss functions, we
remove them individually and show the results in Table 4. As can be seen, the
removal of each component leads to a degradation in performance. When removing
the heatmap condition, a large drop in the CPS can be observed. This shows
that some individual joints cannot be reconstructed without the uncertainty
information of the 2D detector. Providing the condition alone already leads to
a significant improvement of the CPS, indicating that the network can
automatically leverage the information to model ambiguities. Adding
$\mathcal{L}_{\text{HM}}$ further improves all metrics. The importance of the
discriminator becomes especially evident when considering the worst hypothesis
error instead of the best. For example, Protocol 2 computed for the worst
hypothesis deteriorates from $86.8$ mm to $284.1$ mm without the
discriminator. Thus, the adversarial training procedure ensures the
feasibility of the generated poses. Table 4 also shows the influence of the
number of best hypotheses $k$ selected for computing
$\mathcal{L}_{\text{MB}}$. Note that $\mathcal{L}_{\text{MB}}$ with $k=1$ is
equivalent to the typical best-of-M loss.

Method | MPJPE$\downarrow$ | PMPJPE$\downarrow$ | PCK$\uparrow$ | CPS$\uparrow$
---|---|---|---|---
w/o condition | 71.7 | 57.2 | 91.2 | 137.6
w/o $\mathcal{L}_{\text{HM}}$ | 72.4 | 56.2 | 92.2 | 157.3
w/o $\mathcal{L}_{\text{gen}}$ | 73.6 | 58.2 | 92.5 | 165.5
w/o $\mathcal{L}_{\text{MB}}$ | 76.0 | 58.5 | 91.6 | 161.4
$\mathcal{L}_{\text{MB}}(k=1)$ | 71.8 | 54.9 | 92.6 | 167.4
$\mathcal{L}_{\text{MB}}(k=10)$ | $\bm{70.7}$ | 54.9 | 93.3 | 168.3
$\mathcal{L}_{\text{MB}}(k=50)$ | 71.0 | 55.2 | 93.1 | 168.0
Ours (Full) | 71.0 | $\bm{54.2}$ | $\bm{93.4}$ | $\bm{171.0}$

Table 4: Ablation studies on the subset H36MA.

##  5 Conclusion

This paper presents a normalizing flow based method for the ambiguous inverse
problem of 3D human pose estimation from 2D inputs. We exploit the bijectivity
of the normalizing flow by utilizing the known 3D to 2D projection during
training. By incorporating uncertainty information from the heatmaps of a 2D
pose detector, valuable information is maintained which is discarded by
previous approaches. As demonstrated, the generated hypotheses reflect these
uncertainties and additionally show meaningful diversity along the ambiguous
depth of the joints. Furthermore, the introduction of a 3D pose discriminator
ensures the geometrical feasibility of the poses and a proposed generalization
of the best-of-M loss improves the performance. Experimental results show that
our method outperforms all previous multi-hypotheses approaches in most
metrics, especially on a challenging subset of Human3.6M containing highly
ambiguous examples.

#### Acknowledgements.

This work has been supported by the Federal Ministry of Education and Research
(BMBF), Germany, under the project LeibnizKILabor (grant no. 01DD20003), the
Center for Digital Innovations (ZDIN) and the Deutsche Forschungsgemeinschaft
(DFG) under Germany’s Excellence Strategy within the Cluster of Excellence
PhoenixD (EXC 2122).

## References

  * [1] Ijaz Akhter and Michael J. Black.  Pose-conditioned joint angle limits for 3d human pose reconstruction.  In The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2015.
  * [2] Mykhaylo Andriluka, Leonid Pishchulin, Peter Gehler, and Bernt Schiele.  2d human pose estimation: New benchmark and state of the art analysis.  In The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2014.
  * [3] Lynton Ardizzone, Jakob Kruse, Sebastian J. Wirkert, Daniel Rahner, Eric W. Pellegrini, Ralf S. Klessen, Lena Maier-Hein, Carsten Rother, and Ullrich Köthe.  Analyzing inverse problems with invertible neural networks.  International Conference on Learning Representations (ICLR), 2019\.
  * [4] Lynton Ardizzone, Carsten Lüth, Jakob Kruse, Carsten Rother, and Ullrich Köthe.  Guided image generation with conditional invertible neural networks.  arXiv preprint arXiv:1907.02392, 2019.
  * [5] Benjamin Biggs, Sébastien Ehrhadt, Hanbyul Joo, Benjamin Graham, Andrea Vedaldi, and David Novotny.  3d multi-bodies: Fitting sets of plausible 3d human models to ambiguous image data.  In Advances in Neural Information Processing Systems (NeurIPS), 2020\.
  * [6] Christopher M. Bishop.  Mixture density networks.  Technical report, Aston University, 1994.
  * [7] Federica Bogo, Angjoo Kanazawa, Christoph Lassner, Peter Gehler, Javier Romero, and Michael J. Black.  Keep it SMPL: Automatic estimation of 3D human pose and shape from a single image.  In European Conference on Computer Vision (ECCV), 2016.
  * [8] Ching-Hang Chen and Deva Ramanan.  3d human pose estimation = 2d pose estimation + matching.  In The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2017.
  * [9] Hai Ci, Chunyu Wang, Xiaoxuan Ma, and Yizhou Wang.  Optimizing network structure for 3d human pose estimation.  In Proceedings of the IEEE International Conference on Computer Vision (ICCV), 2019.
  * [10] Laurent Dinh, David Krueger, and Yoshua Bengio.  NICE: non-linear independent components estimation.  In International Conference on Learning Representations (ICLR), 2015\.
  * [11] Laurent Dinh, Jascha Sohl-Dickstein, and Samy Bengio.  Density estimation using real NVP.  In International Conference on Learning Representations (ICLR), 2017\.
  * [12] Haoshu Fang, Yuanlu Xu, Wenguan Wang, Xiaobai Liu, and Song-Chun Zhu.  Learning pose grammar to encode human body configuration for 3d pose estimation.  In Proceedings of the AAAI Conference on Artificial Intelligence, 2018.
  * [13] Mathieu Germain, Karol Gregor, Iain Murray, and Hugo Larochelle.  Made: Masked autoencoder for distribution estimation.  In Proceedings of the International Conference on Machine Learning (ICML), 2015.
  * [14] Arthur Gretton, Karsten M. Borgwardt, Malte J. Rasch, Bernhard Schölkopf, and Alexander Smola.  A kernel two-sample test.  Journal of Machine Learning Research (JMLR), 13(25), 2012.
  * [15] Ishaan Gulrajani, Faruk Ahmed, Martin Arjovsky, Vincent Dumoulin, and Aaron C Courville.  Improved training of wasserstein gans.  In Advances in Neural Information Processing Systems (NeurIPS), 2017\.
  * [16] Abner Guzmán-rivera, Dhruv Batra, and Pushmeet Kohli.  Multiple choice learning: Learning to produce multiple structured outputs.  In Advances in Neural Information Processing Systems (NeurIPS), 2012\.
  * [17] Ikhsanul Habibie, Weipeng Xu, Dushyant Mehta, Gerard Pons-Moll, and Christian Theobalt.  In the wild human pose estimation using explicit 2d features and intermediate 3d representations.  In The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2019.
  * [18] Mir Rayat Imtiaz Hossain and James J. Little.  Exploiting temporal information for 3d pose estimation.  European Conference on Computer Vision (ECCV), 2018.
  * [19] Catalin Ionescu, Dragos Papava, Vlad Olaru, and Cristian Sminchisescu.  Human3.6m: Large scale datasets and predictive methods for 3d human sensing in natural environments.  IEEE Transactions on Pattern Analysis and Machine Intelligence (PAMI), 36(7), 2014.
  * [20] Ehsan Jahangiri and Alan L. Yuille.  Generating multiple diverse hypotheses for human 3d pose consistent with 2d joint detections.  International Conference on Computer Vision Workshops (ICCVW), 2017\.
  * [21] Angjoo Kanazawa, Michael J. Black, David W. Jacobs, and Jitendra Malik.  End-to-end recovery of human shape and pose.  In The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2018.
  * [22] Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  In International Conference on Learning Representations (ICLR), 2015\.
  * [23] Durk P Kingma, Tim Salimans, Rafal Jozefowicz, Xi Chen, Ilya Sutskever, and Max Welling.  Improved variational inference with inverse autoregressive flow.  In Advances in Neural Information Processing Systems (NeurIPS), 2016\.
  * [24] Ivan Kobyzev, Simon Prince, and Marcus Brubaker.  Normalizing flows: An introduction and review of current methods.  TPAMI, 2020.
  * [25] Muhammed Kocabas, Nikos Athanasiou, and Michael J. Black.  Vibe: Video inference for human body pose and shape estimation.  In The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), June 2020.
  * [26] Nikos Kolotouros, Georgios Pavlakos, Michael J. Black, and Kostas Daniilidis.  Learning to reconstruct 3d human pose and shape via model-fitting in the loop.  In Proceedings of the IEEE International Conference on Computer Vision (ICCV), 2019.
  * [27] Mun Wai Lee and Isaac Cohen.  Proposal maps driven mcmc for estimating human body pose in static images.  In The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2004.
  * [28] Chen Li and Gim Hee Lee.  Generating multiple hypotheses for 3d human pose estimation with mixture density network.  In The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2019.
  * [29] Chen Li and Gim Hee Lee.  Weakly supervised generative network for multiple 3d human pose hypotheses.  British Machine Vision Conference (BMVC), 2020.
  * [30] Jiefeng Li, Chao Xu, Zhicun Chen, Siyuan Bian, Lixin Yang, and Cewu Lu.  Hybrik: A hybrid analytical-neural inverse kinematics solution for 3d human pose and shape estimation.  In CVPR, 2021.
  * [31] Shichao Li, Lei Ke, Kevin Pratama, Yu-Wing Tai, Chi-Keung Tang, and Kwang-Ting Cheng.  Cascaded deep monocular 3d human pose estimation with evolutionary training data.  In The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2020.
  * [32] Julieta Martinez, Rayat Hossain, Javier Romero, and James J. Little.  A simple yet effective baseline for 3d human pose estimation.  In Proceedings of the IEEE International Conference on Computer Vision (ICCV), 2017.
  * [33] Dushyant Mehta, Helge Rhodin, Dan Casas, Pascal Fua, Oleksandr Sotnychenko, Weipeng Xu, and Christian Theobalt.  Monocular 3d human pose estimation in the wild using improved cnn supervision.  In International Conference on 3D Vision (3DV), 2017.
  * [34] Francesc Moreno-Noguer.  3d human pose estimation from a single image via distance matrix regression.  In The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2017.
  * [35] Tuomas P. Oikarinen, Daniel C. Hannah, and Sohrob Kazerounian.  Graphmdn: Leveraging graph structure and deep learning to solve inverse problems.  arXiv preprint arXiv:2010.13668, 2020.
  * [36] George Papamakarios, Theo Pavlakou, and Iain Murray.  Masked autoregressive flow for density estimation.  In Advances in Neural Information Processing Systems (NeurIPS), 2017\.
  * [37] Georgios Pavlakos, Luyang Zhu, Xiaowei Zhou, and Kostas Daniilidis.  Vibe: Video inference for human body pose and shape estimation.  In The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), June 2018.
  * [38] Danilo Rezende and Shakir Mohamed.  Variational inference with normalizing flows.  Proceedings of Machine Learning Research (PMLR), 2015.
  * [39] Marco Rudolph, Bastian Wandt, and Bodo Rosenhahn.  Same same but differnet: Semi-supervised defect detection with normalizing flows.  In Winter Conference on Applications of Computer Vision (WACV), 2021\.
  * [40] Saurabh Sharma, Pavan Teja Varigonda, Prashast Bindal, Abhishek Sharma, and Arjun Jain.  Monocular 3d human pose estimation by generation and ordinal ranking.  In Proceedings of the IEEE International Conference on Computer Vision (ICCV), 2019.
  * [41] Soshi Shimada, Vladislav Golyanik, Weipeng Xu, and Christian Theobalt.  Physcap: Physically plausible monocular 3d motion capture in real time.  ACM Transactions on Graphics, 39(6), 2020.
  * [42] E. Simo-Serra, A. Ramisa, G. Alenyà, C. Torras, and F. Moreno-Noguer.  Single image 3d human pose estimation from noisy observations.  In The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2012.
  * [43] Cristian Sminchisescu and Bill Triggs.  Covariance scaled sampling for monocular 3d body tracking.  In The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2001.
  * [44] Cristian Sminchisescu and Bill Triggs.  Kinematic jump processes for monocular 3d human tracking.  In The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2003.
  * [45] Ke Sun, Bin Xiao, Dong Liu, and Jingdong Wang.  Deep high-resolution representation learning for human pose estimation.  In The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2019.
  * [46] Esteban G. Tabak and Cristina V. Turner.  A family of nonparametric density estimation algorithms.  Communications on Pure and Applied Mathematics, 66(2), 2013.
  * [47] Esteban G. Tabak and Eric Vanden-Eijnden.  Density estimation by dual ascent of the log-likelihood.  Communications in Mathematical Sciences, 8(1), 2010.
  * [48] Bastian Wandt, Hanno Ackermann, and Bodo Rosenhahn.  A kinematic chain space for monocular motion capture.  In European Conference on Computer Vision Workshops (ECCVW), 2018\.
  * [49] Bastian Wandt and Bodo Rosenhahn.  Repnet: Weakly supervised training of an adversarial reprojection network for 3d human pose estimation.  In The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2019.
  * [50] Bastian Wandt, Marco Rudolph, Petrissa Zell, Helge Rhodin, and Bodo Rosenhahn.  Canonpose: Self-supervised monocular 3d human pose estimation in the wild.  In The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2021.
  * [51] Jue Wang, Shaoli Huang, Xinchao Wang, and Dacheng Tao.  Not all parts are created equal: 3d pose estimation by modelling bi-directional dependencies of body parts.  Proceedings of the IEEE International Conference on Computer Vision (ICCV), 2019.
  * [52] Christina Winkler, Daniel Worrall, Emiel Hoogeboom, and Max Welling.  Learning likelihoods with conditional normalizing flows.  arXiv preprint arXiv:1912.00042, 2019.
  * [53] Hongyi Xu, Eduard Gabriel Bazavan, Andrei Zanfir, William T Freeman, Rahul Sukthankar, and Cristian Sminchisescu.  Ghum & Ghuml: Generative 3d human shape and articulated pose models.  In The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2020.
  * [54] Jingwei Xu, Zhenbo Yu, Bingbing Ni, Jiancheng Yang, Xiaokang Yang, and Wenjun Zhang.  Deep kinematics analysis for monocular 3d human pose estimation.  In The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2020.
  * [55] Andrei Zanfir, Eduard Gabriel Bazavan, Hongyi Xu, Bill Freeman, Rahul Sukthankar, and Cristian Sminchisescu.  Weakly supervised 3d human pose and shape reconstruction with normalizing flows.  European Conference on Computer Vision (ECCV), 2020.
  * [56] Long Zhao, Xi Peng, Yu Tian, Mubbasir Kapadia, and Dimitris N. Metaxas.  Semantic graph convolutional networks for 3d human pose regression.  In The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2019.

## Appendix

##  Appendix A Qualitative Evaluation

###  A.1 Condition Influence

To further show the influence of the heatmap condition and of the loss
$\mathcal{L}_{\text{HM}}$ that forces the network to reflect the 2D detector
uncertainty in the 3D hypotheses, we present several qualitative results in
Fig. 7. Evidently, incorporating the heatmap condition alone already leads to
meaningful diversity along the $x$\- and $y$\- directions. Additionally
optimizing $\mathcal{L}_{\text{HM}}$ further increases the meaningful
diversity of the pose hypotheses such that the uncertainties of the 2D
detector as well as the depth ambiguities are modeled best.

###  A.2 Competitor Comparison

In Fig. 8, we show additional qualitative results comparing our method with
the competing methods [28, 40]. As can be seen, our method achieves
significantly higher diversity mainly for occluded joints. The competing
methods are unable to effectively model occlusions and uncertain detections.
They only achieve significant diversity along the ambiguous depth of the
joints.

###  A.3 High Confidence Detections

If the 2D detector has a high degree of confidence for the 2D pose detection
in a given image, then low variance in the generated 3D hypotheses along the
$x$\- and $y$\- directions is expected. To validate this, we show qualitative
results for images from Human3.6M and MPI-INF-3DHP with low 2D detector
uncertainty in Fig. 9. The generated hypotheses are shown from two
perspectives such that diversity along the image and depth directions can be
seen. Evidently, the hypotheses vary only slightly along the image directions
and thus are all consistent with the input image. They show meaningful
diversity along the ambiguous depth of the joints.

##  Appendix B Captured 2D Detector Uncertainty

In the following, we want to further verify that fitting a Gaussian to the
heatmap can capture the uncertainty of the 2D detector well. Therefore, for
each joint in the test split of Human3.6M, we show the mean of the standard
deviations of the fitted Gaussian together with the 2D error in Fig. 6. As can
be seen, the variances of the Gaussians correlate with the 2D error and thus
are a good surrogate for the uncertainty of the 2D detector.

Figure 6: Computed for all joints in the test split of Human3.6M.

##  Appendix C Performance Lower Number of Samples

To assess the influence of the number of generated hypotheses and make our
approach better comparable to Li  _et al_. [28], we evaluate on Human3.6M
under Protocol 1 (MPJPE) and Protocol 2 (PMPJPE) for lower number of
hypotheses in two different settings. However, we want to emphasize that our
main goal is to model the full posterior distribution, which requires a larger
number of samples. Instead of sampling from their model, Li _et al_. [28] take
the means of the Gaussian kernels as pose predictions. Thus, for better
comparison, we emulate this by running K-Means on our $M=200$ generated
hypotheses. Additionally, we compare the performance when _sampling_ from [28]
in Table 5 (rows marked with *). We outperform them in almost every setting
and metric.

Method | Hypo. |  MPJPE$\downarrow$ |  PMPJPE$\downarrow$ | Method | Hypo. |  MPJPE$\downarrow$ |  PMPJPE$\downarrow$
---|---|---|---|---|---|---|---
Li [28] | 5 | 52.7 | 42.6 | Ours (K-Means) | 5 | ${53.2}$ | ${38.4}$
*Li [28] | 5 | 74.9 | 63.3 | *Ours | 5 | ${59.2}$ | ${42.3}$
Li [28] | 10 | 70.3 | 59.7 | *Ours | 10 | ${55.0}$ | ${39.6}$
Li [28] | 200 | 59.6 | 50.2 | *Ours | 200 | ${44.3}$ | ${32.4}$

Table 5: Results on Human3.6M under Protocol 1 (MPJPE) and Protocol 2
(PMPJPE). The scores for the rows marked with * are computed by sampling from
the models.

##  Appendix D Inference Time

For inference time measurements, we run the code with PyTorch 1.7.1 on a
NVIDIA GeForce RTX 3090 (CUDA 11.4). The majority of the inference time comes
from the 2D detector ($32$ ms) and the Gaussian fitting process (70 ms). Due
to batch processing, generating multiple hypotheses brings nearly no overhead,
with an inference time of $4.6$ ms for a single and $5.1$ ms for 1000 samples.

Figure 7: Qualitative results of our full model, model without
$\mathcal{L}_{\text{HM}}$, and the model without condition. For visualization
purposes, more than three hypotheses are shown only for the most ambiguous
joint. Figure 8: Comparison with competing methods [28, 40]. For visualization
purposes, more than three hypotheses are shown only for the most ambiguous
joint. The model from Li _et al_. [28] can only generate five pose hypotheses.
Figure 9: Qualitative results for images from Human3.6M and MPI-INF-3DHP with
low 2D detector uncertainty. For each image, 50 pose hypotheses are generated
and shown from two perspectives.