# A test for fractal boundaries based on the basin entropy

Andreu Puy  Nonlinear Dynamics, Chaos and Complex Systems Group, Departamento
de Física, Universidad Rey Juan Carlos
Tulipán s/n, 28933 Móstoles, Madrid, Spain  Alvar Daza  Nonlinear Dynamics,
Chaos and Complex Systems Group, Departamento de Física, Universidad Rey Juan
Carlos
Tulipán s/n, 28933 Móstoles, Madrid, Spain  Department of Physics, Harvard
University, Cambridge, MA 02138, U.S.A.  Alexandre Wagemakers  Nonlinear
Dynamics, Chaos and Complex Systems Group, Departamento de Física, Universidad
Rey Juan Carlos
Tulipán s/n, 28933 Móstoles, Madrid, Spain  Miguel A.F. Sanjuán  Nonlinear
Dynamics, Chaos and Complex Systems Group, Departamento de Física, Universidad
Rey Juan Carlos
Tulipán s/n, 28933 Móstoles, Madrid, Spain  Department of Applied Informatics,
Kaunas University of Technology, Studentu 50-415, Kaunas LT-51368, Lithuania

###### Abstract

In dynamical systems, basins of attraction connect a given set of initial
conditions in phase space to their asymptotic states. The basin entropy and
related tools quantify the unpredictability in the final state of a system
when there is an initial perturbation or uncertainty in the initial state.
Based on the basin entropy, the $\ln 2$ criterion allows for efficient testing
of fractal basin boundaries at a fixed resolution. Here, we extend this
criterion into a new test with improved sensitivity that we call the $S_{bb}$
fractality test. Using the same single scale information, the $S_{bb}$
fractality test allows for the detection of fractal boundaries in many more
cases than the $\ln 2$ criterion. The new test is illustrated with the
paradigmatic driven Duffing oscillator, and the results are compared with the
classical approach given by the uncertainty exponent. We believe that this
work can prove particularly useful to study both high-dimensional systems and
experimental basins of attraction.

###### keywords:

basin entropy , fractals , boundaries , basins of attraction , uncertainty
exponent

###### PACS:

05.45.-a , 05.45.Df

††journal: Commun. Nonlinear Sci. Numer. Simul.

##  1 Introduction

Dynamical systems often show multistability with the coexistence of several
asymptotic states, known as _attractors_ [1]. In dissipative systems, the set
of initial conditions that asymptotically approach an attractor is called the
_basin of attraction_ [2]. In open Hamiltonian systems, instead of attractors
we have exits, and consequently the initial conditions leaving the system by
these exits are their _escape basins_ [3]. In both cases, basin boundaries can
either be smooth or fractal curves. Under repeated enlargement, fractal
boundaries reveal new structures at arbitrarily small scales. This leads to
non-integer dimensions and it is often considered as one of the hallmarks of
chaos [4, 5]. The existence of fractal basin boundaries is profoundly
intertwined with the unpredictability under uncertainty in the initial state
of the orbit’s attractor; this is the facet of chaos we aim to study here.

The classical method for studying the lack of predictability of a multistable
system is via the _uncertainty exponent_ $\alpha$ [6]. Basically, it measures
the fractal dimension of the boundaries counting the number of boxes that lie
between basins at different scales. The uncertainty exponent ranges from zero
for the most unpredictable basins to one for smooth boundaries. Nonetheless,
the uncertainty exponent presents some unavoidable numerical difficulties,
such as accessing arbitrarily small scales or exploring multiple box sizes.
Besides, it makes a poor use of the information obtained by sampling the phase
space with the boxes, since it only classifies them as certain (lying in the
interior) or uncertain (lying on a boundary).

Other tools able to quantify the attractors unpredictability given a finite
uncertainty in the initial states are the _basin entropy_ and the _boundary
basin entropy_ [7]. Starting with a similar tiling of phase space, the idea is
to compute the probabilities of going to each attractor within a box and
exploit Shannon’s information entropy. Since their formulation, the basin
entropy and the boundary basin entropy have been applied to experiments with
cold atoms [8], chaotic scattering [9, 10], biological systems [11, 12],
electronic micro/nanodevices [13], oscillators [14] and astrophysical models
[15, 16, 17], among others.

Based on the boundary basin entropy, the _$\ln 2$ criterion_ [7] to detect
fractal boundaries at a fixed resolution was developed. The use of a single
scale to test for fractality makes it both computationally and experimentally
convenient. The $\ln 2$ criterion applies to boundaries separating more than
two different basins, like the Wada basins [18]. The goal of this paper is to
extend the $\ln 2$ criterion, presenting a test as a criterion to ascertain
for fractal structures with an improved sensitivity that also applies for
boundaries separating only two basins.

The article is organized as follows: we start in section 2 with a quick
revision of the basin entropy and its relation to fractal structures on phase
space. We continue in section 3 introducing the fractal test and developing it
for the case when basins are known exactly in two-dimensional phase spaces.
Then, in sec. 4, we study finite resolution effects in discrete phase spaces
given by finite grids. In sec. 5, we give a recipe of the test and illustrate
it with an example. Moreover, in sec. 6, we extend the test to phase spaces in
any dimension. Finally, we summarize and discuss the results in sec. 7.

##  2 Basin entropy and the definition of fractality

Since our main goal is the identification of fractal boundaries, we first need
to consider the definition of _fractal_. Even though fractals are not defined
in the literature in a precise and unambiguous manner, they typically share
some of the following properties [19]: (1) fine or detailed structure at
arbitrarily small scales, (2) local and global irregularity non-describable by
ordinary geometry, (3) some notion of self-similarity, (4) a ‘fractal
dimension’ greater than the topological dimension and (5) simple and perhaps
recursive definitions. Here, we refer to fractal structures as those
fulfilling at least properties (1), (2) and (4). Nonetheless, numerically we
are always able to reach only up to a given scale. This connects with the
experience that physical systems do not exhibit ‘true’ fractal structures, in
the sense that there are only some finite accessible or even defined scales
[20].

The basin entropy provides a way to characterize phase space structures at a
given scale. We consider a region of the phase space with $M$ distinct basins
and we sample it using $N$ boxes of size $\varepsilon$ (see Fig. 1). The box
size $\varepsilon$ dictates the scale to study the boundary structures. These
boxes can be thought as initial states either with uncertainty $\varepsilon$
or an initial perturbation $\varepsilon$. We can construct the discrete
probability distribution of having a basin $k$ for each box, $p(k)$, as the
ratio of the volume of the basin $k$ in the box, $v(k)$, to the phase space
volume comprised by the box, $V$:

$p(k)\equiv\frac{v(k)}{V},$ (1)

such that $\sum_{k=1}^{M}p(k)=1$. Its associated information entropy $s$,

$s\equiv s(p(1),\dots,p(M))=\sum_{k=1}^{M}-p(k)\ln{p(k)},$ (2)

measures the final state unpredictability of an initial condition chosen at
random within the box defined by $\varepsilon$. It increases monotonically
with the number of basins within the box and tends to a maximum value $s=\ln
M$ as the probabilities $p(k)$ tend to the equiprobable conditions
$p(k)=\frac{1}{M}$ for all $k$.

(a)  (b)  (c)

Figure 1: To compute the basin entropy, we draw boxes of size $\varepsilon$
from the phase space. For each box, we calculate the information entropy from
the fractions of each basin within the box; then, averaging over all boxes we
calculate the basin entropy $S_{b}$; and averaging only over boundary boxes
(red dashed circles), the boundary basin entropy $S_{bb}$. These basins are
from the periodically driven Duffing oscillator
$\ddot{x}+0.15\dot{x}-x+x^{3}=F\sin\omega t$, and (a) $F=0.100$,
$\omega=0.200$, (b) $F=0.395$, $\omega=1.617$ and (c) $F=0.128$,
$\omega=1.106$. For $N_{b}=10^{4}$ disks boxes in the boundary of radius
$\varepsilon=0.025$, the figures are ordered from left to right by the
boundary basin entropy: $S_{bb}=0.465\pm 0.002<0.6323\pm 0.0011\leq\ln
2<0.760\pm 0.004$. Using the $\ln 2$ criterion, only basins (c) are tested
with fractal boundaries.

Then, the _basin entropy_ $S_{b}$ is defined as the average of the entropy of
the box $s$ for the total number of boxes $N$,

$S_{b}\equiv\frac{1}{N}\sum_{i=1}^{N}s(i),$ (3)

where $i$ labels the boxes. Therefore, for an initial random position in the
phase space region that is uncertain within a volume of size $\varepsilon$,
the basin entropy quantifies the unpredictability of the orbit’s attractor.
Moreover, the _boundary basin entropy_ $S_{bb}$ is defined as the average of
the entropies $s(i)$ restricted only to the $N_{b}$ boxes falling on basin
boundaries (in Fig. 1, these correspond to the red dashed boxes):

$S_{bb}\equiv\frac{1}{N_{b}}\sum_{i=1}^{N_{b}}s(i).$ (4)

The boundary basin entropy quantifies the unpredictability focusing only on
the unpredictable regions of the phase space: the basin boundaries.

Based on the boundary basin entropy, the $\ln 2$ criterion provides a
sufficient condition to test for fractal boundaries. It is based on the fact
that smooth boundaries separate only two basins, with the possible exception
of a countable number of points that can separate three or more basins at a
time. Therefore, for a sufficient large number of small boxes, for smooth
boundaries $S_{bb}\in[0,\ln 2\simeq 0.693]$. This is the case of the basins in
Fig. 1a, where for a box scale $\varepsilon=0.025$, $S_{bb}=0.465\pm 0.002$.
Equivalently, if $S_{bb}$ is significantly larger than $\ln 2$, this implies
the boundary will not be smooth (i.e., it will be fractal). This occurs for
the basins in Fig. 1c, where $S_{bb}=0.760\pm 0.004$. Nonetheless, the
criterion fails for the case of Fig. 1b with a manifested fractal boundaries
between two basins and $S_{bb}=0.6323\pm 0.0011$. Following this idea, there
is room to make a better test using a single scale with the basin entropy.

##  3  Beyond the ln 2 criterion: the $S_{bb}$ fractality test

Simple smooth boundaries are those which are locally flat (as a straight line
in 2D). According to our previous definition of a fractal, we consider as a
fractal boundary any boundary that is not simple smooth. A simple smooth
boundary could also include a few finite points where boundaries from several
basins intersect. Nevertheless, in the infinitely fine scale these regions
become negligible in front of the other locally flat regions.

Based on this idea and the boundary basin entropy, we can define a statistical
test to identify fractal structures, which we name $S_{bb}$ fractality test.
It consists on comparing the value of the boundary basin entropy $S_{bb}$ of
the boundary under study to the theoretical value $S_{bb}$ of a flat boundary,
at a given small box scale $\varepsilon$. If they are deemed statistically
significantly different, the boundary has a fractal structure. We will define
it formally after studying the boundary basin entropy of a flat boundary.
Since the boundary basin entropy is defined as an average of the box entropies
$s$ and different configuration of the boxes’ entropies can give rise to the
same value, this test is a sufficient but not necessary condition.
Nonetheless, it is more restrictive than the $\ln 2$ criterion and
consequently it can detect fractal boundaries in many more cases, e.g. in
regions with only two basins.

Assuming a perfect knowledge of the phase space structure for a flat boundary
in a two-dimensional phase space, we can derive the theoretical value of the
boundary basin entropy. We start by choosing a disk as box shape, because its
rotational symmetry under any angle allows to compute the boundary basin
entropy independently of the box orientation. Now, having disk boxes in a
phase space with a flat boundary, we can obtain the boundary basin entropy by
sliding the box from side to side along the boundary. If a disk box has a
radius $\varepsilon$ and is centered at the coordinate
$x_{0}\in[-\varepsilon,\varepsilon]$ from the perpendicular direction to the
boundary with origin on the boundary, the probability $p(x_{0})$ to have a
disk point in one of the basins is the fraction of the disk given by the
circular segment

$p(x_{0})=\frac{1}{2}+\frac{1}{\pi}\left[\frac{x_{0}}{\varepsilon}\sqrt{1-\left(\frac{x_{0}}{\varepsilon}\right)^{2}}+\arcsin{\frac{x_{0}}{\varepsilon}}\right].$
(5)

Then, $S_{bb}$ is given numerically by:

$S_{bb}=\frac{1}{2\varepsilon}\int_{-\varepsilon}^{\varepsilon}dx_{0}\
s\left(p\left(x_{0}\right),1-p\left(x_{0}\right)\right)=0.4395093(6).$ (6)

This means that if we were able to compute exactly the boundary basin entropy
of a smooth boundary using small disks we would get that result, and any other
number would correspond to a fractal case. However, in practice we always have
some errors induced by the use of a finite number of trajectories. The effects
that this can introduce into our test are explored in the next section.

##  4 Effects of finite grids

Both numerically and in experiments, basins are commonly calculated using a
finite number of grid points. The number of grid points within a box
determines the observed probabilities $\hat{p}(k)$ of the box. Ultimately,
this lack of detail results in an observed boundary basin entropy
$\hat{S}_{bb}$ with a systematic error or bias $\delta_{\hat{S}_{bb}}$. In
addition, we are generally limited to draw a finite number of boxes too. This
results in an observed boundary basin entropy with a statistical error
$\sigma_{\hat{S}_{bb}}$. Here, we study these effects for a two-dimensional
square grid with a flat boundary and disk boxes, which provides confidence
intervals for our fractal basin boundary test.

The systematic error of the observed boundary basin entropy $\hat{S}_{bb}$
depends on two factors: the extension of the grid and the angular orientation
of the grid. The effect of the grid extension decreases with the number of
grid points per axis. Our numerical simulations showed that this factor is no
longer relevant for larger values than $50\varepsilon_{g}$ grid points per
axis, where $\varepsilon_{g}$ is the disk radius in grid units. The second
factor, the angle between the grid and the boundary, is inherently
unavoidable. We have investigated this effect in Fig. 2, representing in
function of the disk radius $\varepsilon_{g}$ in grid units and for largely
extended grids, $\hat{S}_{bb}$ for different angles (legend) and the exact
$S_{bb}$ (magenta dashed line). The figure inset displays the absolute errors
of $\hat{S}_{bb}$. Even in the worst case scenario (red line), the systematic
error decreases approximately inversely proportional to the disk boxes radius
$\varepsilon_{g}$. Indeed, a power law fit gives an upper bound for the
systematic error:

$\delta_{\hat{S}_{bb}}^{\text{UB}}\simeq A\varepsilon_{g}^{B},$ (7)

with $A=0.224\pm 0.010$ and $B=-1.006\pm 0.014$.

Figure 2: For a grid with a flat boundary and disk boxes of radius
$\varepsilon_{g}$ in grid units, the systematic error of the observed boundary
basin entropy $\hat{S}_{bb}$ depends on the angle between the grid and the
boundary (legend): (outset) $\hat{S}_{bb}$ in function of $\varepsilon_{g}$
and the exact $S_{bb}$ value (magenta dashed line); (inset) absolute errors of
$\hat{S}_{bb}$ in function of $\varepsilon_{g}$. We considered large grid
extensions of around $128\varepsilon_{g}$ points per axis and $N_{b}=10^{5}$
disk boxes in the boundary. The worst systematic error is given by the angle
$0$ in red. In the infinite grid resolution ($\varepsilon_{g}\to\infty$) and
for any angle $\alpha$, we recover the exact boundary basin entropy $S_{bb}$.

On the other hand, there is a statistical error of the observed boundary basin
entropy $\hat{S}_{bb}$, which is due to the finite number $N_{b}$ of boxes in
the boundary. Indeed, $\hat{S}_{bb}$ is the average of the boxes’ entropies in
the boundary. This means that, by the central limit theorem, the statistical
error follows a Gaussian distribution with a standard deviation given by:

$\sigma_{\hat{S}_{bb}}=aN_{b}^{-\frac{1}{2}},$ (8)

where
$a\equiv\sqrt{\frac{1}{N_{b}-1}\sum_{i=1}^{N_{b}}\left(\hat{S}_{bb}-\hat{s}(i)\right)^{2}}$
is the sampling standard deviation and $\hat{s}(i)$ the observed entropy for
the box $i$.

Taking into account all the finite grid effects for 2D phase spaces, we can
formulate our test for fractal boundaries as follows: under an infinitesimal
disk box with radius $\varepsilon_{g}$ in grid units and a finite grid that is
largely extended (with at least $50\varepsilon_{g}$ grid points per axis), if
the observed boundary basin entropy $\hat{S}_{bb}$ with a standard deviation
$\sigma_{\hat{S}_{bb}}$ (Eq. 8) is deemed statistically significant away from
a flat boundary, either below the exact boundary basin entropy value $S_{bb}$
(Eq. 6) or above the upper-bound systematic error value
$S_{bb}+\delta_{\hat{S}_{bb}}^{\text{UB}}$
($\delta_{\hat{S}_{bb}}^{\text{UB}}$ is given by Eq. 7), the boundary has a
fractal structure. Using one standard deviation $\sigma_{\hat{S}_{bb}}$ for
the statistical error, we can express the $S_{bb}$ fractality test by the
following sufficient conditions:

$\displaystyle\hat{S}_{bb}$ $\displaystyle<S_{bb}-\sigma_{\hat{S}_{bb}},$ (9)
$\displaystyle\hat{S}_{bb}$
$\displaystyle>S_{bb}+\delta_{\hat{S}_{bb}}^{\text{UB}}+\sigma_{\hat{S}_{bb}}.$
(10)

##  5 Example of application

The periodically driven Duffing oscillator is a paradigmatic model that
accounts for nonlinear elastic effects in large displacements of a forced
damped elastic structure. It is defined by:

$\ddot{x}+\gamma\dot{x}-x+x^{3}=F\sin\omega t,$ (11)

where $x$ is the displacement of the oscillator at time $t$, $\gamma$ is the
damping coefficient, $F$ is the forcing amplitude and $\omega$ is the
frequency of the driving. Depending on these parameters, the system exhibits a
wide variety of dynamics. Here, we investigate a parameter space region given
by $\gamma=0.15$, $F\in[0.1,0.5]$ and $\omega\in[0.2,2.5]$. We search for
fractal boundaries in the basins of attraction given by a finite grid in the
phase space region
$\Omega=\Omega_{x}\times\Omega_{\dot{x}}=[-2.5,2.5]\times[-2.5,2.5]$, with
$10^{3}$ points per axis.

On the following, we use the example to illustrate an easy-to-follow recipe
for the $S_{bb}$ fractality test for finite grids in two-dimensional phase
spaces:

  1. 1.

Choose the box disk radius $\varepsilon$ appropriately. On the one hand, the
smaller the box, the finer scales we can study. On the other hand, we want to
have sufficient trajectories per box to get good estimates for the
probabilities. For our example, we found that a value of $\varepsilon_{g}=5$ /
$\varepsilon=0.025$ gave good results.

  2. 2.

Verify that we are in the largely extended grid limit, with at least
$50\varepsilon_{g}$ grid points per axis. In our example, we had
$10^{3}=200\varepsilon_{g}$ points per axis.

  3. 3.

Draw $N_{b}$ disk boxes in the boundary, uniformly at random from the phase
space (and not only at grid points). In our example, we took $N_{b}=10^{4}$
boxes in the boundary.

  4. 4.

For all boxes in the boundary, calculate both the observed probabilities of
the basins and the observed box’s entropy.

  5. 5.

Compute the observed boundary basin entropy $\hat{S}_{bb}$ and its standard
deviation $\sigma_{\hat{S}_{bb}}$ (Eq. 8).

  6. 6.

Compute the upper-bound systematic error $\delta_{\hat{S}_{bb}}^{\text{UB}}$
(Eq. 7) for the current value of $\varepsilon_{g}$. In our example, we had
$\delta_{\hat{S}_{bb}}^{\text{UB}}=0.0474$.

  7. 7.

Check if the computed value of $\hat{S}_{bb}$ lies in the interval provided by
the systematic and statistical errors. If the value is outside the interval,
we can affirm that according to our test the boundary is fractal. Otherwise,
it is most likely to be a smooth boundary, although there could be
pathological cases leading to wrong results. It is important to recall that
this fractality test is a sufficient but not necessary condition.

We display the results in the parameter space of Fig. 3a. The color of each
parameter value corresponds to its boundary basin entropy $\hat{S}_{bb}$
value: white is for regions with a single basin, cold colors are for regions
compatible with simple smooth boundaries and hot colors are for fractal basin
boundaries. In particular, red colors are for fractal basin boundaries
detected by the $\ln 2$ criterion. We can observe that this just accounts for
a $27\%$ fraction of all the fractal basin boundaries detected by the single-
scale fractal basin boundary test.

(a)  (b)

Figure 3: For the periodically driven Duffing oscillator
$\ddot{x}+0.15\dot{x}-x+x^{3}=F\sin\omega t$ in the parameter space
$(F,\omega)$, we obtain similar regions with fractal boundaries (hot colors)
and smooth boundaries (cold colors), both for (a) the $S_{bb}$ fractality test
and (b) the uncertainty exponent $\alpha$. In white, we have regions with a
single basin. In (a), the color of each point corresponds to the boundary
basin entropy $\hat{S}_{bb}$, for disk boxes with radius $\varepsilon=0.025$ /
$\varepsilon_{g}=5$ and $N_{b}=10^{4}$ boxes in the boundary: cold colors are
compatible with simple smooth boundaries and hot colors for fractal basin
boundaries (red colors are detected by the $\ln 2$ criterion). In (b), we
represent the uncertainty exponent with a transition point at $\alpha=0.8$
from fractal to smooth boundaries; this value is arbitrary and based on the
observed basins.

To evaluate the performance of the method, we compare it to the uncertainty
exponent $\alpha$ (Fig. 3b). From the observed basins, we have chosen the
arbitrary value of $\alpha=0.8$ as the transition point between smooth and
fractal boundaries: $\alpha\geq 0.8$ correspond to smooth boundaries (cold
colors) and $\alpha<0.8$ to fractal basin boundaries (hot colors). Again,
parameters with a single basin are represented in white. Indeed, both the
uncertainty exponent and the $S_{bb}$ fractality test give qualitatively
similar regions with fractal boundaries. However, there are some quantitative
differences due to numerical errors in both methods.

##  6 Extension to any phase space dimension

We can generalize the test to phase spaces for any dimension. The test again
compares for small box scales $\varepsilon$ the boundary basin entropy
$S_{bb}$ of the boundary under study to the theoretical value of that of a
flat boundary of the corresponding dimension. Here, we only consider the case
with an exact knowledge of the phase space. We believe that the study for
basins with finite grids can be developed similarly to the two-dimensional
case.

For a given dimension, we can derive the theoretical value of the boundary
basin entropy of a phase space with a flat boundary. The box is now an
hyperball of radius $\varepsilon$ centered at a coordinate $x_{0}$ in the
direction that indicates the distance to the flat boundary. Furthermore, the
probability $p(x_{0})$ to have a hyperball point in one basin is the fraction
of the hyperball given by the hyperspherical cap [21]:

$p(x_{0})=-\frac{1}{2}\text{sgn}\left(\frac{x_{0}}{\varepsilon}\right)I_{1-\left(\frac{x_{0}}{\varepsilon}\right)^{2}}\left(\frac{D+1}{2},\frac{1}{2}\right)+\Theta\left(\frac{x_{0}}{\varepsilon}\right),$
(12)

where $D$ is the dimension of the phase space, $\text{sgn}(x)$ is the sign
function, $I_{x}\left(a,b\right)$ is the regularized incomplete beta function
and $\Theta(x)$ is the Heaviside function. We have this quantity plotted for
several dimensions (see colorbar) in Fig. 4a. For $D=1$, the probability
$p(x_{0})$ is linear; for larger values of $D$, it deviates from this
behavior; and in the limiting case $D\to\infty$, $p(x_{0})$ becomes dominated
by the $\Theta\left(\frac{x_{0}}{\varepsilon}\right)$ term and has a switching
behavior. This large dimension behavior is understood because most of the
volume of a high-dimensional hyperball lies within two parallel hyperplanes at
a small distance from its center of order
$\mathcal{O}\left(\frac{\varepsilon}{\sqrt{D-1}}\right)$ [22].

(a)  (b)

Figure 4: For a hyperball box and a flat boundary in the phase space, we plot
(a) for several dimensions of the phase space (colorbar), the probability
$p(x_{0})$ to have a hyperball point in one basin in function of the box
center coordinate $x_{0}$, which indicates the distance to the flat boundary;
and (b) the boundary basin entropy $S_{bb}$ relationship with the dimension
$D$ of the phase space (magenta circles). $S_{bb}$ decreases potentially fast
with a power fit $S_{bb}=AD^{B}$ (purple line): $A=0.898\pm 0.006$ and
$B=-0.4995\pm 0.0012$.

Moreover, we calculate the boundary basin entropy $S_{bb}$ analogously to the
two-dimensional case, following Eq. 6. We have the $S_{bb}$ relationship with
the dimension $D$ in Fig. 4b; we can see how the boundary basin entropy goes
to zero in the infinite dimension limit following the expression $S_{bb}\sim
D^{-\frac{1}{2}}$. We have tabulated the boundary basin entropy $S_{bb}$ for
the first five $D$ dimensions in table 1.

Table 1: Boundary basin entropy $S_{bb}$ for hyperball boxes in a flat boundary of a $D$ dimensional phase space, for the first five $D$ dimensions. $D$ | $S_{bb}$
---|---
1 | 0.499999(9)
2 | 0.4395093(6)
3 | 0.39609176(4)
4 | 0.36319428(1)
5 | 0.33722572

##  7 Conclusions

In this paper, we have applied the theory of the basin entropy and related
tools to characterize fractal basin boundaries using a single scale. Our work
consists on comparing, under a given box scale, the value of the boundary
basin entropy $S_{bb}$ of the boundary to the $S_{bb}$ value of a smooth
boundary. In contrast to the former basin-entropy-based test, the $\ln 2$
criterion, it achieves both an improved sensitivity and a capacity to identify
fractal basin boundaries separating only two basins. We call it $S_{bb}$
fractality test. Nonetheless, this test is still a sufficient but not
necessary condition. One could think about an improvement by means of the
probability density function (PDF) of the probabilities in the boundary boxes.
The configuration of a PDF is more restrictive than an average (the boundary
basin entropy), but it is also more demanding to implement.

On the other hand, testing for fractal basin boundaries with this line of work
has important advantages compared to the classical approach given by the
uncertainty exponent $\alpha$. While both methods hold numerical limitations
in the infinitely fine scale, the uncertainty exponent requires accessing
multiple scales of the system. This is not only numerically inconvenient but,
for some physical systems, it may be impossible. Indeed, our methods only
require accessing a single scale and have a natural formulation for
experimental basins of attraction.

Finally, we believe the extension for basins in any phase space dimension
could prove useful for studying the unpredictability in higher-dimensional
systems, a generally unexplored area in the field; in particular, it could
complement current studies for network systems [23, 24, 25].

## Acknowledgements

This work has been supported by the Spanish State Research Agency (AEI) and
the European Regional Development Fund (ERDF, EU) under Projects No.
FIS2016-76883-P and No. PID2019-105554GB-I00.

## References

  * [1] Feudel U, Complex dynamics in multistable systems, Int J Bifurcat Chaos 2008;18:1607–1626. https://doi.org/10.1142/S0218127408021233
  * [2] Nusse HE, Yorke JA, Basins of attraction, Science 1996;271:1376–1380. https://doi.org/10.1126/science.271.5254.1376
  * [3] Aguirre J, Vallejo JC, Sanjuán MAF, Wada basins and chaotic invariant sets in the Hénon-Heiles system, Phys Rev E 2001;64:066208. https://doi.org/10.1103/PhysRevE.64.066208
  * [4] McDonald SW, Grebogi C, Ott E, Yorke JA, Fractal basin boundaries, Physica D 1985;17:125–153. https://doi.org/10.1016/0167-2789(85)90001-6
  * [5] Aguirre J, Viana RL, Sanjuán MAF, Fractal structures in nonlinear dynamics, Rev Mod Phys 2009;81:333–386. https://doi.org/10.1103/RevModPhys.81.333
  * [6] Grebogi C, McDonald SW, Ott E, Yorke JA, Final state sensitivity: an obstruction to predictability, Phys Lett A 1983;99:415–418. https://doi.org/10.1016/0375-9601(83)90945-3
  * [7] Daza A, Wagemakers A, Georgeot B, Guéry-Odelin D, Sanjuán MAF, Basin entropy: a new tool to analyze uncertainty in dynamical systems, Scientific Reports 2016;6:31416. https://doi.org/10.1038/srep31416
  * [8] Daza A, Georgeot B, Guéry-Odelin D, Wagemakers A, Sanjuán MAF, Chaotic dynamics and fractal structures in experiments with cold atoms, Phys Rev A 2017;95:013629. https://doi.org/10.1103/PhysRevA.95.013629
  * [9] Bernal JD, Seoane JM, Sanjuán MAF, Uncertainty dimension and basin entropy in relativistic chaotic scattering, Phys Rev E 2018;97:042214. https://doi.org/10.1103/PhysRevE.97.042214
  * [10] Nieto AR, Seoane JM, Alvarellos J, Sanjuán MAF, Resonant behavior and unpredictability in forced chaotic scattering, Phys Rev E 2018;98:062206. https://doi.org/10.1103/PhysRevE.98.062206
  * [11] Donepudi R, Ramaswamy R, The collective dynamics of nf- $\kappa$b in cellular ensembles, Eur Phys J-Spec Top 2018;227:851–863. https://doi.org/10.1140/epjst/e2018-800014-7
  * [12] Mugnaine M, Andrade F, Szezech Jr J, Bazeia D, Basin entropy behavior in a cyclic model of the rock-paper-scissors type, Europhys Lett 2019;125:58003.
  * [13] Gusso A, Viana RL, Mathias AC, Caldas IL, Nonlinear dynamics and chaos in micro/nanoelectromechanical beam resonators actuated by two-sided electrodes, Chaos, Soliton Fract 2019;122:6–16. https://doi.org/10.1016/j.chaos.2019.03.004
  * [14] Kong G, Zhang Y, A special type of explosion of basin boundary, Phys Lett A 2019;383:1151–1156. https://doi.org/10.1016/j.physleta.2019.01.017
  * [15] Zotos EE, Suraj MS, Basins of attraction of equilibrium points in the planar circular restricted five-body problem, Astrophys Space Sci 2018;363:20. https://doi.org/10.1007/s10509-017-3240-7
  * [16] Zotos EE, On the Newton–Raphson basins of convergence of the out-of-plane equilibrium points in the Copenhagen problem with oblate primaries, Int J Nonlin Mech 2018;103:93–103. https://doi.org/10.1016/j.ijnonlinmec.2018.05.002
  * [17] Zotos EE, Riaño-Doncel A, Dubeibe FL, Basins of convergence of equilibrium points in the generalized Hénon–Heiles system, Int J Nonlin Mech 2018;99:218–228. https://doi.org/10.1016/j.ijnonlinmec.2017.12.004
  * [18] Kennedy J, Yorke JA, Basins of Wada, Physica D 1991;51:213–225. https://doi.org/10.1016/0167-2789(91)90234-Z
  * [19] Falconer K, Fractal geometry: mathematical foundations and applications, John Wiley & Sons, 2004.
  * [20] Avnir D, Biham O, Lidar D, Malcai O, Is the geometry of nature fractal?, Science 1998;279:39–40. https://doi.org/10.1126/science.279.5347.39
  * [21] Li S, Concise formulas for the area and volume of a hyperspherical cap, Asian Journal of Mathematics and Statistics 2011;4:66–70. DOI: 10.3923/ajms.2011.66.70
  * [22] Hopcroft J, Kannan R, Computer science theory for the information age 2012.
  * [23] Menck PJ, Heitzig J, Marwan N, Kurths J, How basin stability complements the linear-stability paradigm, Nat Phys 2013;9:89–92. https://doi.org/10.1038/nphys2516
  * [24] Menck PJ, Heitzig J, Kurths J, Schellnhuber HJ, How dead ends undermine power grid stability, Nature Commun 2014;5:1–8. https://doi.org/10.1038/ncomms4969
  * [25] Rakshit S, Bera BK, Majhi S, Hens C, Ghosh D, Basin stability measure of different steady states in coupled oscillators, Scientific Reports 2017;7:45909. https://doi.org/10.1038/srep45909


Miguel<EMAIL_ADDRESS>Martin<EMAIL_ADDRESS>Michael<EMAIL_ADDRESS>Timothy
<EMAIL_ADDRESS>University of Edinburgh
Edinburgh, UK Monash University
Melbourne, Australia Vision-based system identification

# Vision-based system identification and 3D keypoint discovery using dynamics
constraints

###### Abstract

This paper introduces V-SysId, a novel method that enables simultaneous
keypoint discovery, 3D system identification, and extrinsic camera calibration
from an unlabeled video taken from a static camera, using only the family of
equations of motion of the object of interest as weak supervision. V-SysId
takes keypoint trajectory proposals and alternates between maximum likelihood
parameter estimation and extrinsic camera calibration, before applying a
suitable selection criterion to identify the track of interest. This is then
used to train a keypoint tracking model using supervised learning. Results on
a range of settings (robotics, physics, physiology) highlight the utility of
this approach.

##  1 Introduction

An understanding of the motion and physics of objects in the real world is a
hallmark of the human visual system. Humans have the ability to identify
objects and their properties (eg. mass, friction, elasticity) as they move and
interact in the world, due to our intuitive understanding of common
trajectories, object interactions, and outcomes. This ability is typically
studied under the umbrella of _intuitive physics_ [(5), (52), (20), (4)], and
often considered a critical component for machines to be able to think more
like humans. In the context of machine learning systems, this ability can be
distilled to a requirement for _unsupervised_ 3D object localization and
physical parameter estimation (also known as system identification) from a
sensory stream, subject to some inductive bias or intuitive physics prior.

Taking inspiration from this view, this paper introduces V-SysId, a novel
method that enables simultaneous keypoint discovery, 3D system identification,
and extrinsic camera calibration from a single unlabeled video taken from a
static camera, using only the family of equations of motion of the object of
interest as weak supervision. Crucially, our approach is able to identify the
correct object(s) in a scene even in the presence of other moving objects or
distractors. This property is key, as it greatly increases applicability to
cluttered real world environments. This allows us to perform queries such as “
_find the 3D location of the bouncing ball, and determine its restitution
coefficient_ ” (Fig. 1).

Figure 1: Given an unlabeled video containing moving objects and an equation
of motion, our V-SysId identifies the trajectory of the object of interest,
along with its physical parameters (e.g. restitution coefficient, initial
height), and 3D pose relative to the camera.

V-SysId follows a 3-stage process of keypoint track proposal, optimization,
and selection (Fig. 2). The optimisation process alternates between maximum
likelihood extrinsic camera calibration and maximum likelihood physical
parameter estimation for motion tracks detected in video. This joint
optimisation can be unstable, which we address through the inclusion of a
curriculum-based optimisation strategy, alongside a maximum entropy criterion
for keypoint identification. A key benefit of V-SysId is that a neural network
is _not_ needed for discovery or system identification in our pipeline. This
means that V-SysId enables keypoint discovery with high-resolution images; and
can also perform system identification in _single_ videos, without the need to
obtain large datasets, which is particularly useful in robotics applications,
where data collection for neural network training can be laborious and time-
consuming. The keypoints discovered by V-SysId can be used as pseudo-labels to
train a supervised keypoint detector, for downstream tracking or control.

These properties provide significant flexibility to V-SysId, enabling its use
in real world environments with important applications for control, physics
understanding, and health monitoring. Specifically, we show that the V-SysId
can be applied to end-effector localization and extrinsic camera calibration,
bouncing ball discovery and physical property estimation, and breathing
frequency estimation from chest videos - all unlabeled and without regions of
interest provided a priori. This is made possible by the fact that V-SysId
identifies keypoints belonging to objects of interest present in scenes, while
ignoring any other moving objects or artifacts that do not follow the expected
dynamical constraints. This alleviates the need for hand-crafted object
segmentation methods or tricks to selectively remove parts of the image that
may contain moving distractors; and allows keypoint discovery at a fraction of
the computational expense of unsupervised neural methods that learn to
identify and model every moving object in an image.

Figure 2: Our V-SysId comprises 3 stages. Stage 1 extracts keypoint tracks
from a video using a grid keypoint detector + KLT tracking. Each of these 2D
tracks is passed to Stage 2, where the physical parameters
$\theta=\\{\bm{\eta},\mathbf{p}_{0},\mathbf{v}_{0}\\}$ of the 3D equation of
motion $f$, and the camera pose parameters $R,\mathbf{t}$ are optimized in
order to minimize the difference between the projected 3D trajectory (black,
Stage 2) and the 2D keypoint track observed (red, Stage 2). Stage 3 chooses
the best trajectory and corresponding parameters as those which maximize the
sum of projected likelihood and a trajectory entropy criterion. Here, a
bouncing ball scene with 2 moving distractors is shown, where the bouncing
ball is correctly discovered as the object that corresponds to the highest
entropy motion that fits the equation of motion $f$.

##  2 Related Work

System identification and physics understanding are key to allow machine
learning agents to interact with the real world. System identification is
typically performed using proprioceptive trajectory data directly, and there
has been extensive research across a range of fields [(29), (7), (8), (57),
(56), (36)] in support of this. Recent contributions include developments in
physical parameter estimation [(6), (9)], simulator learning [(42), (47)],
simulation alignment for robot interaction [(2)], trajectory generation [(27)]
and compositionality [(1), (35)].

Unsupervised system identification from vision is a recent area of research
that removes the requirements for trajectory data, with approaches including
unsupervised physical parameter estimation [(24), (31), (40)], structured
latent space learning [(32), (19), (25)], and Hamiltonian/Lagrangian learning
[(18), (51), (58)]. Unfortunately, these approaches are still relatively
limited in the complexity of scene they can model, and typically restricted to
toy problems and simulated environments. In this work we aim to improve upon
(Jaques2020Physics-as-Inverse-Graphics:Video, ; Kandukuri2020LearningPhysics,
; Murthy2020GradSim:Control, )’s limitation to simulated environments by
performing physical parameter estimation on real dynamical scenes with
distractors.

The seminal GALILEO model Wu2015GalileoLearning  demonstrated physical system
identification and simulation alignment using the Physics101 dataset
Wu2016PhysicsVideos . A key shortcoming of Galileo is that it assumes that the
camera is parallel to the plane of motion, and relies on manually identified
object tracks to lift the visual scenes onto object positions. In contrast
V-SysId is able to simultaneously estimate 3D trajectories and camera pose
relative to the scene from arbitrary camera angles, greatly increasing its
applicability to real world scenes. Furthermore, V-SysId automatically
identifies object tracks from keypoint proposals without needing human
intervention, allowing us to automatically discover the objects of interest in
video that are governed by the relevant equations of motion.

Keypoint discovery Keypoints are a natural representation for object parts,
with keypoint detection and tracking one of the earliest and most studied
areas of computer vision. Approaches like SIFT Lowe2004 , FAST
Rosten2006MachineDetection  and ORB Rublee2011ORB:SURF  are still widely used
to perform SLAM, SFM, VO111Simultaneous Localisation and Mapping, Structure-
from-Motion, Visual Odometry. and other tracking tasks (using, e.g. a KLT
tracker Tomasi1991DetectionFeatures ). Given keypoint trajectories, the
problem of inferring the 3D structure of a 2D trajectory using assumptions
about the dynamics has been coined "trajectory triangulation" by
Avidan2000TrajectorySequence ; Kaminski2002GeneralTriangulation , who assume
that objects follow a straight-line or conic-section trajectory in 3D space,
and that physical parameters can be uniquely identified using multiple
cameras. In contrast, our method assumes only a single static monocular view.
Other approaches to infer moving object structure using motion constraints
include Fitzgibbon2000MultibodyObjects ; Han2003MultipleCameras ;
David2004SoftPOSIT:Determination ; Scaramuzza2009AbsoluteConstraints .

When it comes to 2D keypoint discovery, several recent works have proposed
neural network based methods that use a regularized reconstruction objective
to discover objects of interest in an image Jakab2018UnsupervisedGeneration ;
Jakab2019LearningTranslation ; Kulkarni2019UnsupervisedControl ;
Minderer2019UnsupervisedVideos ; Gopalakrishnan2020UsupervisedPredictability ;
Das2020Model-BasedDemonstrations , which can be used for downstream control
tasks. However, these approaches lack the ability to estimate keypoint depth,
limiting their application in realistic control scenarios. Even though these
approaches obtain semantically meaningful keypoints (and in some instances are
able to ignore scene objects with unpredictable motion
Gopalakrishnan2020UsupervisedPredictability ), they require visual inspection
in order to obtain interpretability. In contrast, V-SysId provides equation-
driven keypoint discovery, ensuring a known semantic meaning for learned
keypoints. A parallel stream of research tackles this from a geometric
perspective, where 3D keypoints are inferred using camera motion cues or
geometric constraints Suwajanakorn2018DiscoveryReasoning ;
Jau2020DeepConstraints ; Vijayanarasimhan2017SfM-Net:Video ;
Wei2020DeepSFM:Adjustment . Even though this approach has been used in complex
real world settings, these keypoints lack semantic meaning, making these
unsuitable for semantic discovery queries (eg. “find the bouncing ball
following these dynamics”).

The use of dynamics as a learning constraint has not been explored in keypoint
discovery literature to date. This work proposes a method to integrate
dynamical inductive biases into the keypoint discovery process, enabling
extrinsic camera calibration and physics-guided discovery of objects of
interest alongside the corresponding physical parameter estimation.

##  3 Method

Our goal is to discover the 3D trajectory of an object of interest in a video
with possibly many moving objects, given only its family of motion dynamics,
$f$. To this end, we must estimate: a) 2D keypoint locations $\mathbf{k}_{t}$
of the object of interest in each frame $\mathbf{I}_{t}$; b) physical
parameters and initial conditions $\bm{\theta}$, of the equation of motion
$f(\bm{\theta})$; and c) camera rotation and translation relative to the scene
$[R,\mathbf{t}]$. Joint estimation of these quantities would be intractable,
so we split the objective into tractable components. Our method, V-SysId, has
3 stages (Fig. 1). We first describe the physical parameter+camera pose
estimation stage.

###  3.1 Physical parameter and camera pose estimation

Setup Let us assume we have a set of N 2D keypoint tracks
$\mathbf{K}=\\{\tilde{\mathbf{k}}^{n}_{1:T}\\}_{n=1}^{N}$ across the video
$\mathbf{I}_{1:T}$, and a family of 3D equations of motion $f$ with unknown
physical parameters $\bm{\eta}$ and initial position and velocity
$\mathbf{p}_{0}$ and $\mathbf{v}_{0}$, respectively. The equation $f$ can be
rolled out over $T$ time steps using a standard integration method in order to
obtain a 3D trajectory $\mathbf{p}_{1:T}=f(\bm{\theta})$, where
$\bm{\theta}=\\{\bm{\eta},\mathbf{p}_{0},\mathbf{v}_{0}\\}$.

Objective Our goal is to maximize the likelihood of the observed keypoint
trajectory $\tilde{\mathbf{k}}_{1:T}$ w.r.t. the physical parameters and
initial conditions, $\bm{\theta}$, and the extrinsic camera rotation and
translation, $[R\;\mathbf{t}]$:

$\bm{\theta}^{*},R^{*},\mathbf{t}^{*}=\operatorname*{arg\,max}_{\bm{\theta},R,\mathbf{t}}\;p(\tilde{\mathbf{k}}_{1:T}|\bm{\theta},R,\mathbf{t}),$
(1)

where we factorize the trajectory likelihood as:

$p(\tilde{\mathbf{k}}_{1:T}|\bm{\theta},R,\mathbf{t})=\prod_{t}p(\tilde{\mathbf{k}}_{t}|\bm{\theta},R,\mathbf{t})=\prod_{t}\mathcal{N}(\tilde{\mathbf{k}}_{t}|\mathbf{k}_{t}(\bm{\theta},R,\mathbf{t}),\sigma^{2}),$

and $\mathbf{k}_{t}(\bm{\theta},R,\mathbf{t})$ are the 2D projection of the
simulated 3D trajectory (given by $f(\theta$)):

$\displaystyle\mathbf{k}_{t}(\bm{\theta},R,\mathbf{t})=[\tilde{\mathbf{p}}_{x,t}/\tilde{\mathbf{p}}_{z,t},\tilde{\mathbf{p}}_{y,t}/\tilde{\mathbf{p}}_{z,t}]\quad;\quad\tilde{\mathbf{p}}_{t}=\mathbf{M}\,[R\;\mathbf{t}]\,\mathbf{p}_{t}$
(2)

with $\mathbf{M}$ being the intrinsic camera matrix. In this work we assume
known camera intrinsics.

In order to reduce the space of possible solutions (and therefore local
minima) of Step 1 above, we restrict the camera rotation matrix $R$ to have
$\text{roll}=0$. This means the camera cannot rotate about its projection
axis, which is the case in the vast majority of settings. Using the projection
plane in camera coordinates as $xy$ and the projection axis as $z$, we
parametrize $R$ as
$R(\alpha,\beta)=\text{EulerRotationMatrix}(\alpha,\beta,0)$, where $\alpha$,
$\beta$ and $\gamma=0$ correspond to the pitch, yaw and roll Euler angles,
respectively. We found that this parametrization greatly improves results and
optimization stability.

Optimization To maximize (3.1) we apply an iterative optimization procedure.
Given an initial estimate for $\bm{\theta}$, $R$ and $\mathbf{t}$, we
alternate the following steps until convergence:

  1. 1.

Keeping $\bm{\theta}$ fixed, maximize (3.1) w.r.t. $R$ and $\mathbf{t}$ using
gradient descent;

  2. 2.

Keeping $R$ and $\mathbf{t}$ fixed, maximize (3.1) wrt $\bm{\theta}$ using
gradient descent (with numerical or analytical Belbute-Peres2018End-to-
EndControl  differentiation) or global optimizer (e.g. CEM
Rubinstein1997OptimizationEvents ; BO Mockus1989BayesianOptimization ).

Estimation of the physical parameters over the full sequence (possibly
hundreds of timesteps) is prone to local minima, as the dependency on the
parameters can be highly non-linear222Global optimizers have a slight
advantage in this case, although they require very many iterations to find a
good minimum.. This is further affected by the use of a non-optimized camera
pose at the first iteration. In order to address this, we start by performing
a step of physical parameter and pose estimation on a small initial trajectory
interval, $T_{0}$, adding $m$ points to the trajectory at each iteration, as
described in the appendix, Algorithm 1.

###  3.2 Trajectory proposal and Selection

Proposal In an unlabeled video, ground-truth 2D keypoints are not available,
but keypoint trajectories are required to maximize the likelihood in (3.1).
Joint estimation of physical parameters with a neural network-based keypoint
detector would be hard to optimize due to the difficulty of backpropagating
through physics rollouts and camera projection into a CNN in a stable manner
(Jaques2020Physics-as-Inverse-Graphics:Video ). Therefore, we propose a
simpler, more robust approach: We extract keypoints from the first frame of
the video using a keypoint detector, and track them using an optical-flow-
based tracker. This produces a set of 2D keypoint tracks
$\tilde{\mathbf{k}}_{1:T}$, and allows physical parameter+pose estimation to
be performed for each track independently.

Selection Once the physical parameters and pose are estimated for each
keypoint track, the best tracklet can be identified by isolating the highest
projection likelihood (3.1). However, in order to prevent trivial keypoint
tracklets from being chosen (since a static keypoint will easily attain
maximal likelihood), we add a temporal entropy term to the likelihood, such as
the temporal standard deviation of the observed trajectory, resulting in the
following selection criterion:

$\displaystyle n^{\text{best}}$ $\displaystyle=\operatorname*{arg\,max}_{n\in
1..N}\;p(\tilde{\mathbf{k}}^{n}_{1:T}|\bm{\theta},R,\mathbf{t})+\text{Stddev}_{t}(\tilde{\mathbf{k}}^{n}_{1:T})$
(3)

This finds the highest entropy trajectory that satisfies the physical motion
constraints.

The full V-SysId procedure is depicted in Fig. 2 and pseudocode is shown in
Algorithm 1 in Appendix B.

Bouncing ball with unknown velocity, initial height, and restitution
coefficient.

Archimedes spiral with unknown radius, radius increase rate, and angular
velocity.

Figure 3: Discovered object and 3D perspective given the only the family of
equations above as weak supervision. Top: Example bouncing ball scene. More
scenes can be found in Fig. 10 in the Appendix. Bottom: Spiral robot arm end-
effector in a real lab setting.

Inference at run-time Once the V-SysId procedure is complete, keypoints are
available for the objects of interest in each frame in the video. These can be
treated as pseudo-ground truth keypoints, and used to train a neural network
(or another visual object detector) by supervised learning, in order to
perform fast keypoint detection at test-time.

##  4 Experiments

Keypoint detection and tracking: We detect keypoints in the first frame by
using taking the locations of a 10x10 grid across the frame, and use the KLT
algorithm to track these across the video. We show comparisons between grid,
ORB, SuperPoint and LF-Net keypoint detectors in Appendix C.

Track filtering: Since the grid keypoint detector extracts hundreds of
keypoints, we remove tracks whose length is less than 60% of the full video,
and whose temporal stddev (3) is less than 10 pixel, prior to optimization.
This reduces computation, as physical parameter + pose estimation is performed
on only the most feasible tracks.

Physical parameter estimation: The gradient-based BFGS
Fletcher2000PracticalOptimization  is used with numerical derivatives for
physical parameter optimization. Although Belbute-Peres2018End-to-EndControl
provides an elegant method for analytical differentiation through contacts, we
found it much harder to implement, and ultimately slower, than simple BFGS.
Since the equations of motion considered here are planar, the $z$ component of
$\mathbf{v}_{0}$ is constrained to $0$. The remaining parameters are
learnable.

On the first iteration, the initial position $\mathbf{p}_{0}$ is set to be the
reprojection of the first 2D keypoint $\tilde{\mathbf{k}_{0}}$ onto the $z=5$
plane in world coordinates. This results in an initial position whose camera
projection is the first keypoint. The initial velocity is
$\mathbf{v}_{0}=[0,0,0]$. We found these settings essential to avoid local
minima in the incremental optimization.

Camera pose estimation: BFGS is also used with finite differencing for the
camera pose optimization step. The parametrization of $R$ on pitch and yaw
provides a smooth objective that is easy to optimize, whereas we found the PnP
algorithm to result in large and not necessarily optimal jumps between steps.
We initialize the camera pose parameters as $\alpha=0$, $\beta=0$, and
$\mathbf{t}=[0,0,0]$.

Curriculum-based optimization: We use 25 input frames to start the
optimization, adding 10 frames per iteration until reaching the full length of
the sequence.

Object of interest: bouncing ball

Object of interest: robot arm end-effector

Figure 4: Left: Keypoint tracks propsed by a grid keypoint detector + KLT
tracker (short or static tracks not shown here for clarity). Right: Subset of
the extracted keypoint tracks (red) and projected fitted trajectories (blue),
with the corresponding projection loglikelihood, entropy, and their sum, over
each plot. See Fig. 9 for a more comprehensive illustration.

###  4.1 Environments

Franka Emika Panda Robot: This sequence consists of a multi-joint robot arm
(Franka Emika Panda) in a laboratory setting, where the goal is to find the
end-effector’s 3D location and the camera pose relative to this. The end-
effector was programmed to follow an archimedes spiral in an unknown 2D plane.
The spiral is described by: $r=a+b\cdot
t\quad;\quad\theta=\theta_{0}+\omega\cdot t$ where $r$, $a$, $b$,
$\theta_{0}$, $\omega$ are unknown parameters, to be learned by V-SysId, and
$t$ is the time in seconds. A sequence of frames for this environment can be
seen on Fig. 3, bottom. The video is 250 frames long, with a resolution of
$640\times 480$.

Simulated bouncing ball: This environment consists of a simulated bouncing
ball with moving distractor objects. The bouncing ball follows the equation of
motion:

$\displaystyle\vspace{-2mm}\begin{cases}a_{y}=-g\;,\;\text{if}\;y>\text{floor}\\\
v_{x}=v_{x_{0}};\quad v_{z}=0;\quad
v_{y}=-\epsilon\,v_{y}\;,\;\text{if}\;y=\text{floor}\\\
\end{cases}\vspace{-2mm}$ (4)

where $a$ is the acceleration, $v$ is the velocity, $y$ is the ball height,
$\epsilon\in[0,1]$ is the restitution coefficient, and $g=9.8$ is the gravity.
The ball moves in the $z=z_{0}$ plane with constant horizontal velocity, with
the pose parameters $R,\mathbf{t}$ being responsible for correctly inferring
the location of this plane relative to the camera. Photorealistic scenes are
rendered in Blender following the Clevr protocol johnson2017clevr , and
trajectories are rolled out using Euler integration.

There are two distractor objects on the floor scene, one moving in a circle,
and another in a straight line. This environment is used to obtain thorough
quantitative results regarding the physical parameter and camera pose
estimation abilities of V-SysId. To this end we generate 108 sequences along
the following factors of variation: initial height; initial horizontal
velocity; restitution coefficient; camera location; moving/static distractor
objects. The physical parameters $y_{0}$, $v_{y_{0}}$, $v_{x_{0}}$, $\eta$,
and floor height are unknown, and discovered by the optimization process of
V-SysId. The sequences are 120 frames long, with a resolution of $320\times
240$.

Figure 5: Visualization of the curriculum-based optimization iterations for
the spiral robot (top) and bouncing ball (bottom) scenes. The red line
corresponds to the extracted keypoint track and the solid blue line
corresponds to the trajectory with parameters estimated so far. The dashed
blue line corresponds to the predicted trajectory over the full length of the
sequence, under the parameters estimated so far. We can see that the
curriculum-based optimization progressively improves the physical parameter
and pose estimates.

###  4.2 Visualizing keypoint proposal and optimization

We start by visually exploring the results obtained by V-SysId on the spiral
robot and bouncing ball datasets. Fig. 3 shows the keypoints discovered for
two of the scenes. These show that V-SysId correctly identifies objects of
interest according to the given equation of motion.

The keypoint proposal and selection process is visualized further in Fig. 4.
Fig. 4 (left) shows the proposed keypoint tracks extracted at the proposal
stage (Sec 3.2), and Fig. 4 (right) shows the results obtained by the
optimization process (Sec 3.1) on a subset of these, ordered by their
selection criterion score (the third number above each plot). The trajectory
chosen by V-SysId according to the maximum entropy criterion is labeled as
“Best”. These figures highlight several important points: Firstly, V-SysId is
successful despite the large number of distractor keypoints from the various
moving parts of the scene (most notable in the robot arm sequence). Secondly
and crucially, the optimization process and the maximum entropy criterion are
able to fit and identify the best trajectory, correctly discovering the object
corresponding to the motion of interest.

In order to further understand the curriculum-based optimization process, we
visualize the optimization iterations of two keypoint tracks selected by
V-SysId in Fig. 5. We can see that upon completion (2nd column), the
orientation of the trajectory in 3D space is correctly identified by the
model, and that each iteration progressively adjusts both the trajectory’s
shape (parametrized by the physical parameters) and the camera pose. This
leads to a stable optimization procedure where both physical parameters and
camera pose are identified.

###  4.3 Evaluating parameter estimation

Even though the scale is generally unidentifiable (this and other limitations
are discussed in Sec. A), in the case of a bouncing ball both the initial
height and the restitution coefficient are exactly identifiable. This allows
us to compare their learned values to the ground truth values used for the
simulations. In addition, we can compare the camera angles identified to those
used in simulation in order to evaluate the quality of the extrinsic camera
calibration.

The percentage error in restitution coefficient, initial height (distance to
floor), and camera angle relative to the simulation ground-truth can be seen
in Table 1. We can see that all parameters are found with a good degree of
accuracy, with physical parameters being slightly more accurate than the
camera pose. Notably, the errors are similar with and without moving
distractors (within 95% confidence intervals), showing that V-SysId is able to
correctly identify the object of interest even in the presence of distractor
objects.

In order to highlight the importance of the curriculum-based optimization
strategy, we compare the projection likelihood using our incremental alternate
optimization with alternate optimization using the full sequence at every
step. Averaging over the bouncing ball scenes, we obtain projection RMSE
(pixels) of $-9.31$ and $-109.35$, respectively. A similar decrease in
performance was observed when using CEM and BO optimizers. This shows that
gradually increasing sequence length and using a gradient-based optimizer is
key to the convergence of V-SysId.

Distractors |  Restitution coefficient (%) |  Initial height in 3D (%) |  Camera angle (°)
---|---|---|---
With | $3.8\pm 1.5$ | $9.7\pm 4.0$ | $8.0\pm 1.8$
Without | $2.7\pm 0.8$ | $6.7\pm 3.0$ | $9.9\pm 2.6$
Table 1: Relative error (percentage) between the ground-truth simulation
physical parameters and camera pose, and those estimated by V-SysId, for the
bouncing ball scene. Error bounds correspond to a 95% confidence interval.

###  4.4 Tracking by supervised keypoint detection

Once detected, the keypoints discovered by V-SysId can be used as pseudo-
ground-truth to train a supervised keypoint detector. For the bouncing ball
dataset, the training set consists of 2838 pseudo-labeled frames, and the test
set consists of 948 hand-labeled frames from unseen scene configurations. For
the robot dataset, the training set consists of 250 pseudo-labeled frames, and
the test set consists of 150 hand-labeled frames from unseen end-effector
positions. For the supervised keypoint detector, we use a fully convolutional
neural network with 6 ReLU layers with 32 channels, with stride 2 on the 3rd
layer, and 2 output channels with 2D softmax activation. These maps are
converted to $[x,y]$ coordinates by taking the softmax-weighted mean over the
output coordinate grid, as per Jakab2018UnsupervisedGeneration . The input
images have a downsampling factor of 4 relative to the original frame
resolutions, but we report the keypoint error in the original image space. We
train the networks for 20 epochs with batch size 16, and Adam Kingma
(learning rate $3\times 10^{-4}$).

Results are shown in Table 2. The supervised keypoint detector produces highly
accurate detections, confirming the quality and usability of the keypoints
discovered by V-SysId even on small datasets of high-resolution scenes.

Environment | RMSE (pixel distance)
---|---
Simulated bouncing ball ($240\times 320$) | $8.41\pm 1.50$
Spiral robot ($480\times 640$) | $3.89\pm 0.45$
Table 2: Detection error on the held-out test set of the keypoints extracted
by the inference neural network, after training using the keypoints discovered
by V-SysId as supervision. Bounds correspond to 95% confidence interval.

We recommend reading of the supplementary materials for additional ablations
experiments, and an application of V-SysId to breathing rate estimation from
video.

###  4.5 ROI discovery in breathing videos using RANSAC

Setup To further demonstrate the applicability of V-SysId to real world
scenarios, we collected 8 videos of people breathing under different pose,
lightning, clothing and distractor settings, with the goal of discovering the
relevant region region of the image and using it for breathing rate
identification. The true breathing rate was obtained by manual annotation.
Videos contain between 150 and 300 frames, at 30 fps and 480x640 resolution.

Unlike seminal work in video-based physiology and plethysmography boer2010slp
, V-SysId does not require careful hand selection of the regions of interest
and is robust to the existence of distractor motions in the scene. V-SysId
simultaneously identifies the region of interest (here, the set of relevant
keypoints, rather than a single one) corresponding to sinusoidal motion, and
the underlying breathing rate.

Results We have seen how single keypoint discovery can be achieved using
V-SysId, but the algorithm can be easily modified to allow discovery of sets
of keypoints constituting a region-of-interest. We use the chest video dataset
as a prototypical application. The goal is to discover the keypoints in the
video corresponding to sinusoidal motion. We start by extracting keypoint
tracks as in Stage 1 of V-SysId (filtering out any tracks with a temporal
stddev less than 0.7), and transform these 2D tracks into 1D timeseries by
taking the projection onto the 1st PCA component of the timeseries (i.e. the
2D direction of highest variance). Each timeseries is standardised, and fit to
a sinusoid as per Stage 2 of V-SysId (without the 3D component). In order to
identify the best set of tracks, we use a RANSAC inlier count, by measuring
the error between a track’s fitted sinusoid and all the other extracted
tracks, and considering a track an inlier if the MSE is below 0.75. The best
track is chosen according to a modified maximum entropy criterion in Stage 3,
where the likelihood term is replaced by the inlier count. The ROI is defined
as the set of inlier tracks of the best track.

Fig. 6 (top) shows the keypoints discovered for the 8 videos, with Fig. 6
(bottom) showing the timeseries and its sinusoidal fit for one of the
keypoints in the ROI. The model correctly identifies keypoints corresponding
to the chest area, while ignoring distractor and lower-body keypoints.
Comparing the respiratory periods identified with V-SysId with the annotated
values results in an MSE of 0.016 (in seconds/breath). In contrast, a baseline
that uses the mean of the true rates for all videos obtains an MSE of 0.085.
These results demonstrate the accuracy of V-SysId for physical parameter
estimation from an unknown region of interest, using only the knowledge that
the motion of interest is sinusoidal as supervision.

Figure 6: Top: Green dots correspond to keypoints identified by V-SysId as
relevant for determining the breathing rate. The red dots are discarded
keypoints. Note that some the videos contain distractors that move in the
scene (rollouts of scenes are shown in Fig. 11 in Appendix). V-SysId with
RANSAC is able to automatically discover regions of interest. Bottom:
Timeseries (blue) and sinusoidal fit (orange) of one keypoint in the ROI for
each of the scenes (same position in the $2\times 4$ grid)

##  5 Conclusion and future work

This paper has introduced V-SysId, a 3-stage method for dynamics-constrained
keypoint discovery and system identification, which alternates between maximum
likelihood extrinsic camera calibration and maximum likelihood physical
parameter estimation for motion tracks detected in video. We enhance the
stability of this optimization through the inclusion of a curriculum-based
optimisation strategy, alongside a maximum entropy selection criterion for
keypoint identification. Future avenues of work include extensions to multiple
interacting objects, rigid or fluid body dynamics from video, and
incorporation with a neural network for material and volume inference from
vision.

## References

  * (1) Ian Abraham, Gerardo De, La Torre, and Todd D Murphey.  Model-Based Control Using Koopman Operators.  In RSS, 2017.
  * (2) Martin Asenov, Michael Burke, Daniel Angelov, Todor Davchev, Kartic Subr, and Subramanian Ramamoorthy.  Vid2Param: Modelling of Dynamics Parameters from Video.  In ICRA, 2019.
  * (3) Shai Avidan and Amnon Shashua.  Trajectory triangulation: 3D reconstruction of moving points from a monocular image sequence.  PAMI, 2000.
  * (4) Chris L. Baker, Julian Jara-Ettinger, Rebecca Saxe, and Joshua B. Tenenbaum.  Rational quantitative attribution of beliefs, desires and percepts in human mentalizing.  Nature Human Behaviour, 2017.
  * (5) Peter W. Battaglia, Jessica B. Hamrick, and Joshua B. Tenenbaum.  Simulation as an engine of physical scene understanding.  PNAS, 2013.
  * (6) Filipe De A Belbute-Peres, Kevin A Smith, Kelsey R Allen, Joshua B Tenenbaum, and J Zico Kolter.  End-to-End Differentiable Physics for Learning and Control.  In NIPS, 2018.
  * (7) Rune Brincker, Lingmi Zhang, and Palle Andersen.  Modal identification of output-only systems using frequency domain decomposition.  Smart materials and structures, 2001.
  * (8) Steven L Brunton, Joshua L Proctor, and J Nathan Kutz.  Discovering governing equations from data by sparse identification of nonlinear dynamical systems.  Proceedings of the national academy of sciences, 2016.
  * (9) Miles Cranmer, Alvaro Sanchez-Gonzalez, Peter Battaglia, Rui Xu, Kyle Cranmer, David Spergel, and Shirley Ho.  Discovering Symbolic Models from Deep Learning with Inductive Biases.  In NeurIPS, 2020.
  * (10) Neha Das, Sarah Bechtle, Todor Davchev, Dinesh Jayaraman, Akshara Rai, and Franziska Meier.  Model-Based Inverse Reinforcement Learning from Visual Demonstrations.  In CoRL, 2020.
  * (11) Philip David, Daniel Dementhon, Ramani Duraiswami, and Hanan Samet.  SoftPOSIT: Simultaneous Pose and Correspondence Determination.  IJCV, 2004.
  * (12) Willem de Boer, Joan Lasenby, Jonathan Cameron, Rich Wareham, Shiraz Ahmad, Charlotte Roach, Ward Hills, and Richard Iles.  Slp: A zero-contact non-invasive method for pulmonary function testing.  In BMVC, 2010.
  * (13) Daniel DeTone, Tomasz Malisiewicz, and Andrew Rabinovich.  SuperPoint: Self-Supervised Interest Point Detection and Description.  In CVPR, 2018.
  * (14) Philipp Fischer, Alexey Dosovitskiy, Eddy Ilg, Philip Hausser, Caner Hazırbas, Vladimir Golkov, Patrick Van Der Smagt, Daniel Cremers, and Thomas Brox.  FlowNet: Learning Optical Flow with Convolutional Networks.  In ICCV, 2015.
  * (15) Andrew W Fitzgibbon and Andrew Zisserman.  Multibody Structure and Motion: 3-D Reconstruction of Independently Moving Objects.  2000\.
  * (16) R. Fletcher.  Practical Methods of Optimization.  John Wiley & Sons, Ltd, 2000.
  * (17) Anand Gopalakrishnan, Sjoerd Van Steenkiste, and J Urgen Schmidhuber.  Usupervised Object Keypoint Learning using Local Spatial Predictability.  arxiv.org/abs/2011.12930, 2020.
  * (18) Sam Greydanus, Misko Dzamba, and Jason Yosinski.  Hamiltonian Neural Networks.  In NeurIPS, 2019.
  * (19) Vincent Le Guen and Nicolas Thome.  Disentangling Physical Dynamics from Unknown Factors for Unsupervised Video Prediction.  In CVPR, 2020.
  * (20) Jessica B Hamrick, Peter W Battaglia, Thomas L Griffiths, and Joshua B Tenenbaum.  Inferring mass in complex scenes by mental simulation.  Cognition, 2016.
  * (21) Mei Han and Takeo Kanade.  Multiple Motion Scene Reconstruction with Uncalibrated Cameras.  In PAMI, 2003.
  * (22) Tomas Jakab, Ankush Gupta, Hakan Bilen, and Andrea Vedaldi.  Unsupervised Learning of Object Landmarks through Conditional Image Generation.  In NIPS, 2018.
  * (23) Tomas Jakab, Ankush Gupta, Hakan Bilen, and Andrea Vedaldi.  Learning Landmarks from Unaligned Data using Image Translation.  7 2019.
  * (24) Miguel Jaques, Michael Burke, and Timothy Hospedales.  Physics-as-Inverse-Graphics: Unsupervised Physical Parameter Estimation from Video.  In ICLR, 2020.
  * (25) Miguel Jaques, Michael Burke, and Timothy Hospedales.  NewtonianVAE: Proportional Control and Goal Identification from Pixels via Physical Latent Spaces.  In CVPR, 2021.
  * (26) You-Yi Jau, Rui Zhu, Hao Su, and Manmohan Chandraker.  Deep Keypoint-Based Camera Pose Estimation with Geometric Constraints.  In IROS, 2020.
  * (27) Marija Jegorova, Joshua Smith, Michael Mistry, and Timothy Hospedales.  Adversarial generation of informative trajectories for dynamics system identification.  In IROS, 2020.
  * (28) Justin Johnson, Bharath Hariharan, Laurens Van Der Maaten, Li Fei-Fei, C Lawrence Zitnick, and Ross Girshick.  Clevr: A diagnostic dataset for compositional language and elementary visual reasoning.  In CVPR, 2017.
  * (29) Jer-Nan Juang and Richard S Pappa.  An eigensystem realization algorithm for modal parameter identification and model reduction.  JGCD, 1985.
  * (30) Jeremy Yirmeyahu Kaminski and Mina Teicher.  General trajectory triangulation.  In ECCV. 2002.
  * (31) Rama Kandukuri, Jan Achterhold, Michael Moeller, and Joerg Stueckler.  Learning to Identify Physical Parameters from Video Using Differentiable Physics.  In GCPR, 2020.
  * (32) Maximilian Karl, Maximilian Soelch, Justin Bayer, and Patrick Van Der Smagt.  Deep Variational Bayes Filters: Unsupervised Learning of State Space Models from Raw Data.  In ICLR, 2017.
  * (33) Diederik Kingma and Jimmy Ba.  Adam: A Method for Stochastic Optimization.  In ICLR, 2014.
  * (34) Tejas Kulkarni, Ankush Gupta, Catalin Ionescu, Sebastian Borgeaud, Malcolm Reynolds, Andrew Zisserman, and Volodymyr Mnih.  Unsupervised Learning of Object Keypoints for Perception and Control.  In NIPS, 2019.
  * (35) Yunzhu Li, Hao He, Jiajun Wu, Dina Katabi, and Antonio Torralba.  Learning compositional koopman operators for model-based control.  ICLR, 2020.
  * (36) Yunzhu Li, Toru Lin, Kexin Yi, Daniel M Bear, Daniel L K Yamins, Jiajun Wu, Joshua B Tenenbaum, and Antonio Torralba.  Visual Grounding of Learned Physical Models.  In ICML, 2020.
  * (37) David G Lowe.  Distinctive Image Features from Scale-Invariant Keypoints.  IJCV, 2004.
  * (38) Matthias Minderer, Chen Sun, Ruben Villegas, Forrester Cole, Kevin Murphy, and Honglak Lee.  Unsupervised Learning of Object Structure and Dynamics from Videos.  In NIPS, 2019.
  * (39) Jonas Mockus.  Bayesian Approach to Global Optimization.  Springer, 1989.
  * (40) J. Krishna Murthy, Miles Macklin, Florian Golemo, Vikram Voleti, Linda Petrini, Martin Weiss, Breandan Considine, Jérôme Parent-Lévesque, Kevin Xie, Kenny Erleben, Liam Paull, Florian Shkurti, Derek Nowrouzezahrai, and Sanja Fidler.  gradSim: Differentiable simulation for system identification and visuomotor control, 9 2020.
  * (41) Yuki Ono, Eduard Trulls, Pascal Fua, and Kwang Moo Yi.  LF-Net: Learning Local Features from Images.  In NeurIPS, 2018.
  * (42) Yi-Ling Qiao, Junbang Liang, Vladlen Koltun, and Ming C Lin.  Scalable Differentiable Physics for Learning and Control.  In ICML, 2020.
  * (43) Santhosh K Ramakrishnan, Swarna Kamlam Ravindran, Anurag Mittal, and Iit Madras.  CoMaL Tracking: Tracking Points at the Object Boundaries.  In CVPR, 2016.
  * (44) Edward Rosten and Tom Drummond.  Machine learning for high-speed corner detection.  In ECCV, 2006.
  * (45) Reuven Y. Rubinstein.  Optimization of computer simulation models with rare events.  EJOR, 1997.
  * (46) Ethan Rublee, Vincent Rabaud, Kurt Konolige, and Gary Bradski.  ORB: An efficient alternative to SIFT or SURF.  In ICCV, 2011.
  * (47) Alvaro Sanchez-Gonzalez, Jonathan Godwin, Tobias Pfaff, Rex Ying, Jure Leskovec, and Peter W Battaglia.  Learning to Simulate Complex Physics with Graph Networks.  In ICML, 2020.
  * (48) Davide Scaramuzza, Friedrich Fraundorfer, Marc Pollefeys, and Roland Siegwart.  Absolute scale in structure from motion from a single vehicle mounted camera by exploiting nonholonomic constraints.  In ICCV, 2009.
  * (49) Supasorn Suwajanakorn, Noah Snavely, Jonathan Tompson, and Mohammad Norouzi.  Discovery of Latent 3D Keypoints via End-to-end Geometric Reasoning.  In NeurIPS, 2018.
  * (50) Carlo Tomasi and Takeo Kanade.  Detection and Tracking of Point Features.  Technical report, 1991.
  * (51) Peter Toth, Danilo J Rezende, Andrew Jaegle, Sébastien Racanière, Aleksandar Botev, and Irina Higgins.  Hamiltonian Generative Networks.  In ICLR, 2020.
  * (52) Tomer Ullman, Andreas Stuhlmuller, Noah Goodman, and Joshua Tenenbaum.  Learning Physics from Dynamical Scenes.  In CogSci, 2014.
  * (53) Sudheendra Vijayanarasimhan, Susanna Ricco, Cordelia Schmid, Rahul Sukthankar, and Katerina Fragkiadaki.  SfM-Net: Learning of Structure and Motion from Video.  arXiv preprint arXiv:1704.07804, 2017.
  * (54) Xingkui Wei, Yinda Zhang, Zhuwen Li, Yanwei Fu, and Xiangyang Xue.  DeepSFM: Structure From Motion Via Deep Bundle Adjustment.  In ECCV, 2020.
  * (55) Jiajun Wu, Joseph J Lim, Hongyi Zhang, Joshua B Tenenbaum, and William T Freeman.  Physics 101: Learning physical object properties from unlabeled videos.  In BMVC, 2016.
  * (56) Jiajun Wu, Erika Lu, Pushmeet Kohli, William T Freeman, and Joshua B Tenenbaum.  Learning to See Physics via Visual De-animation.  In NIPS, 2017.
  * (57) Jiajun Wu, Ilker Yildirim, J.J. Lim, W.T. Freeman, and J.B. Tenenbaum.  Galileo : Perceiving Physical Object Properties by Integrating a Physics Engine with Deep Learning.  In NIPS, 2015.
  * (58) Yaofeng Desmond Zhong and Naomi Ehrich Leonard.  Unsupervised Learning of Lagrangian Dynamics from Images for Prediction and Control.  In NeurIPS, 2020.

##  Appendix A Discussion: Challenges and limitations

Scale unidentifiability Due to the projection operation
$\mathbf{k}_{t}(\bm{\theta},R,\mathbf{t})=[\tilde{\mathbf{p}}_{x,t}/\tilde{\mathbf{p}}_{z,t},\tilde{\mathbf{p}}_{y,t}/\tilde{\mathbf{p}}_{z,t}]$,
the 3D trajectory $\mathbf{p}_{1:T}$ can only be determined up to a scale
parameter. For this reason, we evaluate the correlation between the true- and
learned parameters, not the error. This is also the metric used by GALILEO and
Physics101 when doing physical parameter estimation from visual trajectories.
Although scale unidentifiability leads to the existence of infinitely many
solutions for $\bm{\theta}$ and $\mathbf{t}$, and therefore instability with
joint optimization, our use of alternate optimization steps guarantees that
the algorithm converges to a single solution, as $\bm{\theta}$ and
$\mathbf{t}$ are optimized conditioned on one another, not jointly (this is
akin to the Expectation-Maximization algorithm, where a marginal distribution
is maximized via alternate optimization of conditionals).

Broken trajectories and occlusions: In settings where classical keypoint
detectors are unreliable, one can either use state-of-the-art pretrained
keypoint networks, like SuperPoint DeTone2018SuperPoint:Description  and LF-
Net Ono2018LF-Net:Images , or pretrain an unsupervised keypoint discovery
network Jakab2018UnsupervisedGeneration ; Minderer2019UnsupervisedVideos ;
Kulkarni2019UnsupervisedControl . However, we found show that a simple grid
keypoint detector yielded more reliable tracks than using classic (SIFT, ORB,
FAST), or modern (SuperPoint, LF-Net) keypoint detectors.

In settings where standard optical flow computation is unreliable, more recent
models (eg. FlowNet Fischer2015FlowNet:Networks ) could be used to provide
flow estimates to the KLT tracker. More recent improvements to the KLT tracker
(eg. CoMaL Ramakrishnan2016CoMaLBoundaries ) could also be used. The multi-
stage, modular nature of our pipeline allows for the easy replacement of
individual components, although we found standard optical flow computation to
work very well in practice.

Camera roll set to zero: We found that setting the camera roll angle to zero
greatly stabilized the optimization procedure. While this might be perceived
as too strong of a constraint on the model, in the vast majority of real
settings the camera has zero roll (i.e. it’s rare for the camera to rotate
around its projection axis). Therefore, imposing this constraint does not
reduce the applicability of our method in the vast majority of cases, while
providing improved results. Naturally, allowing roll optimization would make
the model more general, but this is left as future work.

##  Appendix B V-SysId pseudocode

Algorithm 1 V-SysId

0: Video $V$ of length $T$

0: Equation of motion $f$ of the object of interest

0: KeypointTrackExtractor # function that outputs a set of keypoint track
proposals

0: Trajectory, physical parameters, and camera pose of the object of interest
# Get $N$ keypoint track proposals

Tracks $\leftarrow$ KeypointTrackExtractor($V$) # Fit physical parameters and
camera pose to trajectory

$\text{SelectionCriterion}\leftarrow[\,]$

$\text{Params}\leftarrow[\,]$

for $n\in\\{1...N\\}$ do

$\tilde{\mathbf{k}}_{1:T}\leftarrow\text{Tracks}[n]$

Initialize $\alpha\leftarrow 0$, $\beta\leftarrow 0$,
$\mathbf{t}\leftarrow[0,0,0]$, $\mathbf{v}_{0}\leftarrow[0,0,0]$;

Initialize $\mathbf{p}_{0}$ as the projection of $\tilde{\mathbf{k}}_{0}$ onto
the $z=5$ plane in world coordinates;

Initialize $\bm{\eta}$ to some sensible initial values (setting dependent);

for $t\in\\{1...T\\}$ do

$\bm{\theta}\leftarrow\operatorname*{arg\,max}_{\bm{\theta}}\;p(\tilde{\mathbf{k}}_{1:t}|\bm{\theta},R,\mathbf{t})$

$R,\mathbf{t}\leftarrow\operatorname*{arg\,max}_{R,\mathbf{t}}\;p(\tilde{\mathbf{k}}_{1:t}|\bm{\theta},R,\mathbf{t})$

end for

Append $\\{\bm{\theta},R,\mathbf{t}\\}$ to Params

Append the scalar
$p(\tilde{\mathbf{k}}_{1:t}|\bm{\theta},R,\mathbf{t})+H(\tilde{\mathbf{k}}_{1:T})$
to SelectionCriterion

end for # Trajectory selection

$n^{*}\leftarrow\operatorname*{arg\,max}\text{SelectionCriterion}$

$\tilde{\mathbf{k}}^{*}=\text{Tracks}[n^{*}]$

$\bm{\theta}^{*},R^{*},\mathbf{t}^{*}=\text{Params}[n^{*}]$

return $\tilde{\mathbf{k}}^{*},\bm{\theta}^{*},R^{*},\mathbf{t}^{*}$

###  B.1 Evaluating future trajectory prediction

We now evaluate the ability of the optimization process of V-SysId to perform
accurate tracking and prediction given a sequence of correct keypoints. At
each optimization iteration, we rollout the trajectory under the current
parameters (i.e. those learned with the fraction of the sequence up to that
iteration), and measure the prediction error relative to the observed keypoint
track. This is done in hindsight only for the trajectory chosen by V-SysId,
although it could equally be done for a keypoint sequence inferred by a test-
time inference neural network. The results for the bouncing ball and spiral
robot are shown in Fig. 7. The curves show that the optimization process
quickly converges to correct system identification, leading to correct
trajectory prediction after only 2 seconds of input.

Figure 7: Left: Future trajectory prediction error under estimated parameters
as a function of input length.

##  Appendix C Comparison of keypoint detectors

Here we provide a visual comparison of the trajectory proposals obtained using
grid, ORB, LF-Net and SuperPoint keypoint extractors, in conjunction with a
KLT tracker. Fig. 8 shows this comparison for the bouncing ball dataset (after
filtering for short and static tracks). It can be seen that despite its
simplicity, the grid extractor performs just as well as the more modern
keypoint detectors, while running over an order of magnitude faster.

Figure 8: More visualization of the discovered object in various bouncing ball
scenes.

##  Appendix D Further visualizations

See next page.

Bouncing ball

Robot arm end-effector

\begin{overpic}[width=433.62pt]{figures/spiral_trajs_and_stats.pdf}
\put(15.0,24.0){\small{\rotatebox{10.0}{\noindent\hbox{}\hfill{{\hbox{\begin{tabular}[c]{@{}c@{}}Best\end{tabular}}}}\hfill\hbox{}}}}
\put(60.0,5.0){\scriptsize{\rotatebox{10.0}{\noindent\hbox{}\hfill{{\hbox{\begin{tabular}[c]{@{}c@{}}Distractor\end{tabular}}}}\hfill\hbox{}}}}
\put(85.0,5.0){\scriptsize{\rotatebox{10.0}{\noindent\hbox{}\hfill{{\hbox{\begin{tabular}[c]{@{}c@{}}Distractor\end{tabular}}}}\hfill\hbox{}}}}
\put(15.0,67.0){\small{\rotatebox{10.0}{\noindent\hbox{}\hfill{{\hbox{\begin{tabular}[c]{@{}c@{}}Best\end{tabular}}}}\hfill\hbox{}}}}
\put(60.0,67.0){\scriptsize{\rotatebox{10.0}{\noindent\hbox{}\hfill{{\hbox{\begin{tabular}[c]{@{}c@{}}Distractor\end{tabular}}}}\hfill\hbox{}}}}
\put(85.0,67.0){\scriptsize{\rotatebox{10.0}{\noindent\hbox{}\hfill{{\hbox{\begin{tabular}[c]{@{}c@{}}Distractor\end{tabular}}}}\hfill\hbox{}}}}
\put(11.0,47.0){\scriptsize{\rotatebox{10.0}{\noindent\hbox{}\hfill{{\hbox{\begin{tabular}[c]{@{}c@{}}Distractor\end{tabular}}}}\hfill\hbox{}}}}
\put(60.0,54.0){\scriptsize{\rotatebox{10.0}{\noindent\hbox{}\hfill{{\hbox{\begin{tabular}[c]{@{}c@{}}Distractor\end{tabular}}}}\hfill\hbox{}}}}
\put(85.0,54.0){\scriptsize{\rotatebox{10.0}{\noindent\hbox{}\hfill{{\hbox{\begin{tabular}[c]{@{}c@{}}Distractor\end{tabular}}}}\hfill\hbox{}}}}
\end{overpic}

Figure 9: Left: Keypoint tracks propsed by a grid keypoint detector + KLT
tracker (short or static tracks not shown here for improved visualization).
Right: Subset of the extracted keypoint tracks (red) and projected fitted
trajectories (blue), with the corresponding projection loglikelihood, entropy,
and their sum, over each plot.

Figure 10: Comparison of various keypoints extractor and trackers on a
bouncing ball scene.

Figure 11: Frames of breathing scenes containing distractors.