Hamiltonian formulation for the theory of gravity
and canonical transformations in extended phase space

T. P. Shestakova

Department of Theoretical and Computational Physics, Southern Federal
University,

Sorge St. 5, Rostov-on-Don 344090, Russia

E-mail<EMAIL_ADDRESS>

Abstract

A starting point for the present work was the statement recently discussed in
the literature that two Hamiltonian formulations for the theory of gravity,
the one proposed by Dirac and the other by Arnowitt – Deser – Misner, may not
be related by a canonical transformation. In its turn, it raises a question
about the equivalence of these two Hamiltonian formulations and their
equivalence to the original formulation of General Relativity. We argue that,
since the transformation from components of metric tensor to the ADM variables
touches gauge degrees of freedom, which are non-canonical from the point of
view of Dirac, the problem cannot be resolved in the limits of the Dirac
approach. The proposed solution requires the extension of phase space by
treating gauge degrees of freedom on an equal footing with other variables and
introducing missing velocities into the Lagrangian by means of gauge
conditions in differential form. We illustrate with a simple cosmological
model the features of Hamiltonian dynamics in extended phase space. Then, we
give a clear proof for the full gravitational theory that the ADM-like
transformation is canonical in extended phase space in a wide enough class of
possible parametrizations.

##  1\. Introduction

It is generally accepted that the problem of formulating Hamiltonian dynamics
for systems with constraints has been solved by Dirac in his seminal papers
[1, 2]. It was Dirac who pointed to the importance of Hamiltonian formulation
for any dynamical theory before its quantization [3]. Other approaches, such
as the Batalin – Fradkin – Vilkovisky (BFV) path integral approach [4, 5, 6]
follow the Dirac one in what concerns the rule of constructing a Hamiltonian
and the role of constraints as generators of transformations in phase space.
It is believed that Dirac generalized Hamiltonian dynamics is equivalent to
Lagrangian dynamics of original theory. However, even for electrodynamics the
constraints do not generate a correct transformation for zero component of
vector potential, $A_{0}$. The same situation we face in General Relativity,
since gravitational constraints cannot produce correct transformations for
$g_{00}$, $g_{0\mu}$ components of metric tensor. In fact, it means that the
group of transformations generated by constraints differs from the group of
gauge transformations of the original theory. Some authors have tried to
remedy this shortcoming by modifying the Dirac approach and proposing some
special prescriptions how the generator should be constructed (see, for
example, [7, 8]). Until now this problem has not attracted much attention
mainly because that it touches only transformations of gauge variables which,
according to conventional viewpoint, are redundant and must not affect the
physical content of the theory. It will be demonstrated in this paper that the
role of gauge degrees of freedom may be more significant that it is usually
thought, and the difference in the groups of transformations is the first
indication to the inconsistence of the theory.

Historically, while constructing Hamiltonian dynamics for gravitational field
theorists used various parametrizations of gravitational variables. Dirac
dealt with original variables, which are components of metric tensor [3],
whereas the most famous parametrization is probably that of Arnowitt – Deser –
Misner (ADM) [9], who expressed $g_{00}$, $g_{0\mu}$ through the lapse and
shift functions. To give another example, let us mention the work by Faddeev
[10] where quite specific variables $\lambda^{0}=1/h^{00}+1$,
$\lambda^{i}=h^{0i}/h^{00}$, $q^{ij}=h^{0i}h^{0j}-h^{00}h^{ij}$,
$h^{\mu\nu}=\sqrt{-g}g^{\mu\nu}$ were introduced. From the point of view of
Lagrangian formalism, all the parametrizations are rightful, and the
correspondent formulations are equivalent. Meanwhile, it has been shown in
[11] that components of metric tensor and the ADM variables are not related by
a canonical transformation. In other words, it implies that the Dirac
Hamiltonian formulation for gravitation and the ADM one are not equivalent,
though it is believed that each of them is equivalent to the Einstein
(Lagrangian) formulation. There exists the contradiction that again witnesses
about the incompleteness of the theoretical foundation.

The purpose of the present paper is to demonstrate that this contradiction can
be resolved if one treats gauge gravitational degrees of freedom on an equal
footing with physical variables in extended phase space. The idea of extended
phase space was put forward by Batalin, Fradkin and Vilkovisky [4, 5, 6] who
included integration over gauge and ghost degrees of freedom in their
definition of path integral. However, in their approach gauge variables were
still considered as non-physical, secondary degrees of freedom playing just an
auxiliary role in the theory. To construct Hamiltonian dynamics for a
constrained system which would be completely equivalent to Lagrangian
formulation, we need to take yet another step: we should introduce into the
Lagrangian missing velocities corresponding to gauge variables by means of
special (differential) gauge conditions. It actually extends the phase space
of physical degrees of freedom.

In Section 2 a mathematical formulation of the problem will be given. We shall
see that non-equivalence of Hamiltonian formulations for different
parametrizations prevents from constructing a generator of transformation in
phase space which would produce correct transformations for any
parametrizations. These ideas will be illustrated in Section 3 for a simple
model with finite number of degrees of freedom. The mentioned above algorithms
[7, 8] work correctly only for some parametrizations. One possible point of
view (advocated, in particular, in [11]) is that only these parametrizations
should be allowed while all other, not related with the first ones by
canonical transformations, should be prohibited, including the ADM
parametrization. However, imposing any limitations on admissible
parametrizations or transformations does not seem to be a true solution to the
problem. In Section 4 the outline of Hamiltonian dynamics in extended phase
space will be presented, and in Section 5 it will be demonstrated for the full
gravitational theory that different parametrizations from a wide enough class
are related by canonical transformations. In particular, it will restore a
legitimate status of the ADM parametrization. We shall discuss the results and
future problems in Section 6.

##  2\. Canonical transformations in phase space

It is generally known that for a system without constraints Lagrangian as well
as Hamiltonian equations maintain their form under transformations to a new
set of generalized coordinates

$q^{a}=v^{a}(Q),$ (1)

where $v^{a}(Q)$ are invertible functions of their arguments. It is easy to
see that any transformation (2.1) correspond to a canonical transformation in
phase space. Indeed, consider a quadratic in velocities Lagrangian

$L=\frac{1}{2}\;\Gamma_{ab}(q)\dot{q}^{a}\dot{q}^{b}-U(q).$ (2)

After the transformation (2.1) the Lagrangian (2.2) would read

$L=\frac{1}{2}\;\Gamma_{cd}(Q)\frac{\partial v^{c}}{\partial
Q^{a}}\frac{\partial v^{d}}{\partial Q^{b}}\dot{Q}^{a}\dot{Q}^{b}-U(Q).$ (3)

New momenta $\\{P_{a}\\}$ are expressed through old momenta $\\{p_{a}\\}$ by
relations

$P_{a}=p_{b}\frac{\partial v^{b}}{\partial Q^{a}}.$ (4)

The transformation (2.1), (2.4) is canonical with the generating function
which depends on new coordinates and old momenta,

$\Phi(Q,\,p)=-p_{a}v^{a}(Q).$ (5)

The equations

$q^{a}=-\frac{\partial\Phi}{\partial p_{a}};\qquad
P^{a}=-\frac{\partial\Phi}{\partial Q^{a}}$ (6)

reproduce exactly the transformation (2.1), (2.4). It is also easy to check
that the transformation (2.1), (2.4) maintains the Poisson brackets

$\\{Q^{a},\,Q^{b}\\}=0,\qquad\\{P_{a},\,P_{b}\\}=0,\qquad\\{Q^{a},\,P_{b}\\}=\delta^{a}_{b}.$
(7)

For a system with constraints, gauge variables (i.e. the variables whose
velocities cannot be expressed in terms of conjugate momenta) do not enter
into the Lagrangian quadratically, and a general transformation like (2.1) may
not be canonical. An example can be found in the theory of gravity by the
transformation from components of metric tensor to the ADM variables,

$g_{00}=\gamma_{ij}N^{i}N^{j}-N^{2},\qquad g_{0i}=\gamma_{ij}N^{j},\qquad
g_{ij}=\gamma_{ij}.$ (8)

This transformation concerns gauge degrees of freedom which, from the
viewpoint of Dirac, are not canonical variables at all. To pose the question,
if the transformation (2.8) is canonical, one should formally extend the
original phase space including into it gauge degrees of freedom and their
momenta. In order to prove non-canonicity of (2.8) it is enough to check that
some of the relations (2.7) are broken. Using the transformation inverted to
(2.8), one can see that $\\{N,\,\Pi^{ij}\\}\neq 0$, where $\Pi^{ij}$ are the
momenta conjugate to $\gamma_{ij}$ (see Equation (152) in [11]). More
generally, let us consider the ADM-like transformation

$N_{\mu}=V_{\mu}(g_{0\nu},\,g_{ij}),\qquad\gamma_{ij}=g_{ij}.$ (9)

Here $V_{\mu}$ are some functions of components of metric tensor (but
$N_{\mu}$ ought not to form 4-vector). A feature of this transformation is
that space components of metric tensor remain unchanged, and so do their
conjugate momenta: $\Pi^{ij}=p^{ij}$. Then

$\left.\\{N_{\mu},\,\Pi^{ij}\\}\right|_{g_{\nu\lambda},p^{\rho\sigma}}=\frac{\partial
V_{\mu}}{\partial g_{ij}}.$ (10)

It is equal to zero if only the functions $V_{\mu}$ do not depend on $g_{ij}$.
This is quite a trivial case when old gauge variables are expressed through
some new gauge variables only, and the ADM transformation (2.8) does not
belong to this class.

One could pose the question: Is it worth considering the equivalence of
different formulations in extended phase space? Would not it better to
restrict ourself by transformations in phase space of original canonical
variables in the sense of Dirac? In the second case, we can prove the
equivalence of equation of motion in Lagrangian and Hamiltonian formalism,
however, we have to fix a form of gravitational constraints by forbidding any
reparametrizations of gauge variables. Determination of the constraints’ form
is of importance for a subsequent procedure of quantization which gives rise
to the problem of parametrization noninvariance (see, for example, [12]). From
the viewpoint of subsequent quantization, the ADM parametrization is more
preferable, since the constraints do not depend on gauge variables in this
case. I would like to emphasize that there are no solid grounds for fixing the
form of the constraints, and, as we shall see in this paper, extension of
phase space enables us to solve the problem of equivalence of Lagrangian and
Hamiltonian formalism for gravity without any restriction on parametrizations.

As it has been already mentioned, the constraints, being considered as
generators of transformations in phase space, do not produce correct
transformation for all gravitational variables. To ensure the full equivalence
of two formulations one has to modify the Dirac prescription, according to
which the generator must be a linear combination of constraints, and replace
it by a more sophisticated algorithm. The known algorithms, firstly, are
relied upon algebra of constraints and, secondly, require extension of phase
space. Indeed, a transformation for a variable $q^{a}$ produced by any
generator $G$ in phase space reads

$\delta q^{a}=\\{q^{a},\,G\\}.$ (11)

So, to generate correct transformations for gauge variables the Poisson
brackets should be defined in extended phase space. Again, the dependence of
the algorithm on the algebra of constraints together with non-canonicity of
the transformations like (2.9) leads to the fact that the algorithm would work
only for a limited class of parametrizations. Thus, non-equivalence of
Hamiltonian formulations for different parametrizations, resulting in
different algebra of constraints, prevents from constructing the generator
which would produce correct transformations for any parametrizations. In the
next section we shall illustrate it making use of the algorithm [7], for a
simple model with finite number of degrees of freedom.

##  3\. The generator of gauge transformation: a simple example

Now we shall consider a closed isotropic cosmological model with the
Lagrangian

$L_{1}=-\frac{1}{2}\frac{a\dot{a}^{2}}{N}+\frac{1}{2}Na.$ (1)

This model is traditionally described in the ADM variables ($N$ is the lapse
function, $a$ is the scale factor). For our purpose, it is more convenient to
go to a new variable $\mu=N^{2}$ which corresponds to $g_{00}$. So the
Lagrangian is

$L_{2}=-\frac{1}{2}\frac{a\dot{a}^{2}}{\sqrt{\mu}}+\frac{1}{2}\sqrt{\mu}\,a.$
(2)

The canonical Hamiltonian constructed according to the rule
$H=p_{a}\dot{q}^{a}-L$, where $\\{p_{a},\;q^{a}\\}$ are pairs of variables
called canonical in the sense that all the velocities $\dot{q}^{a}$ can be
expressed through conjugate momenta, for our model is

$H_{C}=p\dot{a}-L_{2}=-\frac{1}{2}\frac{\sqrt{\mu}}{a}\;p^{2}-\frac{1}{2}\sqrt{\mu}\,a$
(3)

($p$ is the momentum conjugate to the scale factor). However, some authors
include into the form $p_{a}\dot{q}^{a}$ also gauge variables and their
momenta which are non-canonical variables in the above sense. Then we have the
so-called total Hamiltonian which for our model takes the form

$H_{T}=\pi\dot{\mu}+p\dot{a}-L_{2}=\pi\dot{\mu}-\frac{1}{2}\frac{\sqrt{\mu}}{a}\;p^{2}-\frac{1}{2}\sqrt{\mu}\,a$
(4)

($\pi$ is the momentum conjugate to the gauge variable $\mu$). Making use of
the total Hamiltonian implies a mixed formalism in which the Hamiltonian is
written in terms of canonical coordinates and momenta but as well of
velocities that cannot be expressed through the momenta. Nevertheless, this
very Hamiltonian plays the central role in the algorithm suggested in [7]
while the canonical Hamiltonian (3.3) will not lead to the correct result.

In [7] the generator of gauge transformations is sought in the form

$G=\sum\limits_{n}\theta_{\mu}^{(n)}G_{n}^{\mu},$ (5)

where $G_{n}^{\mu}$ are first class constraints, $\theta_{\mu}^{(n)}$ are the
$n$th order time derivatives of the gauge parameters $\theta_{\mu}$. In the
theory of gravity the variations of $g_{\mu\nu}$ involve first order
derivatives of gauge parameters, thus the generator is

$G=\theta_{\mu}G_{0}^{\mu}+\dot{\theta}_{\mu}G_{1}^{\mu}.$ (6)

$G_{n}^{\mu}$ satisfy the following conditions that were derived from the
requirement of invariance of motion equations under transformations in phase
space:

$G_{1}^{\mu}\quad{\rm are\;primary\;constraints};$ (7)
$G_{0}^{\mu}+\left\\{G_{1}^{\mu},\;H\right\\}\quad{\rm
are\;primary\;constraints};$ (8) $\left\\{G_{0}^{\mu},\;H\right\\}\quad{\rm
are\;primary\;constraints}.$ (9)

In our case $\pi=0$ is the only primary constraint of the model, so that
$G_{1}=\pi$. The secondary constraint is

$\dot{\pi}=\left\\{\pi,\;H_{T}\right\\}=-\frac{\partial
H_{T}}{\partial\mu}=\frac{1}{4}\frac{1}{a\sqrt{\mu}}\;p^{2}+\frac{1}{4}\frac{a}{\sqrt{\mu}}=T.$
(10)

The canonical Hamiltonian (3.3) appears to be proportional to the secondary
constraint $T$, $H_{C}=-2\mu T$.

The condition (3.8) becomes

$G_{0}+\left\\{\pi,\;H_{T}\right\\}=\alpha\pi;$ (11) $G_{0}=-T+\alpha\pi,$
(12)

$\alpha$ is a coefficient that can be found from the requirement (3.9):

$\left\\{G_{0},\;H_{T}\right\\}=\beta\pi;$ (13)
$\displaystyle\left\\{G_{0},\;H_{T}\right\\}$ $\displaystyle=$
$\displaystyle-\left\\{T,\;H_{T}\right\\}+\alpha\left\\{\pi,\;H_{T}\right\\}=-\left\\{T,\;\pi\dot{\mu}-2\mu
T\right\\}+\alpha T$ (14) $\displaystyle=$
$\displaystyle-\left\\{T,\;\pi\right\\}\dot{\mu}+\alpha
T=\frac{1}{2\mu}\;\dot{\mu}T+\alpha T;$
$\beta=0;\qquad\alpha=-\frac{1}{2\mu}\;\dot{\mu};$ (15)
$G_{0}=-\frac{1}{2\mu}\;\dot{\mu}\pi-T.$ (16)

The full generator $G$ (3.6) can be written as

$G=\left(-\frac{1}{2\mu}\;\dot{\mu}\pi-T\right)\theta+\pi\dot{\theta}.$ (17)

The transformation of the variable $\mu$ is

$\delta\mu=\left\\{\mu,\;G\right\\}=-\frac{1}{2\mu}\;\dot{\mu}\theta+\dot{\theta}.$
(18)

The same expression (up to the multiplier being equal to 2) can be obtained
from general transformations of the metric tensor,

$\delta
g_{\mu\nu}=\theta^{\lambda}\partial_{\lambda}g_{\mu\nu}+g_{\mu\lambda}\partial_{\nu}\theta^{\lambda}+g_{\nu\lambda}\partial_{\mu}\theta^{\lambda};$
(19) $\delta g_{00}=\dot{g}_{00}\theta^{0}+2g_{00}\dot{\theta}^{0},$ (20)

if one keeps in mind that $g_{00}=\mu$ and in the above formulas
$\theta=\theta_{0}=g_{00}\theta^{0}$.

Ir is easy to see that the correct expression (3.18) is entirely due to the
replacement of the canonical Hamiltonian (3.3) by the total Hamiltonian (3.4),
otherwise one would miss the contribution from the Poisson bracket
$\\{T,\;\pi\\}$ to the generator (3.17) (see the second line of (3.14)).

On the other hand, making use of the total Hamiltonian may not lead to a
correct result for another parametrization. Let us return to the Lagrangian
(3.1). Now the total Hamiltonian is

$H_{T}=\pi\dot{N}-\frac{1}{2}\frac{N}{a}\;p^{2}-\frac{1}{2}\;N\,a$ (21)

Again, $\pi$ is the momentum conjugate to the gauge variable $N$, and $\pi=0$
is the only primary constraint. The secondary constraint does not depend on
$N$ in this case:

$\dot{\pi}=\left\\{\pi,\;H_{T}\right\\}=-\frac{\partial H_{T}}{\partial
N}=\frac{1}{2a}\;p^{2}+\frac{1}{2}\;a=T,$ (22)

therefore, the Poisson bracket $\left\\{T,\;\pi\right\\}$ in (3.14) is equal
to zero, and one would obtain an incorrect expression for the generator,

$G=-T\theta+\pi\dot{\theta}.$ (23)

It cannot produce the correct variation of $N$, that reads

$\delta N=-\dot{N}\theta-N\dot{\theta}.$ (24)

As we can see, this algorithm fails to produce correct results for an
arbitrary parametrization. In the next section we shall construct Hamiltonian
dynamics in extended phase space and discuss its features and advantages.

##  4\. Extended phase space: the isotropic model

We shall consider the effective action including gauge and ghost sectors as it
appears in the path integral approach to gauge field theories,

$S=\int dt\left(L_{(grav)}+L_{(gauge)}+L_{(ghost)}\right)$ (1)

As was mentioned above, it is not enough just to extend pase space by
including formally gauge degrees of freedom in it. One should also introduce
missing velocities into the Lagrangian. It can be done by means of special
(differential) gauge conditions that actually extends the phase space and
enables one to avoid the mixed formalism. For our model (3.1) the equation
$N=f(a)$ determines in a general form a relation between the only gauge
variable $N$ and the scale factor $a$. The differential form of this relation
is

$\dot{N}=\frac{df}{da}\;\dot{a}.$ (2)

The ghost sector of the model reads

$L_{(ghost)}=\dot{\bar{\theta}}N\dot{\theta}+\dot{\bar{\theta}}\left(\dot{N}-\frac{df}{da}\;\dot{a}\right)\theta,$
(3)

so that

$\displaystyle L$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\frac{a\dot{a}^{2}}{N}+\frac{1}{2}Na+\lambda\left(\dot{N}-\frac{df}{da}\;\dot{a}\right)+\dot{\bar{\theta}}\left(\dot{N}-\frac{df}{da}\;\dot{a}\right)\theta+\dot{\bar{\theta}}N\dot{\theta}=$
(4) $\displaystyle=$
$\displaystyle-\frac{1}{2}\frac{a\dot{a}^{2}}{N}+\frac{1}{2}Na+\pi\left(\dot{N}-\frac{df}{da}\;\dot{a}\right)+\dot{\bar{\theta}}N\dot{\theta}.$

The conjugate momenta are:

$\pi=\lambda+\dot{\bar{\theta}}\theta;\quad
p=-\frac{a\dot{a}}{N}-\pi\frac{df}{da};\quad\bar{\cal
P}=N\dot{\bar{\theta}};\quad{\cal P}=N\dot{\theta}.$ (5)

Let us now go to a new variable

$N=v(\tilde{N},a).$ (6)

At the same time, the rest variables are unchanged:

$a=\tilde{a};\quad\theta=\tilde{\theta};\quad\bar{\theta}=\tilde{\bar{\theta}}.$
(7)

It is the analog of the transformation from the original gravitational
variables $g_{\mu\nu}$ to the ADM variables. Indeed, in the both cases only
gauge variables are transformed while the rest variables remain unchanged.
After the change (4.6) the Lagrangian is written as (below we shall omit the
tilde over $a$ and ghost variables which remain unchanged)

$L=-\frac{1}{2}\;\frac{a\dot{a}^{2}}{v(\tilde{N},a)}+\frac{1}{2}\;v(\tilde{N},a)\;a+\pi\left(\frac{\partial
v}{\partial\tilde{N}}\;\dot{\tilde{N}}+\frac{\partial v}{\partial
a}\;\dot{a}-\frac{df}{da}\;\dot{a}\right)+v(\tilde{N},a)\;\dot{\bar{\theta}}\dot{\theta}.$
(8)

The new momenta are:

$\tilde{\pi}=\pi\frac{\partial
v}{\partial\tilde{N}};\qquad\tilde{p}=-\frac{a\dot{a}}{v(\tilde{N},a)}+\pi\frac{\partial
v}{\partial a}-\pi\frac{df}{da}=p+\pi\frac{\partial v}{\partial a};$ (9)
$\tilde{\bar{\cal P}}=v(\tilde{N},a)\;\dot{\bar{\theta}}=\bar{\cal
P};\qquad\tilde{\cal P}=v(\tilde{N},a)\;\dot{\theta}={\cal P}.$ (10)

It is easy to demonstrate that the transformations (4.6), (4.7), (4.9), (4.10)
are canonical in extended phase space. The generating function will depend on
new coordinates and old momenta,

$\Phi\left(\tilde{N},\;\tilde{a},\;\tilde{\bar{\theta}},\;\tilde{\theta},\;\pi,\;p,\;\bar{\cal
P},\;{\cal P}\right)=-\pi\,v(\tilde{N},\tilde{a})-p\,\tilde{a}-\bar{\cal
P}\,\tilde{\theta}-\tilde{\bar{\theta}}\,{\cal P}.$ (11)

One can see that the generating function has the same form as in (2.5). The
relations

$N=-\frac{\partial\Phi}{\partial\pi};\qquad a=-\frac{\partial\Phi}{\partial
p};\qquad\tilde{\pi}=-\frac{\partial\Phi}{\partial\tilde{N}};\qquad\tilde{p}=-\frac{\partial\Phi}{\partial\tilde{a}};$
(12) $\theta=-\frac{\partial\Phi}{\partial\bar{\cal
P}\vphantom{\sqrt{N}}};\qquad\bar{\theta}=-\frac{\partial\Phi}{\partial{\cal
P}};\qquad\tilde{\cal
P}=-\frac{\partial\Phi}{\partial\tilde{\bar{\theta}}};\qquad\tilde{\bar{\cal
P}}=-\frac{\partial\Phi}{\partial\tilde{\theta}}$ (13)

give exactly the transformation (4.6), (4.7), (4.9), (4.10). On the other
hand, one can check that Poisson brackets among all phase variables maintain
their canonical form.

Now we are going to write down equations of motion in extended phase space.
Firstly, we rewrite the Lagrangian (4.8) through the momentum $\tilde{\pi}$.

$\displaystyle L$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\;\frac{a\dot{a}^{2}}{v(\tilde{N},a)}+\frac{1}{2}\;v(\tilde{N},a)\;a$
(14) $\displaystyle+$
$\displaystyle\tilde{\pi}\left[\dot{\tilde{N}}+\left(\frac{\partial
v}{\partial\tilde{N}}\right)^{-1}\frac{\partial v}{\partial
a}\;\dot{a}-\left(\frac{\partial
v}{\partial\tilde{N}}\right)^{-1}\frac{df}{da}\;\dot{a}\right]+v(\tilde{N},a)\;\dot{\bar{\theta}}\dot{\theta}.$

The variation of (4.14) gives, accordingly, the equation of motion (4.15), the
constraint (4.16), the gauge condition (4.17) and the ghost equations (4.18) –
(4.19):

$\displaystyle\frac{a\ddot{a}}{v(\tilde{N},a)}$ $\displaystyle+$
$\displaystyle\frac{1}{2}\;\frac{\dot{a}^{2}}{v(\tilde{N},a)}-\frac{1}{2}\;\frac{a\dot{a}^{2}}{v^{2}(\tilde{N},a)}\;\frac{\partial
v}{\partial a}-\frac{a\dot{a}}{v^{2}(\tilde{N},a)}\;\frac{\partial
v}{\partial\tilde{N}}\dot{\tilde{N}}$ (15) $\displaystyle+$
$\displaystyle\frac{1}{2}\;\frac{\partial v}{\partial
a}\;a+\frac{1}{2}v(\tilde{N},a)-\dot{\tilde{\pi}}\left(\frac{\partial
v}{\partial\tilde{N}}\right)^{-1}\frac{\partial v}{\partial
a}+\dot{\tilde{\pi}}\left(\frac{\partial
v}{\partial\tilde{N}}\right)^{-1}\frac{df}{da}$ $\displaystyle+$
$\displaystyle\tilde{\pi}\left(\frac{\partial
v}{\partial\tilde{N}}\right)^{-2}\frac{\partial^{2}v}{\partial\tilde{N}^{2}}\;\frac{\partial
v}{\partial a}\;\dot{\tilde{N}}-\tilde{\pi}\left(\frac{\partial
v}{\partial\tilde{N}}\right)^{-1}\frac{\partial^{2}v}{\partial\tilde{N}\partial
a}\;\dot{\tilde{N}}$ $\displaystyle-$
$\displaystyle\tilde{\pi}\left(\frac{\partial
v}{\partial\tilde{N}}\right)^{-2}\frac{\partial^{2}v}{\partial\tilde{N}^{2}}\;\frac{df}{da}\;\dot{\tilde{N}}+\frac{\partial
v}{\partial a}\;\dot{\bar{\theta}}\dot{\theta}=0;$
$\displaystyle\frac{1}{2}\;\frac{a\dot{a}^{2}}{v^{2}(\tilde{N},a)}\;\frac{\partial
v}{\partial\tilde{N}}$ $\displaystyle+$
$\displaystyle\frac{1}{2}\;\frac{\partial
v}{\partial\tilde{N}}\;a-\dot{\tilde{\pi}}-\tilde{\pi}\left(\frac{\partial
v}{\partial\tilde{N}}\right)^{-2}\frac{\partial^{2}v}{\partial\tilde{N}^{2}}\;\frac{\partial
v}{\partial a}\;\dot{a}+\tilde{\pi}\left(\frac{\partial
v}{\partial\tilde{N}}\right)^{-1}\frac{\partial^{2}v}{\partial\tilde{N}\partial
a}\;\dot{a}$ (16) $\displaystyle+$
$\displaystyle\tilde{\pi}\left(\frac{\partial
v}{\partial\tilde{N}}\right)^{-2}\frac{\partial^{2}v}{\partial\tilde{N}^{2}}\;\frac{df}{da}\;\dot{a}+\frac{\partial
v}{\partial\tilde{N}}\;\dot{\bar{\theta}}\dot{\theta}=0;$ $\frac{\partial
v}{\partial\tilde{N}}\;\dot{\tilde{N}}+\frac{\partial v}{\partial
a}\;\dot{a}-\frac{df}{da}\;\dot{a}=0;$ (17)
$v(\tilde{N},\;a)\;\ddot{\theta}+\frac{\partial
v}{\partial\tilde{N}}\;\dot{\tilde{N}}\dot{\theta}+\frac{\partial v}{\partial
a}\;\dot{a}\dot{\theta}=0;$ (18)
$v(\tilde{N},\;a)\;\ddot{\bar{\theta}}+\frac{\partial
v}{\partial\tilde{N}}\;\dot{\tilde{N}}\dot{\bar{\theta}}+\frac{\partial
v}{\partial a}\;\dot{a}\dot{\bar{\theta}}=0.$ (19)

The Hamiltonian in extended phase space looks like

$\displaystyle H$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\;\frac{v(\tilde{N},a)}{a}\left[\tilde{p}^{2}+2\tilde{p}\tilde{\pi}\left(\frac{\partial
v}{\partial\tilde{N}}\right)^{-1}\frac{df}{da}+\tilde{\pi}^{2}\left(\frac{\partial
v}{\partial\tilde{N}}\right)^{-2}\left(\frac{df}{da}\right)^{2}\right.$ (20)
$\displaystyle-$ $\displaystyle\left.2\tilde{p}\tilde{\pi}\left(\frac{\partial
v}{\partial\tilde{N}}\right)^{-1}\frac{\partial v}{\partial
a}-2\tilde{\pi}^{2}\left(\frac{\partial
v}{\partial\tilde{N}}\right)^{-2}\frac{\partial v}{\partial
a}\;\frac{df}{da}+\tilde{\pi}^{2}\left(\frac{\partial
v}{\partial\tilde{N}}\right)^{-2}\left(\frac{\partial v}{\partial
a}\right)^{2}\right]$ $\displaystyle-$
$\displaystyle\frac{1}{2}\;v(\tilde{N},a)\;a+\frac{1}{v(\tilde{N},a)}\;\bar{\cal
P}{\cal P}.$

The Hamiltonian equations in extended phase space are:

$\displaystyle\dot{\tilde{p}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left[\frac{1}{a}\frac{\partial v}{\partial
a}-\frac{v(\tilde{N},a)}{a^{2}}\right]\left[\tilde{p}+\tilde{\pi}\left(\frac{\partial
v}{\partial\tilde{N}}\right)^{-1}\frac{df}{da}-\tilde{\pi}\left(\frac{\partial
v}{\partial\tilde{N}}\right)^{-1}\frac{\partial v}{\partial a}\right]^{2}$
(21) $\displaystyle-$
$\displaystyle\frac{v(\tilde{N},a)}{a}\left[\tilde{\pi}\left(\frac{\partial
v}{\partial\tilde{N}}\right)^{-2}\frac{\partial^{2}v}{\partial\tilde{N}\partial
a}\;\frac{df}{da}-\tilde{\pi}\left(\frac{\partial
v}{\partial\tilde{N}}\right)^{-1}\frac{d^{2}f}{da^{2}}\right.$
$\displaystyle-$ $\displaystyle\left.\tilde{\pi}\left(\frac{\partial
v}{\partial\tilde{N}}\right)^{-2}\frac{\partial^{2}v}{\partial\tilde{N}\partial
a}\;\frac{\partial v}{\partial a}+\tilde{\pi}\left(\frac{\partial
v}{\partial\tilde{N}}\right)^{-1}\frac{\partial^{2}v}{\partial a^{2}}\right]$
$\displaystyle\times$
$\displaystyle\left[\tilde{p}+\tilde{\pi}\left(\frac{\partial
v}{\partial\tilde{N}}\right)^{-1}\frac{df}{da}-\tilde{\pi}\left(\frac{\partial
v}{\partial\tilde{N}}\right)^{-1}\frac{\partial v}{\partial a}\right]$
$\displaystyle+$ $\displaystyle\frac{1}{2}\;\frac{\partial v}{\partial
a}\;a+\frac{1}{2}\;v(\tilde{N},a)+\frac{1}{v^{2}(\tilde{N},a)}\;\bar{\cal
P}{\cal P};$ $\displaystyle\dot{a}$ $\displaystyle=$
$\displaystyle-\frac{v(\tilde{N},a)}{a}\left[\tilde{p}+\tilde{\pi}\left(\frac{\partial
v}{\partial\tilde{N}}\right)^{-1}\frac{df}{da}-\tilde{\pi}\left(\frac{\partial
v}{\partial\tilde{N}}\right)^{-1}\frac{\partial v}{\partial a}\right];$ (22)
$\displaystyle\dot{\tilde{\pi}}$ $\displaystyle=$
$\displaystyle\frac{1}{2a}\;\frac{\partial
v}{\partial\tilde{N}}\left[\tilde{p}+\tilde{\pi}\left(\frac{\partial
v}{\partial\tilde{N}}\right)^{-1}\frac{df}{da}-\tilde{\pi}\left(\frac{\partial
v}{\partial\tilde{N}}\right)^{-1}\frac{\partial v}{\partial a}\right]^{2}$
(23) $\displaystyle-$
$\displaystyle\frac{v(\tilde{N},a)}{a}\left[\tilde{\pi}\left(\frac{\partial
v}{\partial\tilde{N}}\right)^{-2}\frac{\partial^{2}v}{\partial\tilde{N}^{2}}\;\frac{df}{da}-\tilde{\pi}\left(\frac{\partial
v}{\partial\tilde{N}}\right)^{-2}\frac{\partial^{2}v}{\partial\tilde{N}^{2}}\;\frac{\partial
v}{\partial a}\right.$ $\displaystyle+$
$\displaystyle\left.\tilde{\pi}\left(\frac{\partial
v}{\partial\tilde{N}}\right)^{-1}\frac{\partial^{2}v}{\partial\tilde{N}\partial
a}\right]\left[\tilde{p}+\tilde{\pi}\left(\frac{\partial
v}{\partial\tilde{N}}\right)^{-1}\frac{df}{da}-\tilde{\pi}\left(\frac{\partial
v}{\partial\tilde{N}}\right)^{-1}\frac{\partial v}{\partial a}\right]$
$\displaystyle+$ $\displaystyle\frac{1}{2}\;\frac{\partial
v}{\partial\tilde{N}}\;a+\frac{1}{v^{2}(\tilde{N},a)}\;\frac{\partial
v}{\partial\tilde{N}}\;\bar{\cal P}{\cal P};$ $\displaystyle\dot{\tilde{N}}$
$\displaystyle=$
$\displaystyle-\frac{v(\tilde{N},a)}{a}\left[\left(\frac{\partial
v}{\partial\tilde{N}}\right)^{-1}\frac{df}{da}-\left(\frac{\partial
v}{\partial\tilde{N}}\right)^{-1}\frac{\partial v}{\partial a}\right]$ (24)
$\displaystyle\times$
$\displaystyle\left[\tilde{p}+\tilde{\pi}\left(\frac{\partial
v}{\partial\tilde{N}}\right)^{-1}\frac{df}{da}-\tilde{\pi}\left(\frac{\partial
v}{\partial\tilde{N}}\right)^{-1}\frac{\partial v}{\partial a}\right];$
$\displaystyle\dot{\bar{\cal P}}$ $\displaystyle=$ $\displaystyle 0;$ (25)
$\displaystyle\dot{\theta}$ $\displaystyle=$
$\displaystyle\frac{1}{v(\tilde{N},a)}\;{\cal P};$ (26)
$\displaystyle\dot{\cal P}$ $\displaystyle=$ $\displaystyle 0;$ (27)
$\displaystyle\dot{\bar{\theta}}$ $\displaystyle=$
$\displaystyle\frac{1}{v(\tilde{N},a)}\;\bar{\cal P}.$ (28)

One can check that the Hamiltonian equations (4.21) – (4.28) are completely
equivalent to the Lagrangian equations (4.15) – (4.19), the constraint (4.23)
and the gauge condition (4.24) being true Hamiltonian equations.

The Hamiltonian equations (4.21) – (4.28) in extended phase space, as well as
the equations (4.15) – (4.19), include gauge-dependent terms. In this
connection one can object that the equations are not equivalent to the
original Einstein equation, which are known to be gauge-invariant. However, we
remember that any solution to the gauge-invariant Einstein equation is
determined up to arbitrary functions which have to be fix by a choice of a
reference frame (a state of the observer). It is usually done on the final
stage of solving the Einstein equations. It is important that one cannot avoid
fixing a reference frame to obtain a final form of the solution. By varying
the gauged action (4.1), in fact, we deal with a generalized mathematical
problem, its generalization has come from the development of quantum field
theory.

In the case of the extended set of equations (4.21) – (4.28) (or, (4.15) –
(4.19)) one can keep the function $f(a)$ non-fixed up to the final stage of
their resolution. Further, under the conditions $\bar{\pi}=0$, $\theta=0$,
$\bar{\theta}=0$ all gauge-dependent terms are excluded, and the extended set
of equations is reduced to gauge-invariant equations, therefore, any solution
of the Einstein equations can be found among solutions of the extended set.
Solutions with non-trivial $\bar{\pi}$, $\theta$, $\bar{\theta}$ should be
considered and physically interpreted separately.

One can also reveal that there exists a quantity conserved if the Hamiltonian
(or, equivalently, Lagrangian) equations hold. It plays the role of the BRST
generator for our model:

$\Omega=-H\theta-\left(\frac{\partial
v}{\partial\tilde{N}}\right)^{-1}\tilde{\pi}{\cal P}.$ (29)

It generates correct transformations for the variables $a$, $\theta$,
$\bar{\theta}$ and for any gauge variable $\tilde{N}$ given by the relation
(4.6),

$\delta\tilde{N}=-\frac{\partial
H}{\partial\tilde{\pi}}\;\theta-\left(\frac{\partial
v}{\partial\tilde{N}}\right)^{-1}{\cal
P}=-\dot{\tilde{N}}\theta-\left(\frac{\partial
v}{\partial\tilde{N}}\right)^{-1}v(\tilde{N},a)\;\dot{\theta}.$ (30)

In particular, for the original variable $N$ one gets the transformation
(3.24).

##  5\. The canonicity of transformations in extended phase space
for the full gravitational theory

In this section we shall demonstrate for the full gravitational theory that
different parametrizations from a wide enough class (2.9) are related by
canonical transformations. Again, we shall start from the gauged action

$S=\int d^{4}x\left({\cal L}_{(grav)}+{\cal L}_{(gauge)}+{\cal
L}_{(ghost)}\right)$ (1)

We shall use a gauge condition in a general form, $f^{\mu}(g_{\nu\lambda})=0$.
The differential form of this gauge condition introduces the missing
velocities and actually extends phase space,

$\frac{d}{dt}f^{\mu}(g_{\nu\lambda})=0,\qquad\frac{\partial f^{\mu}}{\partial
g_{00}}\dot{g}_{00}+2\frac{\partial f^{\mu}}{\partial
g_{0i}}\dot{g}_{0i}+\frac{\partial f^{\mu}}{\partial g_{ij}}\dot{g}_{ij}=0.$
(2)

Then,

${\cal L}_{(gauge)}=\lambda_{\mu}\left(\frac{\partial f^{\mu}}{\partial
g_{00}}\dot{g}_{00}+2\frac{\partial f^{\mu}}{\partial
g_{0i}}\dot{g}_{0i}+\frac{\partial f^{\mu}}{\partial
g_{ij}}\dot{g}_{ij}\right).$ (3)

Taking into account the gauge transformations,

$\delta
g_{\mu\nu}=\partial_{\lambda}g_{\mu\nu}\theta^{\lambda}+g_{\mu\lambda}\partial_{\nu}\theta^{\lambda}+g_{\nu\lambda}\partial_{\mu}\theta^{\lambda},$
(4)

one can write the ghost sector:

${\cal L}_{(ghost)}=\bar{\theta}_{\mu}\frac{d}{dt}\left[\frac{\partial
f^{\mu}}{\partial
g_{\nu\lambda}}\left(\partial_{\rho}g_{\nu\lambda}\theta^{\rho}+g_{\lambda\rho}\partial_{\nu}\theta^{\rho}+g_{\nu\rho}\partial_{\lambda}\theta^{\rho}\right)\right].$
(5)

It is convenient to write down the action (5.1), (5.3), (5.5) in the form

$\displaystyle S$ $\displaystyle=$ $\displaystyle\int d^{4}x\left[{\cal
L}_{(grav)}+\Lambda_{\mu}\left(\frac{\partial f^{\mu}}{\partial
g_{00}}\dot{g}_{00}+2\frac{\partial f^{\mu}}{\partial
g_{0i}}\dot{g}_{0i}+\frac{\partial f^{\mu}}{\partial
g_{ij}}\dot{g}_{ij}\right)\right.$ (6) $\displaystyle-$
$\displaystyle\dot{\bar{\theta_{\mu}}}\left(\frac{\partial f^{\mu}}{\partial
g_{00}}\left(\partial_{i}g_{00}\theta^{i}+2g_{0\nu}\dot{\theta}^{\nu}\right)+2\frac{\partial
f^{\mu}}{\partial
g_{0i}}\left(\partial_{j}g_{0i}\theta^{j}+g_{0\nu}\partial_{i}\theta^{\nu}+g_{i\nu}\dot{\theta}^{\nu}\right)\right.$
$\displaystyle+$ $\displaystyle\left.\left.\frac{\partial f^{\mu}}{\partial
g_{ij}}\left(\partial_{k}g_{ij}\theta^{k}+g_{i\nu}\partial_{j}\theta^{\nu}+g_{j\nu}\partial_{i}\theta^{\nu}\right)\right)\right].$

Here $\Lambda_{\mu}=\lambda_{\mu}-\dot{\bar{\theta_{\mu}}}\theta^{0}$. One can
see that the generalized velocities enter into the bracket multiplied by
$\Lambda_{\mu}$, in addition to the gravitational part ${\cal L}_{(grav)}$.
This very circumstance will ensure the canonicity of the transformation to new
variables.

Our goal now is to introduce new variables by

$g_{0\mu}=v_{\mu}\left(N_{\nu},g_{ij}\right).\qquad
g_{ij}=\gamma_{ij};\qquad\theta^{\mu}=\tilde{\theta}^{\mu};\qquad\bar{\theta}_{\mu}=\tilde{\bar{\theta}_{\mu}}.$
(7)

This is the inverse transformation for (2.9) and concerns only $g_{0\mu}$
metric components. After the transformation the action will read

$\displaystyle S$ $\displaystyle=$ $\displaystyle\int d^{4}x\left[{\cal
L^{\prime}}_{(grav)}+\Lambda_{\mu}\left(\frac{\partial f^{\mu}}{\partial
g_{00}}\;\frac{\partial v_{0}}{\partial N_{\nu}}\;\dot{N}_{\nu}+\frac{\partial
f^{\mu}}{\partial g_{00}}\;\frac{\partial v_{0}}{\partial
g_{ij}}\;\dot{g}_{ij}\right.\right.$ (8) $\displaystyle+$
$\displaystyle\left.2\;\frac{\partial f^{\mu}}{\partial
g_{0i}}\;\frac{\partial v_{i}}{\partial
N_{\nu}}\;\dot{N}_{\nu}+2\;\frac{\partial f^{\mu}}{\partial
g_{0k}}\;\frac{\partial v_{k}}{\partial g_{ij}}\;\dot{g}_{ij}+\frac{\partial
f^{\mu}}{\partial g_{ij}}\;\dot{g}_{ij}\right)$ $\displaystyle-$
$\displaystyle\dot{\bar{\theta_{\mu}}}\left(\frac{\partial f^{\mu}}{\partial
g_{00}}\;\frac{\partial v_{0}}{\partial
N_{\nu}}\;\partial_{i}N_{\nu}\theta^{i}+\frac{\partial f^{\mu}}{\partial
g_{00}}\;\frac{\partial v_{0}}{\partial
g_{ij}}\;\partial_{k}g_{ij}\theta^{k}+2\;\frac{\partial f^{\mu}}{\partial
g_{00}}\;v_{\nu}(N_{\lambda},g_{ij})\;\dot{\theta}^{\nu}\right.$
$\displaystyle+$ $\displaystyle 2\;\frac{\partial f^{\mu}}{\partial
g_{0i}}\;\frac{\partial v_{i}}{\partial
N_{\nu}}\;\partial_{j}N_{\nu}\theta^{j}+2\;\frac{\partial f^{\mu}}{\partial
g_{0i}}\;\frac{\partial v_{i}}{\partial g_{jk}}\;\partial_{l}g_{jk}\theta^{l}$
$\displaystyle+$ $\displaystyle 2\;\frac{\partial f^{\mu}}{\partial
g_{0i}}\left[v_{\nu}(N_{\lambda},g_{jk})\;\partial_{i}\theta^{\nu}+v_{i}(N_{\lambda},g_{jk})\;\dot{\theta}^{0}+g_{ij}\dot{\theta}^{j}\right]$
$\displaystyle+$ $\displaystyle\frac{\partial f^{\mu}}{\partial
g_{ij}}\left[\partial_{k}g_{ij}\theta^{k}+v_{i}(N_{\lambda},g_{kl})\;\partial_{j}\theta^{0}+g_{ik}\partial_{j}\theta^{k}\right.$
$\displaystyle+$
$\displaystyle\left.\left.\left.v_{j}(N_{\lambda},g_{kl})\;\partial_{i}\theta^{0}+g_{jk}\partial_{i}\theta^{k}\right]\right)\right]$

We can write down the “old” momenta,

$\displaystyle\pi^{ij}$ $\displaystyle=$ $\displaystyle\frac{\partial{\cal
L}_{(grav)}}{\partial\dot{g}_{ij}}+\Lambda_{\mu}\;\frac{\partial
f^{\mu}}{\partial g_{ij}};$ (9) $\displaystyle\pi^{0}$ $\displaystyle=$
$\displaystyle\frac{\partial{\cal
L}_{(grav)}}{\partial\dot{g}_{00}}+\Lambda_{\mu}\;\frac{\partial
f^{\mu}}{\partial g_{00}};$ (10) $\displaystyle\pi^{i}$ $\displaystyle=$
$\displaystyle\frac{\partial{\cal
L}_{(grav)}}{\partial\dot{g}_{0i}}+2\Lambda_{\mu}\;\frac{\partial
f^{\mu}}{\partial g_{0i}},$ (11)

and the “new” momenta are:

$\displaystyle\Pi^{ij}$ $\displaystyle=$ $\displaystyle\frac{\partial{\cal
L^{\prime}}_{(grav)}}{\partial\dot{g}_{ij}}+\Lambda_{\mu}\left(\frac{\partial
f^{\mu}}{\partial g_{00}}\;\frac{\partial v_{0}}{\partial
g_{ij}}+2\;\frac{\partial f^{\mu}}{\partial g_{0k}}\;\frac{\partial
v_{k}}{\partial g_{ij}}+\frac{\partial f^{\mu}}{\partial g_{ij}}\right);$ (12)
$\displaystyle\Pi^{0}$ $\displaystyle=$ $\displaystyle\frac{\partial{\cal
L^{\prime}}_{(grav)}}{\partial\dot{N}_{0}}+\Lambda_{\mu}\left(\frac{\partial
f^{\mu}}{\partial g_{00}}\;\frac{\partial v_{0}}{\partial
N_{0}}+2\;\frac{\partial f^{\mu}}{\partial g_{0i}}\;\frac{\partial
v_{i}}{\partial N_{0}}\right);$ (13) $\displaystyle\Pi^{i}$ $\displaystyle=$
$\displaystyle\frac{\partial{\cal
L^{\prime}}_{(grav)}}{\partial\dot{N}_{i}}+\Lambda_{\mu}\left(\frac{\partial
f^{\mu}}{\partial g_{00}}\;\frac{\partial v_{0}}{\partial
N_{i}}+2\;\frac{\partial f^{\mu}}{\partial g_{0j}}\;\frac{\partial
v_{j}}{\partial N_{i}}\right).$ (14)

The relations between the “old” and “new” momenta:

$\displaystyle\Pi^{ij}$ $\displaystyle=$
$\displaystyle\pi^{ij}+\left(\pi^{\mu}-\frac{\partial{\cal
L}_{(grav)}}{\partial\dot{g}_{0\mu}}\right)\frac{\partial v_{\mu}}{\partial
g_{ij}};$ (15) $\displaystyle\Pi^{\mu}$ $\displaystyle=$
$\displaystyle\frac{\partial{\cal
L^{\prime}}_{(grav)}}{\partial\dot{N}_{\mu}}+\left(\pi^{\nu}-\frac{\partial{\cal
L}_{(grav)}}{\partial\dot{g}_{0\nu}}\right)\frac{\partial v_{\nu}}{\partial
N_{\mu}}.$ (16)

It is easy to check that the momenta conjugate to ghosts remain unchanged,
$\tilde{\cal P}^{\mu}={\cal P}^{\mu}$, $\tilde{\bar{\cal P}}_{\mu}=\bar{\cal
P}_{\mu}$.

As any Lagrangian is determined up to total derivatives, the gravitational
Lagrangian density ${\cal L}_{(grav)}$ can be modified in such a way for the
primary constraints to take the form $\pi^{\mu}=0$, where $\pi^{\mu}$ are the
momenta conjugate to gauge variables $g_{0\mu}$. This change of the Lagrangian
density does not affect the equation of motion. It was made by Dirac [3] to
simplify the calculations. A similar change of the Lagrangian density by
omitting a divergence and a total time derivative was made also in the ADM
paper [9]. Therefore, one can put

$\frac{\partial{\cal
L}_{(grav)}}{\partial\dot{g}_{0\mu}}=0,\qquad\frac{\partial{\cal
L^{\prime}}_{(grav)}}{\partial\dot{N}_{\mu}}=0.$ (17)

Then, the relations (5.15) – (5.16) would become simpler and take the form

$\Pi^{ij}=\pi^{ij}+\pi^{\mu}\frac{\partial v_{\mu}}{\partial
g_{ij}};\qquad\Pi^{\mu}=\pi^{\nu}\frac{\partial v_{\nu}}{\partial N_{\mu}}.$
(18)

It is easy to demonstrate that the transformations (2.9), (5.18) are canonical
in extended phase space. The generating function again depends on new
coordinates and old momenta and has the same form as for a non-constrained
system (see (2.5), compare also with (4.11)),

$\Phi\left(N_{\mu},\;\gamma_{ij},\;\tilde{\theta}^{\mu},\;\tilde{\bar{\theta}}_{\mu},\;\pi^{\mu},\;\pi^{ij},\;\bar{\cal
P}_{\mu},\;{\cal
P}^{\mu}\right)=-\pi^{\mu}v_{\mu}(N_{\nu},\gamma_{ij})-\pi^{ij}\gamma_{ij}-\bar{\cal
P}_{\mu}\tilde{\theta}^{\mu}-\tilde{\bar{\theta}}_{\mu}{\cal P}^{\mu}.$ (19)

Then the following relations take place

$g_{0\mu}=-\frac{\partial\Phi}{\partial\pi^{\mu}};\qquad
g_{ij}=-\frac{\partial\Phi}{\partial\pi^{ij}};\qquad\theta^{\mu}=-\frac{\partial\Phi}{\partial\bar{\cal
P}\vphantom{\sqrt{N}}_{\mu}};\qquad\bar{\theta}_{\mu}=-\frac{\partial\Phi}{\partial{\cal
P}^{\mu}};$ (20) $\Pi^{\mu}=-\frac{\partial\Phi}{\partial
N_{\mu}};\qquad\Pi^{ij}=-\frac{\partial\Phi}{\partial\gamma_{ij}};\qquad\tilde{\bar{\cal
P}}_{\mu}=-\frac{\partial\Phi}{\partial\tilde{\theta}^{\mu}};\qquad\tilde{\cal
P}^{\mu}=-\frac{\partial\Phi}{\partial\tilde{\bar{\theta}}\vphantom{\sqrt{N}}_{\mu}},$
(21)

that give exactly the transformations

$g_{0\mu}=v_{\mu}(N_{\nu},\gamma_{ij});\qquad
g_{ij}=\gamma_{ij};\qquad\qquad\qquad\qquad\theta^{\mu}=\tilde{\theta}^{\mu};\qquad\bar{\theta}_{\mu}=\tilde{\bar{\theta}}_{\mu};$
(22) $\Pi^{\mu}=\pi^{\nu}\frac{\partial v_{\nu}}{\partial
N_{\mu}};\qquad\qquad\Pi^{ij}=\pi^{ij}+\pi^{\mu}\frac{\partial
v_{\mu}}{\partial g_{ij}};\qquad\quad\tilde{\bar{\cal P}}_{\mu}=\bar{\cal
P}_{\mu};\qquad\tilde{\cal P}^{\mu}={\cal P}^{\mu}.$ (23)

We can now check if the Poisson brackets maintain their form. Differentiating
the first relation in (2.9) with respect to $g_{ij}$ one gets

$\frac{\partial V_{\mu}}{\partial g_{ij}}+\frac{\partial V_{\mu}}{\partial
g_{0\lambda}}\frac{\partial v_{\lambda}}{\partial g_{ij}}=0,$ (24)

Similarly, differentiating the same relation with respect to $N_{\nu}$ gives

$\delta_{\mu}^{\nu}-\frac{\partial V_{\mu}}{\partial
g_{0\lambda}}\frac{\partial v_{\lambda}}{\partial N_{\nu}}=0.$ (25)

Making use of (5.24), (5.25), it is not difficult to calculate the Poisson
brackets. So, we can recalculate the bracket (2.10) to see that it will be
zero in our extended phase space formalism,

$\displaystyle\left.\\{N_{\mu},\,\Pi^{ij}\\}\right|_{g_{\nu\lambda},p^{\rho\sigma}}$
$\displaystyle=$ $\displaystyle\frac{\partial N_{\mu}}{\partial
g_{0\rho}}\;\frac{\partial\Pi^{ij}}{\partial\pi^{\rho}}+\frac{\partial
N_{\mu}}{\partial
g_{kl}}\;\frac{\partial\Pi^{ij}}{\partial\pi^{kl}}=\left\\{V_{\mu}(g_{0\nu},g_{kl}),\;\pi^{ij}+\pi^{\lambda}\frac{\partial
v_{\lambda}}{\partial g_{ij}}\right\\}$ (26) $\displaystyle=$
$\displaystyle\frac{\partial V_{\mu}}{\partial g_{0\rho}}\;\frac{\partial
v_{\lambda}}{\partial g_{ij}}\delta_{\rho}^{\lambda}+\frac{\partial
V_{\mu}}{\partial
g_{kl}}\;\frac{1}{2}\left(\delta_{k}^{i}\delta_{l}^{j}+\delta_{l}^{j}\delta_{k}^{i}\right)=\frac{\partial
V_{\mu}}{\partial g_{0\lambda}}\;\frac{\partial v_{\lambda}}{\partial
g_{ij}}+\frac{\partial V_{\mu}}{\partial g_{ij}}=0.$

To give another example, let us check the following bracket:

$\left.\\{N_{\mu},\,\Pi^{\nu}\\}\right|_{g_{\lambda\rho},p^{\sigma\tau}}=\frac{\partial
N_{\mu}}{\partial
g_{0\rho}}\;\frac{\partial\Pi^{\nu}}{\partial\pi^{\rho}}=\left\\{V_{\mu}(g_{0\rho},g_{ij}),\;\pi^{\lambda}\frac{\partial
v_{\lambda}}{\partial N_{\nu}}\right\\}=\frac{\partial V_{\mu}}{\partial
g_{0\rho}}\;\frac{\partial v_{\rho}}{\partial N_{\nu}}=\delta_{\mu}^{\nu}.$
(27)

The rest of the brackets can be checked by analogy. This completes the proof
of canonicity of the transformation (2.9) for the full gravitational theory.

##  6\. Discussion

A starting point for the present investigation was the paper [11] and the
statement made by its authors that components of metric tensor and the ADM
variables are not related by a canonical transformation. However, it is
misunderstanding to pose the question about canonicity of the transformation
(2.8) which involves, from the viewpoint of the Dirac approach, non-canonical
variables. Let us remind that Dirac himself consider these variables,
$g_{0\mu}$ (along with the zero component of vector potential of
electromagnetic field $A_{0}$) as playing the role of Lagrange multipliers
while the phase space in his approach includes pairs of generalized
coordinates and momenta for which corresponding velocities can be expressed
through the momenta.

We should remember also that the Einstein equations were originally formulated
in Lagrangian formalism. Dirac’s Hamiltonian formulation for gravity is
equivalent to Einstein’s formulation at the level of equations. It means that
Hamiltonian equations for canonical variables (in Dirac’s sense) are
equivalent to the ($ij$) Einstein equations, and the gravitational constraints
are equivalent to the ($0\mu$) Einstein equations. On the other hand, it
implies that a group of transformations in Hamiltonian formalism must involve
the full group of gauge transformations of the original theory. However, in
the limits of the Dirac approach we fail to construct a generator that would
produce correct transformations for all variables. We inevitably have to
modify the Dirac scheme, and attempts to do it were presented yet in [7, 8].
Therefore, we cannot consider the Dirac approach as fundamental and undoubted.

The ADM formulation of Hamiltonian dynamics for gravity is, first of all, the
choice of parametrization, which is preferable because of its geometrical
interpretation. There is no any special ”ADM procedure”: Arnowitt, Deser and
Misner constructed the Hamiltonian dynamics following exactly the Dirac
scheme, just making use of another variables. The fact that two Hamiltonian
formulations (both according to the Dirac scheme, but the one for original
variables and the other for the ADM variables) are not related by canonical
transformations, should not lead to any bad-grounded conclusions like the one
made in [11], p. 68, that the gravitational Lagrangian used by Dirac and the
ADM Lagrangian are not equivalent. At the Lagrangian level, the transition to
the ADM variables is nothing more as a change of variables in the Einstein
equations, and there are no mathematical rules that would prohibit such change
of variables. It is the Lagrangian formulation of General Relativity which is
original and fundamental while its Hamiltonian formulation still remains
questionable, in spite of fifty years passed after Dirac’s paper [3]. The
extended phase space approach treating all degrees of freedom on an equal
footing may be a real alternative to the Dirac generalization of Hamiltonian
dynamics.

The example considered in Section 4 shows that the BRST charge can play the
role of a sought generator in extended phase space. Nevertheless, the
algorithm suggested by BFV for constructing the BRST charge again relies upon
the algebra of constraints. Even for the model from Section 4 it would not
lead to the correct result (4.30). Another way is to construct the BRST charge
as a conserved quantity based on BRST-invariance of the action and making use
of the first Noether theorem. This method works satisfactory for simple models
with a given symmetry. Below we mentioned that the gravitational Lagrangian
density can be modified for the primary constraints to take the simplest form
$\pi^{\mu}=0$ without affecting the equation of motion. However, after this
modification the full action may not be BRST-invariant. Some authors (see, for
example, [12, 13]) use some boundary conditions to exclude total derivatives
and ensure BRST-invariance. The boundary conditions (as a rule, these are
trivial boundary conditions for ghosts and $\pi^{\mu}$) correspond to
asymptotic states and are well-grounded in ordinary quantum field theory. This
way does not seem to be general enough, and for gravitational field the
justification of the boundary conditions, as well as the control of BRST-
invariance of the action, requires special study.

In [3] Dirac pointed out that “any dynamical theory must first be put in the
Hamiltonian form before one can quantize it”. Based upon Hamiltonian dynamics
in extended phase space, a new approach to quantum theory of gravity has been
proposed in [14, 15]. Ir was argued that it is impossible to construct a
mathematically consistent quantum theory of gravity without taking into
account the role of gauge degrees of freedom in description of quantum
gravitational phenomena from the point of view of different observers. The
present paper show that even at the classical level gauge degrees of freedom
cannot be excluded from consideration. As we have seen, the extension of phase
space by introducing the missing velocities changes the relations between the
“old” and “new” momenta (see (5.18)). As a consequence, the transformation
(2.9) is canonical. In that way, we consider extended phase space not just as
an auxiliary construction which enables one to compensate residual degrees of
freedom and regularize a path integral, as it was in the Batalin – Fradkin –
Vilkovisky approach [4, 5, 6], but rather as a structure that ensures
equivalence of Hamiltonian dynamics for a constrained system and Lagrangian
formulation of the original theory.

## Acknowledgements

I would like to thank Giovanni Montani and Francesco Cianfrani for attracting
my attention to the paper [11] and discussions.

## References

  * [1] P. A. M. Dirac, Can. J. Math. 2 (1950), P. 129–148.
  * [2] P. A. M. Dirac, Proc. Roy. Soc. A246 (1958), P. 326–332.
  * [3] P. A. M. Dirac, Proc. Roy. Soc. A246 (1958), P. 333–343.
  * [4] E. S. Fradkin and G. A. Vilkovisky, Phys. Lett B55 (1975), P. 224–226.
  * [5] I. A. Batalin and G. A. Vilkovisky, Phys. Lett B69 (1977), P. 309–312.
  * [6] E. S. Fradkin and T. E. Fradkina, Phys. Lett B72 (1978), P. 343–348.
  * [7] L. Castellani, Ann. Phys. 143 (1982), P. 357–371.
  * [8] R. Banerjee, H. J. Rothe and K. D. Rothe, Phys. Lett. B463 (1999), P. 248–251.
  * [9] R. Arnowitt, S. Deser and C. W. Misner, “The Dynamics of General Relativity”, in: Gravitation: an Introduction to Current Research, ed. by L. Witten, John Wiley & Sons, New York (1962), P. 227–265.
  * [10] L. D. Faddeev, Usp. Fiz. Nauk 136 (1982), P. 435–457 [Sov. Phys. Usp. 25 (1982). P. 130–142].
  * [11] N. Kiriushcheva and S. V. Kuzmin, “The Hamiltonian formulation of General Relativity: myth and reality”, E-print arXiv: gr-qc/0809.0097.
  * [12] J. J. Halliwell, Phys. Rev. D38 (1988), P. 2468–2481.
  * [13] M. Hennaux, Phys. Rep. 126 (1985), P. 1–66.
  * [14] V. A. Savchenko, T. P. Shestakova and G. M. Vereshkov, Gravitation & Cosmology 7 (2001), P. 18–28.
  * [15] V. A. Savchenko, T. P. Shestakova and G. M. Vereshkov, Gravitation & Cosmology 7 (2001), P. 102–116.


#  DiffSRL: Learning Dynamic-aware State Representation for Deformable Object
Control with Differentiable Simulation

Sirui Chen∗, Yunhao Liu∗, Jialong Li, Shang Wen Yao, Tingxiang Fan, Jia Pan† *
denotes equal contribution. † denotes the corresponding author.S. Chen, Y.
Liu, J. Li, S. Yao, T. Fan, J. Pan are with the University of Hong Kong.

###### Abstract

Dynamic state representation learning is an important task in robot learning.
Latent space that can capture dynamics related information has wide
application in areas such as accelerating model free reinforcement learning,
closing the simulation to reality gap, as well as reducing the motion planning
complexity. However, current dynamic state representation learning methods
scale poorly on complex dynamic systems such as deformable objects, and cannot
directly embed well defined simulation function into the training pipeline. We
propose DiffSRL, a dynamic state representation learning pipeline utilizing
differentiable simulation that can embed complex dynamics models as part of
the end-to-end training. We also integrate differentiable dynamic constraints
as part of the pipeline which provide incentives for the latent state to be
aware of dynamical constraints. We further establish a state representation
learning benchmark on a soft-body simulation system, PlasticineLab, and our
model demonstrates superior performance in terms of capturing long-term
dynamics as well as reward prediction. The source code and more experiments
results is available at https://ericcsr.github.io/DiffSRL/.

##  I Introduction

Deep neural networks have become a powerful tool in representing high
dimensional data using fewer dimensions. Its rapid development and application
nowadays have significantly enhanced and accelerated the machines to process
complex data such as images, words, sentences, and graphs. Among the
applications, one that has recently aroused great research interest is
representing real-world objects in dynamical environments from high
dimensional point clouds or images for control. Low dimension latent state can
benefit learning algorithms by reducing the trainable parameters, accelerating
the training, as well as achieving better control result. The efficiency of
this paradigm on rigid-body objects has been demonstrated in recent research
[1, 2].

Figure 1: Similar states but expecting distinct action trajectories.

However, representing more commonly seen deformable objects in low dimension
remains to be a challenging problem [3] and is essential for deformable object
manipulation. Representing deformable objects is difficult for multiple
reasons: a) The state of objects tends to be extremely large due to their high
degrees of freedom (DoFs) during deformation [3]. For example, a piece of
cloth can be translated integrally like an ordinary rigid object. But it can
also be rolled, crumpled, and folded with no predetermined pattern. Due to the
complexity resulting from high DoFs, soft-body objects lack uniform and
convenient representations but rely on surface particles or meshes to record
their current states. b) Different state elements of the same deformable
object can have very different correlations when undergoing different
deformation [3]. For instance, surface points on a piece of unfolded cloth,
when being translated forward, all obey the same movement. However, once being
folded, the cloth divides into upper and lower parts where surface points move
in opposite directions. c) Visually similar states can generate distinct
outcomes, which can confuse most static state representation learning methods.
Fig. 1 shows one such example where a robotic arm is removing a rod from a
piece of plasticine and then placing it at a target position. The initial
states of the given two situations are similar but with subtle differences: on
the left, the rod is completely surrounded by plasticine while on the right
there is a large enough notch in the plasticine. The difference between them
can be overlooked by static representation learning approaches because the
notch occupies merely a tiny number of volume. As a result, the controller
will generate resemblant actions while we prefer the robotic arm to first drag
the rod out (action  1) and then translate it (action  2) in the first case,
and a direct translation towards the target is more desirable for the second
one. Similar challenges are ubiquitous in soft-body object manipulation,
calling for the development of new state representation learning algorithms
that are dynamic-aware, i.e., the representation shall not only examine the
current state but also understand the dynamics to reflect different futures
bifurcated from similar current states.

Existing state representation learning methods, such as AutoEncoder [2] and
Embed to Control (E2C) [1], lack specific designs for dealing with the above
challenges brought by the dynamics of deformable objects. Newly emerged
solutions focusing on deformable objects, such as CFM [3] or G-Doom [4], as
will be elaborated in Sec. II, cannot directly embed simulation into their
training pipelines, leading to the performance bottleneck on soft-body object
manipulation tasks [5]. To improve the performance of state representation
learning on deformable objects, we hereby propose a new pipeline DiffSRL that
utilizes a differentiable simulator to better encode dynamic- and constraint-
aware information.

Our main contributions are threefold.

  * •

We propose an end-to-end dynamic state representation learning pipeline
utilizing differentiable simulation to capture multistep dynamic of deformable
objects. To the best of our knowledge, this is the first time that state
representation learning directly uses the ground truth simulation function as
part of the end-to-end training pipeline.

  * •

Our method explicitly considers dynamics constraints related to the training
pipeline and demonstrates awareness of physical constraints.

  * •

We establish a benchmark on PlasticineLab [6], a non-Newtonian deformable
object simulation environment. As baselines, we also implement a set of state-
of-the-art state representation learning algorithms, such as CFM [3], inverse
model learning [7], AutoEncoder [2] and E2C [1].

We have compared our proposed pipeline and baselines according to the reward
prediction accuracy and the trajectory reconstruction errors, and investigated
the performance of the learned encoder on important downstream tasks such as
model-free reinforcement learning and model-based policy optimization. We have
demonstrated that our approach outperforms all baselines on all tasks most of
the time.

##  II Related Work

Figure 2: Overview of our DiffSRL model pipeline.

Dynamic state representation learning is an emerging topic together with
reinforcement learning for the acceleration and enhancement it offers. Auto-
Encoder [2] trains a CNN encoder-decoder pair to extract features from images
to achieve better performances in deep reinforcement learning than training
all parameters from scratch. This method, however, only uses static state
information and will fail to distinguish visually similar but dynamically
different states [5] as previously shown. As a solution to this issue, Embed
to Control (E2C) [1] proposes a pipeline including a learned forward model, an
auto-encoder, and a linear dynamic transfer constraint between each current
state and the next state. E2C demonstrates its strength in enhancing deep RL’s
performance in target-reaching problems with image observations, but its
linear model cannot well describe the complex dynamics of deformable objects
especially for multiple steps. Replacing the linear model in E2C with a neural
network is not trivial because the neural net with significantly more
parameters may be hard to converge and prevent the E2C variant from achieving
desirable outcomes as well as failed to rollout accurate long run
trajectories. Alternatively, the Inverse-model based state representation
learning [8] is based on the heuristics that good state representation should
be able to predict action leading to the transfer between two consecutive
states. Thus, the encoder can be trained with an action predictor by
minimizing the MSE loss between real and predicted action. Inverse model based
method has been further applied in model-free RL tasks [9, 10], with more
variants and applications elucidated in [5], but existing designs do not take
the complexity of deformable objects into consideration either. CFM [3]
explored the area of state representation learning on deformable objects such
as clothes and ropes. It utilizes contrastive learning to jointly train both
the forward model and the encoder to represent images of deformable objects
with distinguishable latent states. The trained forward model and latent space
has been used in some MPC-based manipulation control. G-Doom [4] applies a
similar paradigm as CFM but with a different architecture in terms of using
latent graphs as the representation in the latent space and graph dynamic
models as the forward model. However, despite their achievement in
manipulating deformable objects, both CFM and G-Doom assume no access to the
ground-truth analytical model of the dynamic system and didn’t take physical
constraints into consideration, which may lead to their forward models’
incapability in accurate long-term planning.

Recent advances of physically-based simulation, such as Material Point Method
(MPM) [11], allow us to accurate and efficiently simulate deformable objects
such as cloth and plasticine. What’s more, the differentiable simulator [12,
13, 14] allows each simulation step to be integrated as a layer in deep neural
networks, making end-to-end training possible.This provide opportunity to
overcome the difficulty of capture complex dynamic with a neural net faced by
E2C, CFM and G-DOOM. In end-to-end training, gradients from the upstream
neural networks are propagated through physics engines and affect network
parameters accordingly. For instance, DiffTaichi [14], a newly-emerged
differentiable physics engine based on Taichi programming language [15], has
included many MPM simulation examples and established a soft-body control
benchmark with differentiable simulation in PlasticineLab [6]. We will also
perform our experiments on PlasticineLab for a fair comparison among different
methods.

As for the states to be encoded, common methods take 2D images as the input
and thus adopt the convolution methods. But they may encounter difficulties
when being applied in point clouds, due to the intrinsic permutation
irregularity of point clouds. PointNet [16] is a pioneering work of applying
deep learning on point cloud data, where the permutation invariance is
achieved using point-wise symmetric functions, including point-wise MLP and
max-pooling layers. Moreover, being trained on point clouds allows the deep
learning algorithms to tolerate the variety of input shapes, regardless of how
many points are contained in the incoming point clouds. Thus, our proposed
DiffSRL method will utilize point clouds to extract samples from the
PlasticineLab environment as the training dataset.

##  III Problem Formulation

A general dynamic system consists of three components: a state space
$\mathcal{S}$, an action space $\mathcal{A}$, and a transfer function:

$\displaystyle s_{t+1}=f_{\text{sim}}(s_{t},a_{t}),\
s_{t},s_{t+1}\in\mathcal{S};a_{t}\in\mathcal{A}.$ (1)

The state $s_{t}$, obtained from the differentiable simulator, is a mixture of
observable parts (such as the particles’ coordinates) and imperceptible parts
(including the velocity fields and constraints). The observable part of the
state $s_{t}$ is denoted as $s^{\text{obs}}_{t}$. The sensors of the robotic
system may only observe part of the entire observable surroundings (denoted as
$\mathcal{O}$) due to its limited visual angles or occlusion of obstacles.
This can be modeled as an observation function
$o_{t}=g_{\text{obs}}(s^{\text{obs}}_{t})$, which selects accessible
observation from observable states. Our goal is to establish a dynamic system
state representation learning method capable of computing a low dimensional
latent state $h_{t}$ mapped from observation:
$h_{t}=f^{\theta}(o_{t}),o\in\mathcal{O}$. The mapping function $f^{\theta}$
is expected to capture sufficient information from the current state
observation related to dynamics so that it can figure out whether two latent
states are dynamically equivalent or not. In particular, two latent states
$h_{t}$ and $h_{t}^{\prime}$ are regarded as dynamically equivalent if given
an arbitrary action sequence $a_{t},a_{t+1}...,a_{t+k}$, their future states
$s_{t+k},s^{\prime}_{t+k}$ will be similar.

##  IV Approach

To learn the state representation that is aware of dynamic-equivalence, we
propose the DiffSRL framework, which consists of three main components: 1) a
state observation Autoencoder, which compresses the dynamic state to a latent
state and reconstructs the state from latent state; 2) a constraint regulator,
which regulates the decoded state against system dynamical constraints, and 3)
a differentiable simulator, which rolls out the trajectory starting from the
decoded state. The rolling-out process is differentiable so that gradients can
be propagated from the end of the trajectory all the way back to the encoder.
The overall pipeline of DiffSRL is shown in Fig. 2.

In this section, we will first introduce some preliminary knowledge of the
Particle-In-Cell paradigm and the differentiability relevant to the simulator
component in our design, and then move on to design details of the three major
components of our model as well as our loss function.

###  IV-A Preliminaries

Particle-In-Cell paradigm: Lagrangian methods and Euler methods are two major
branches of simulation algorithms [11]. The former method carries physical
property with particles and is convenient for numerical integration; while the
latter one stores physical properties on fixed grids and can be used for fast
interactive simulation. The Particle-In-Cell (PIC) method is a combination of
both methods and is commonly applied in fluid and deformable object
simulation. The deformable object simulation in PlasticineLab is based on
Material Point Method [11], a variant of Particle-In-Cell. Our constraint
regulator also utilizes PIC paradigm to efficiently detect penetration.

Differentiable Simulation: Physically-based simulation in Eq. 1 maps the
current state $s_{t}$ and action $a_{t}$ to the next state $s_{t+1}$.
Traditional simulators such as Bullet [17], Dart [18], Mujoco [19] only
provide non-differentiable forward functions. Thanks to the recent advances in
differentiable simulation Sec. II, now the gradient of simulation function can
be also acquired:

$\displaystyle\nabla{s}_{t}$
$\displaystyle=\nabla{s}_{t+1}\text{Jac}^{s_{t+1}}_{s_{t}}f(s_{t},a_{t})$ (2)
$\displaystyle\nabla{a}_{t}$
$\displaystyle=\nabla{s}_{t+1}\text{Jac}^{s_{t+1}}_{a_{t}}f(s_{t},a_{t})$

which allows us to embed the simulation function into a deep learning pipeline
such as Pytorch and Tensorflow as part of the end-to-end training pipeline. In
DiffSRL’s architecture, the trajectory rolling-out module is based on
differentiable simulation.

Distance Metric on Point Clouds: Point clouds usually are expressed as an
unordered set whose distance cannot be directly measured by common metrics
such as MSE or MAE. Instead, here we use the Earth Mover’s Distance (EMD) [20]
to measure the difference between two point clouds in terms of the minimum
total pairwise distance between all the points. The EMD distance can be
expressed as:

$\displaystyle\text{Dist}_{\text{emd}}(A,B)=\min_{\phi:A\rightarrow
B}\sum_{a\in A}d(a,\phi(a))\vspace{-5pt}$ (3)

where $A$ and $B$ are two unordered sets, $d(\cdot,\cdot)$ is the distance
between two points from these two sets, and $\phi(\cdot)$ is the
correspondence between $a$ and $b$ to be optimized. The minimization is
implemented using Iterative Closest Point (ICP) [21] algorithm. Due to the
high computational cost when solving this optimization, in practice, a greedy
point matching based variant of EMD, Chamfer Distance [22], is also commonly
used, which averages all the nearest neighbor distance by nearest neighbors:

$\displaystyle\text{Dist}_{\text{chamfer}}(A,B)=\sum_{a\in A}\min_{b\in
B}d(a,b)+\sum_{b\in B}\min_{x\in A}d(a,b).$ (4)

We will use the losses based on Chamfer Distance and ICP algorithm as the
constraint regulator.

###  IV-B State Observation Autoencoder

To obtain the latent states of a dynamic system, we use a deep neural network
with learnable parameters $\theta_{\text{encoder}}$ as the encoder, where the
specific architecture depends on the type of system’s observation. In
particular, we will use Multi-layer Perceptron (MLP), Convolution Neural
Network (CNN), and permutation invariant encoders such as PointNet [16] for
observations in form of an ordered vector, image, or an unordered set of
particles such as point clouds, respectively. Meanwhile, a decoder with
parameters $\theta_{\text{decoder}}$ is trained simultaneously for recovering
the observable state from the latent state. Formally, the Autoencoder can be
expressed as:

$\displaystyle h_{t}$ $\displaystyle=f^{\theta_{\text{encoder}}}(o_{t}),$ (5)
$\displaystyle\tilde{s}^{\text{obs}}_{t}$
$\displaystyle=f^{\theta_{\text{decoder}}}(h_{t}),$

where $o_{t},s^{obs}_{t},h_{t}$ are as defined in Sec. III and
$\tilde{s}^{\text{obs}}_{t}$ is the reconstructed observation from $h_{t}$.

###  IV-C Differentiable Constraint Regulator

After the reconstructed observable state $\tilde{s}^{\text{obs}}_{t}$ is
obtained from the decoder, the reconstructed state $\tilde{s}_{t}$ can be
obtained by replacing observable part $s^{\text{obs}}_{t}$ from $s_{t}$ with
$\tilde{s}^{\text{obs}}_{t}$. However, the reconstructed state may violate
hard dynamics constraints, such as non-interpenetration constraints, joint
limit constraints, as well as continuity constraints within continuous
deformable objects. Constraint violation usually causes simulator failure or
significant numerical errors. Hence, it is necessary to maintain the encoder’s
awareness of the constraints from observation throughout the training
procedure to enforce the gradient to be accurate and meaningful. Therefore,
the constraints regulator is designed to find the feasible state closest to
the generated state and, meanwhile, punish the Autoencoder for reconstructing
unrealistic states. Formally, the reconstructed state, the constraint
regulator loss, and the regulated state are computed as

$\displaystyle\tilde{s}_{t}$ $\displaystyle=s_{t}\ominus
s^{\text{obs}}_{t}\oplus\tilde{s}^{\text{obs}}_{t},$ (6)
$\displaystyle\mathcal{L}_{\text{constraint}}$
$\displaystyle=\min_{s\in\mathcal{S}}\text{Dist}(s,\tilde{s}_{t}),$ (7)
$\displaystyle\hat{s}_{t}$
$\displaystyle=\arg\min_{s\in\mathcal{S}}\text{Dist}(s,\tilde{s}_{t}),$ (8)

where $\ominus$ and $\oplus$ represent removing and adding back the observable
parts from and to the full state, respectively. For ordered state vectors, the
distance Dist can be measured directly as the weighted MSE between the two
vectors; as for unordered states such as point clouds, we use EMD as the
distance metric.

###  IV-D Differentiable Simulator and loss design

By using the reconstructed state, the simulation in Eq. 1 then becomes

$\displaystyle\hat{s}_{t+1}$ $\displaystyle=f_{\text{sim}}(\hat{s}_{t},a_{t})$
(9)

where $\hat{s}_{t}$ is reconstructed state after regulation, onto which
$f_{\text{sim}}$ executes the input action $a_{t}$ in the simulation world and
reaches the corresponding result state $\hat{s}_{t+1}$. When being applied to
an trajectory of length $k$, we obtain a sequence:

$\displaystyle\hat{s}_{t:t+k}$
$\displaystyle=\text{DiffRollout}(\hat{s}_{t},a_{t:t+k})$ (10)

where DiffRollout refers to the successive execution of the simulator,
starting from $\hat{s}_{t}$, along the trajectory consisting of an action
sequence $a_{t},a_{t+1},...,a_{t+k}$.

The model is trained with both the constraints violation loss
$\mathcal{L}_{\text{constraint}}$ (Eq. 7) and the multi-step reconstruction
loss $\mathcal{L}_{\text{multi-step}}$ defined as

$\displaystyle\mathcal{L}_{\text{multi-
step}}=\sum_{i=1}^{k}\gamma\text{Dist}(s_{t+i},\hat{s}_{t+i}),$ (11)

which punishes the distance between each state $s_{t+i}$ in the ground truth
trajectory and its correspondence in the rolling-out trajectory starting from
the reconstructed state. A decaying factor $\gamma$ is used to mitigate the
gradient instability effect when performing back propagation through a long
horizon. In this way, the total loss to be optimized is

$\displaystyle\mathcal{L}_{\text{total}}=(1-\beta)\mathcal{L}_{\text{multi-
step}}+\beta\mathcal{L}_{\text{constraints}}$ (12)

where the weight factor $\beta$ is introduced to trade off among different
terms through the entire training scheme. Initially, the state reconstructed
by the decoder violates constraints significantly, so $\beta$ is set to $1$ to
encourage the Autoencoder to respect the constraints in the physics system.
Progressively, $\beta$ decays exponentially at rate $\lambda$ so that the
multi-step reconstruction loss dominates.

##  V Experiments and Results

To evaluate effectiveness of DiffSRL, we evaluated our model with three
experiments: 1) trajectory reconstruction, 2) reward prediction, and 3) model-
free reinforcement learning on Chopsticks and Rope tasks from the
PlasticineLab benchmark [6]. Our proposed DiffSRL has achieved state-of-the-
art performance in these experiments most of the time.

###  V-A Baselines for Comparison

We compared our approach against state-of-the-art method CFM [3] and four
commonly used dynamic state representation learning methods: E2C [1], Forward
[23] that learns a dynamic forward model on latent states, Inverse [7] that
conducts inverse model learning for action prediction, as well as AutoEncoder
directly reconstructing the state [2]. In Tab. I, we summarize characteristics
of different approaches.

In terms of details of baseline implementation, similar to our method, we use
Chamfer Distance [22] as loss function to train autoencoder in the baseline.
The forward model used in E2C and Forward
$f^{\theta}_{\text{fwd}}(h_{t},a_{t})$ is 3 layer MLP on latent space, the
action predictor
$a^{\text{pred}}_{t}=f^{\phi}_{\text{inverse}}(h_{t},h_{t+1})$ in inverse
model [7] uses concatenation of two latent states as input and predicts the
action with a 3 layer MLP. The contrastive loss in [3] uses the same formation
and negative sample selection method as in the original paper. All MLPs
mentioned here use 256 as the size of a hidden layer.

|  predictive dynamic model |  predict action |  observation reconstruction |  multistep information |  physical constraints
---|---|---|---|---|---
DiffSRL (Our) |  $\surd$ |  |  $\surd$ |  $\surd$ |  $\surd$
E2C[1] |  $\surd$ |  |  $\surd$ |  |
Forward[23] |  $\surd$ |  |  |  |
Inverse[7] |  |  $\surd$ |  |  |
CFM[3] |  $\surd$ |  |  |  |
Autoencoder[2] |  |  |  $\surd$ |  |

TABLE I: Comparison of different components used in different state
representation learning pipelines. Figure 3: Collision-free and smooth
velocity field constraint regulator, where the dots and the capsule represent
the point clouds and the rigid-body object which may collide with it,
respectively.

###  V-B Implementation Details

Since PlasticineLab is a deformable object simulator based on Material Point
Method (MPM) [11], the particles used for simulation can be naturally
considered as an unordered point cloud. For simplicity, we choose each
observation $o_{t}$ as the position of all particles at time $t$. An
observation, together with velocity and other particle wise properties, forms
a full state $s_{t}$. A variant of PCN [24] is then used to implement the
Autoencoder for the point cloud $o_{t}$, with the detailed architecture
presented in Fig. 2.

As for the constraint regulator, the MPM simulation disallows penetration
between particles and rigid bodies, and therefore requires smooth velocity
fields within continuous materials. Hence, our regulator is constructed in two
stages: 1) The signed distance field (SDF) algorithm with MPM is used to
enforce non-penetration and detect collisions efficiently. For those particles
$p_{\text{collide}}$ resided inside rigid bodies, we resolve the collision by
finding the minimum distance $d_{p_{\text{collide}}}$ to exit the colliding
object. For simplicity, it is assumed that each rigid-body object is a
geometric primitive and has analytical collision detection routine, so that
the minimum exit distance can be solved efficiently using a simple linear
programming solver. The loss function for this objective is defined as the sum
of all $d_{p_{\text{collide}}}$. 2) The Iterative Cloest Point (ICP) algorithm
is applied to find the best pairwise matches between particles in ground truth
observation and those in the reconstructed one where some particles have been
translated outside rigid-body objects to avoid collisions. This step is
necessary because the decoded observation $\tilde{o}_{t}$ are expected to
satisfy soft-body dynamic properties from $s_{t}$ such as continuous velocity
fields. After pairwise relationships are established by ICP, the velocity and
smoothness constraints of each ground-true particle are exerted to the matched
pairs. The details are described in Algorithm 1 and also illustrated in Fig.
3.

Algorithm 1 Constraints Regulator

Input: Rigid bodies $\mathcal{R}$, the ground-truth particle observable state
$s^{\text{obs}}_{t}$ (i.e., position of particles), and decoded particles
position $\tilde{s}^{\text{obs}}_{t}$

Initialize: Grids to store interpenetrating rigid bodies information
$c_{\text{grid}}^{\text{col}}$, loss of penetration $l_{\text{penetration}}$,
and loss reconstruction $l_{\text{rec}}$.

$\hat{s}_{t}^{\text{obs, no-
penetration}}\leftarrow\text{ResolvePenetration}(\mathcal{R},\tilde{s}^{\text{obs}}_{t})$

$\triangleleft$ Compute grid mass using particle-in-cell

$m_{\text{grid}}\leftarrow\text{ComputeGridMass}(\tilde{s}^{\text{obs}}_{t})$

for $r$ in $\mathcal{R}$ do

$\triangleleft$ Check penetration using signed distance field

$c_{\text{grid}}^{\text{col}}\leftarrow\text{SDF}(r,m_{\text{grid}})$

end for

$\triangleleft$ Compute Interpenetration related information of each particle
to all rigid bodies

for $p$ in $\hat{s}_{t}^{obs\text{no-penetration}}$ do

if p.grid is not empty then

$\triangleleft$ Solve linear programming problem to find minimum displacement
for a particle to stay out of penetrating rigid body

$\text{p.pos}\leftarrow\text{LPSolver}(\mathcal{R},p)$

$\triangleleft$ Sum up minimum displacement costs as loss

$l_{\text{penetration}}\mathrel{+}=\text{LPCost}(r,p)$

end if

end for

$\triangleleft$
$\hat{s}^{\text{obs}}_{t}\leftarrow\text{ResolveSmoothness}(\hat{s}_{t}^{\text{obs
no-penetration}})$

$\text{order}\leftarrow\text{ICPSolver}(s^{\text{obs}}_{t},\hat{s}^{\text{obs}}_{t})$

$\hat{s}^{\text{obs}}_{t}\leftarrow\text{perm}(\hat{s}_{t}^{\text{obs no-
penetration}},\text{perm(order)})$

return $\hat{s}^{\text{obs}}_{t}$

The distance between the regulated observation and the corresponding ground
truth is the sum of pairwise squared distance, i.e., the EMD [20] distance
between $\tilde{o}_{t}$ and $o_{t}$.

###  V-C Data Collection

We collect 6,000 trajectories of length 8 from benchmarks Chopsticks and Rope.
To make the datasets cover as many states as possible, 30% of data are sampled
from random policies while others are collected when optimizing trajectories
for manipulating different targets. This is because random actions frequently
cause the manipulator to detach from the deformable object too early to
capture sufficient deformation information. We use the ratio of 10:1:1 to
split collected samples into training, validation and testing sets.

Datasets | Rope | Chopsticks
---|---|---
Metric | ChamferDist | MSE | ChamferDist | MSE
AutoEncoder | 0.048 | 0.236 | 0.044 | 1.102
Inverse | 0.042 | 0.276 | 0.048 | 1.19
Forward | 0.083 | 5.435 | 0.124 | 6.599
E2C | 0.048 | 0.223 | 0.064 | 1.244
CFM | 0.037 | 0.236 | 0.037 | 1.010
DiffSRL | 0.027 | 0.218 | 0.031 | 0.840
DiffSRL w/o Reg | 0.039 | 0.292 | 0.036 | 1.226

TABLE II: Comparison of the average and standard deviation of trajectory
rollout chamfer distance and mean squared error of reward prediction. The
smaller the MSE and Chamfer Distance, the better the model. The best result is
highlighted in bold and the second best is highlighted with underline. DiffSRL
w/o Reg is our proposed method trained without non-penetration constraint
regulator.

###  V-D Comparison in Trajectory Reconstruction

A straight forward evaluation of dynamical equivalence is the similarity
between the real trajectory and the trajectory roll-out from the reconstructed
state given the same action sequence, which is computed by accumulating the
Chamfer Distance [22], i.e.,
$\sum_{s_{t:t+k}}D_{\text{chamfer}}(s_{t},\hat{s}_{t})$, between each pairs of
states in both trajectories. Some baseline methods such as CFM, Forward and
Inverse only train the encoder without decoder and thus cannot reconstruct
states. Thus, we trained a decoder for each fixed encoder using the
Autoencoder training pipeline until convergence. Tab. II shows that DiffSRL
achieves the best performance among all models. Moreover, to evaluate the
effectiveness of the constraint regulator, we conduct an ablation study by
removing non-penetration regulation modules during the training phase. It
shows the removal of non-penetration regulator will result in higher
trajectory reconstruction error as well as reward prediction error. We also
tried to remove the smoothness regulator during training, but this will
results in NaN in gradient from differentiable simulator since the MPM’s
requirement that particles in the same grid have consistent properties can be
violated. One thing worth noting is that the Forward model is the worst among
all models, which may be because it tends to be trapped in a trivial minimum
as mentioned in [5].

###  V-E Comparison in Reward Prediction

Figure 4: The architecture of policy network that we used in model-based
policy optimization and model-free reinforcement learning.

For many downstream control tasks, it is essential that the encoded latent
states can capture the state’s main characters for accurate reward prediction.
To predict rewards from latent states, we train a 3-layer MLP for latent
states obtained from each encoder. The reward that is related to poses of the
deformable object is the only part we predicted. The error is measured in mean
square error:
$\frac{1}{|\mathcal{D}|}\sum_{s\in|\mathcal{D}|}\|\hat{r}_{\text{pred}}-r_{\text{real}}\|^{2}$.
As summarized in Tab. II, our model achieves the best performance on both
datasets. The reward prediction error also increases when non-penetration
regulator is absent.

###  V-F Model-Based Policy Optimization (MBPO)

Figure 5: TD3 Reward in different environments. Figure 6: Keyframes of
different approaches in Model-Based Policy Optimization (MBPO) task.

To investigate adaptability as well as effectiveness on various more complex
downstream tasks, we first evaluated the encoders trained from different state
representation learning methods on end-to-end policy optimization with the
differentiable simulator. Our experiment setting is similar to the end-to-end
policy optimization in DiffTaichi [14, 25]. The trained policy is expected to
control manipulators to move the plasticine to different target positions
within finite steps. Detailed target and reward design are the same as
described in [6]. The architecture of the policy network is shown in Fig. 4,
where the encoder’s weights are fixed during policy optimization. For
comparison, we add another baseline of the default down-sampling based method
provided by PlasticineLab, which samples 400 points from each point cloud and
keeps track of those points in the sampling order throughout the
trajectory111The down-sampling based method is not feasible in practice since
tracking 400 points all the way on deformable object accurately is difficult..
The policy network is initialized as a simple MLP and its inputs are obtained
by concatenating features. We evaluated performance using the per-epoch
average reward. Each experiment is repeated 5 times and the means and standard
deviations are plotted in Fig. 5. The first row illustrates the reward of
different methods in Chopstick and Rope environments, where most methods using
state representation learning can improve their policies faster and achieve
better final rewards, compared with the down-sampling based method. Moreover,
our approach achieves the best result among all state representation learning
methods including the AutoEncoder that is trained using a similar loss
function. Surprisingly, the Forward method achieves competitive performance
compared with other methods even though it has the worst performance in
previous two tasks, which may indicate that the performance correlation
between MBPO and reward / trajectory reconstruction is not significant and we
will leave the study to the future work. Some key frames from the trajectory
are shown in Fig. 6.

Figure 7: Model based policy optimization reward in rope environment when
using different number of point as encoder input.

###  V-G Model-Free Reinforcement Learning (MFRL)

Finally, we tested all models on model-free reinforcement learning on the same
target-reaching tasks as described in MBPO. Since the trajectory information
has been used in MBPO experiments similar to on-policy reinforcement learning,
to fully evaluate different state representation learning methods, our
experiment uses the state-of-the-art off-policy algorithm TD3 [26], which uses
latent state for both state value function and policy function. The
architecture of policy function is shown in Fig. 4. We use similar experiments
and evaluation criteria as MBPO but increase the epoch number from 200 to
1,000 since the model-free method requires more data. The results are shown in
Fig. 5, from which we can observe that most state representation learning
methods have higher maximum rewards and faster convergence speeds than the
sample-based method, and DiffSRL achieves the best performance while CFM’s
performance in two environments are very different, implying its limitation in
general situations.

###  V-H Ablation Study

We conduct ablation study to show the effectiveness of constraints regulator.
Since the removal of smoothness constraint regulator will crash the
representation learning training procedure, we retrained the encoders using
our pipeline without non-penetration constraints regulator and test them using
MFRL and MBPO training to verify the effectiveness. Fig. 8 shows the result,
we can justify the effectiveness of our design since the reward drop
significantly if we remove the non-penetration regulator.

Figure 8: Ablation Study: without non-penetration regulator

###  V-I Robustness Analysis

The number of points observed by a sensor may vary greatly in real
applications. Since the point cloud encoder can be directly applied on
different point number without re-training, we use different numbers of
observed points to test the robustness of our model. We trained a DiffSRL
model using point clouds consisting of 8,192 points, and then investigated the
deployment performance using 6,000, 4,000, 2,000 particles, by computing the
reward curves when training MBPO from latent states. As shown in Fig. 7,
although using fewer observable particles (from 8,192 to 6,000) than training
will decrease the performance, further decreasing of the particle number (from
6,000 to 4,000 or 2,000) will not significantly downgrade the performance,
which demonstrates that our method is relatively robust against disturbance in
observation.

###  V-J Hyperparameters Setting

Figure 9: 10-step Chamfer loss on both rope and chopstick benchmarks when
using different number of training steps.

Parameters | Meaning | Value
---|---|---
N | Number of particles | 8192
I | ICP iterations | 3000
$d_{\text{latent}}$ | Latent Space dimension | 1024
$\alpha$ | Learning rate | 1e-5
$\gamma$ |  weight decay rate in 11 | 0.99
$\beta$ |  Initial weight between in 12 | 0.99
$\lambda$ |  Decay rate of $\beta$ per epoch | 0.9
E | Number of epochs | 20

TABLE III: Hyperparameters

Tab. III summarizes the hyperparameters used in our model. One important
hyperparameter of our method is $k$, the number of trajectory rollout steps
used when training the state representation. We choose $k$ by evaluating the
10-step trajectory reconstruction loss on the validation set, with the result
shown in Fig. 9. For both environments, 7-step achieves the best performance
and thus we choose $k=7$. Notice that the loss curve demonstrates non
monotonic behavior as we increase the training step size, and the possible
reason is that as the step size increases, the gradient back propagation
though time tends to be more noisy due to numerical issues and the
optimization landscape becomes more wiggly. We leave further investigation as
future work.

##  VI Conclusion and Future Work

We have presented a novel dynamic state representation learning model and its
sample efficient training scheme. Moreover, we evaluated our model using
multiple tasks on soft body benchmarks in PlasticineLab [6] and set up
benchmarks to compare with other state-of-the-art models. Currently our model
uses part of states directly from the simulator as observation, while in the
real world, robots may not necessarily have sensors to access such
information. In future, it is essential to extract the observation using
physical sensors such as Lidar or cameras for deploying our method on real
robots. One possible follow-up might include the differentiable rendering
method [27] in our pipeline, which enables using easily accessible images as
observation while maintaining the overall differentiability.

## References

  * [1] M. Watter, J. T. Springenberg, J. Boedecker, and M. A. Riedmiller, “Embed to control: A locally linear latent dynamics model for control from raw images,” in _NeurIPS_ , 2015, pp. 2746–2754.
  * [2] C. Finn, X. Y. Tan, Y. Duan, T. Darrell, S. Levine, and P. Abbeel, “Deep spatial autoencoders for visuomotor learning,” in _ICRA_ , 2016, pp. 512–519.
  * [3] W. Yan, A. Vangipuram, P. Abbeel, and L. Pinto, “Learning predictive representations for deformable objects using contrastive estimation,” _CoRR_ , vol. abs/2003.05436, 2020.
  * [4] X. Ma, D. Hsu, and W. S. Lee, “Learning latent graph dynamics for deformable object manipulation,” _CoRR_ , vol. abs/2104.12149, 2021.
  * [5] T. Lesort, N. D. Rodríguez, J. Goudou, and D. Filliat, “State representation learning for control: An overview,” _Neural Networks_ , vol. 108, pp. 379–392, 2018.
  * [6] Z. Huang, Y. Hu, T. Du, S. Zhou, H. Su, J. B. Tenenbaum, and C. Gan, “Plasticinelab: A soft-body manipulation benchmark with differentiable physics,” in _ICLR_ , 2021.
  * [7] P. Agrawal, A. Nair, P. Abbeel, J. Malik, and S. Levine, “Learning to poke by poking: Experiential learning of intuitive physics,” _CoRR_ , vol. abs/1606.07419, 2016.
  * [8] W. Duan and K. Jens, “Learning state representations for robotic control: Information disentangling and multi-modal learning,” Master’s thesis, TU Delft, 2017.
  * [9] E. Shelhamer, P. Mahmoudieh, M. Argus, and T. Darrell, “Loss is its own reward: Self-supervision for reinforcement learning,” in _ICLR_ , 2017\.
  * [10] D. Pathak, P. Agrawal, A. A. Efros, and T. Darrell, “Curiosity-driven exploration by self-supervised prediction,” in _CVPR Workshops_ , 2017, pp. 488–489.
  * [11] C. Jiang, C. A. Schroeder, J. Teran, A. Stomakhin, and A. Selle, “The material point method for simulating continuum materials,” in _SIGGRAPH_ , 2016, pp. 24:1–24:52.
  * [12] K. Werling, D. Omens, J. Lee, I. Exarchos, and C. K. Liu, “Fast and feature-complete differentiable physics engine for articulated rigid bodies with contact constraints,” in _RSS_ , 2021.
  * [13] J. Xu, T. Chen, L. Zlokapa, M. Foshey, W. Matusik, S. Sueda, and P. Agrawal, “An end-to-end differentiable framework for contact-aware robot design,” in _RSS_ , 2021.
  * [14] Y. Hu, L. Anderson, T. Li, Q. Sun, N. Carr, J. Ragan-Kelley, and F. Durand, “Difftaichi: Differentiable programming for physical simulation,” in _ICLR_ , 2020.
  * [15] Y. Hu, T. Li, L. Anderson, J. Ragan-Kelley, and F. Durand, “Taichi: a language for high-performance computation on spatially sparse data structures,” _ACM Trans. Graph._ , vol. 38, no. 6, pp. 201:1–201:16, 2019\.
  * [16] C. R. Qi, H. Su, K. Mo, and L. J. Guibas, “Pointnet: Deep learning on point sets for 3d classification and segmentation,” in _CVPR_ , 2017, pp. 77–85.
  * [17] E. Coumans _et al._ , “Bullet physics library,” _Open source: bulletphysics. org_ , vol. 15, no. 49, p. 5, 2013.
  * [18] J. Lee, M. X. Grey, S. Ha, T. Kunz, S. Jain, Y. Ye, S. S. Srinivasa, M. Stilman, and C. K. Liu, “DART: Dynamic animation and robotics toolkit,” _The Journal of Open Source Software_ , vol. 3, no. 22, p. 500, 2018.
  * [19] E. Todorov, T. Erez, and Y. Tassa, “Mujoco: A physics engine for model-based control,” in _IROS_ , 2012, pp. 5026–5033.
  * [20] Y. Rubner, C. Tomasi, and L. J. Guibas, “A metric for distributions with applications to image databases,” in _ICCV_ , 1998, pp. 59–66.
  * [21] K. S. Arun, T. S. Huang, and S. D. Blostein, “Least-squares fitting of two 3-d point sets,” _IEEE Trans. Pattern Anal. Mach. Intell._ , vol. 9, no. 5, pp. 698–700, 1987.
  * [22] M. Liu, O. Tuzel, A. Veeraraghavan, and R. Chellappa, “Fast directional chamfer matching,” in _CVPR_ , 2010, pp. 1696–1703.
  * [23] L. Kaiser, M. Babaeizadeh, P. Milos, B. Osinski, R. H. Campbell, K. Czechowski, D. Erhan, C. Finn, P. Kozakowski, S. Levine, A. Mohiuddin, R. Sepassi, G. Tucker, and H. Michalewski, “Model based reinforcement learning for atari,” in _ICLR_ , 2020.
  * [24] W. Yuan, T. Khot, D. Held, C. Mertz, and M. Hebert, “PCN: point completion network,” in _3DV_ , 2018, pp. 728–737.
  * [25] Y. Qiao, J. Liang, V. Koltun, and M. C. Lin, “Scalable differentiable physics for learning and control,” in _ICML_ , vol. 119, 2020, pp. 7847–7856.
  * [26] S. Fujimoto, H. van Hoof, and D. Meger, “Addressing function approximation error in actor-critic methods,” in _ICML_ , vol. 80, 2018, pp. 1582–1591.
  * [27] H. Kato, D. Beker, M. Morariu, T. Ando, T. Matsuoka, W. Kehl, and A. Gaidon, “Differentiable rendering: A survey,” _CoRR_ , vol. abs/2006.12057, 2020\.