Hamiltonian Monte Carlo On Lie Groups and Constrained Mechanics on Homogeneous Manifolds

TL;DR

This paper unifies Hamiltonian Monte Carlo methods on Lie groups with geometric mechanics, extending algorithms to non-compact groups and formulating mechanics on homogeneous spaces as constrained systems.

Contribution

It establishes a geometric framework linking HMC on Lie groups with classical mechanics, enabling extensions to non-compact groups and constrained systems on homogeneous spaces.

Findings

Recovered HMC on compact Lie groups from geometric mechanics.

Extended HMC algorithms to non-compact Lie groups.

Formulated mechanics on homogeneous spaces as constrained systems.

Abstract

In this paper we show that the Hamiltonian Monte Carlo method for compact Lie groups constructed in \cite{kennedy88b} using a symplectic structure can be recovered from canonical geometric mechanics with a bi-invariant metric. Hence we obtain the correspondence between the various formulations of Hamiltonian mechanics on Lie groups, and their induced HMC algorithms. Working on $\G\times \g$ we recover the Euler-Arnold formulation of geodesic motion, and construct explicit HMC schemes that extend \cite{kennedy88b,Kennedy:2012} to non-compact Lie groups by choosing metrics with appropriate invariances. Finally we explain how mechanics on homogeneous spaces can be formulated as a constrained system over their associated Lie groups, and how in some important cases the constraints can be naturally handled by the symmetries of the Hamiltonian.

Equations

$$\omega \equiv -\mathrm{d}\eta =-\mathrm{d}(e_i{\pi }^{\ast }{\theta }^i)={\pi }^{\ast }{\theta }^i\wedge \mathrm{d}{e}_i+\frac{1}{2}{e}_i{c}_{jk}^i{\pi }^{\ast }{\theta }^j\wedge {\pi }^{\ast }{\theta }^k,$$

$${\omega }^n\propto \mathrm{d}{e}_1\wedge \cdots \wedge \mathrm{d}{e}_n\wedge {\pi }^{\ast }{\theta }^1\wedge \cdots \wedge {\pi }^{\ast }{\theta }^n$$

$${\omega }_{\mathcal{L}}\equiv \mathbb{F}{\mathcal{L}}^{\ast }\omega =-\mathrm{d}(g_{ij}{v}^j{\pi }^{\ast }{\theta }^i)=g_{ij}{\pi }^{\ast }{\theta }^i\wedge \mathrm{d}{v}^j+\frac{1}{2}{g}_{ir}{v}^r{c}_{jk}^i{\pi }^{\ast }{\theta }^j\wedge {\pi }^{\ast }{\theta }^k.$$

$$\hat{f}=g^{jk}\xi_{k}(f)e_{j}+\Big{(}g^{ik}g^{lj}v^{r}c_{rji}\xi_{l}(f)-g^{jk}e_{j}(f)\Big{)}\xi_{k}.$$

$$\hat{T}=v^k{e}_k+g^{ik}{v}^r{v}^s{c}_{rsi}{\xi }_k\hat{V}=-g^{jk}{e}_j(V){\xi }_k$$

$$q(t)=q(0)v(t)=v(0)-t{g}^{jk}{e}_j(V){|}_{q(0)}{\xi }_k.$$

$$⟨{\mathrm{ad}}_v^{\ast }v,{\xi }_k⟩$$

$$=v^i{v}^j⟨{\xi }_j,c_{ik}^r{\xi }_r⟩=c_{ik}^r{v}^i{v}^j{g}_{jr}=v^i{v}^j{c}_{jik}.$$

$$\big{(}\dot{q},\dot{v}\big{)}=\big{(}v^{k}e_{k},\mathop{\rm ad}\nolimits^{*}_{v}v\big{)}|_{\gamma_{T}(t)}.$$

$$\dot{q}=v^j{e}_j(q)={\partial }_1{L}_{q}v\dot{v}={\mathrm{ad}}_v^{\ast }v-g^{jk}{e}_j(V)(q){\xi }_k,$$

$$q(t)=q(0)v(t)=v(0)-t{\nabla }_{q(0)}V,$$

$$\big{(}q(t),v(t)\big{)}=\big{(}\exp_{q(0)}(v(0)),\mathbb{P}_{q(t)}v(0)\big{)},$$

$${\nabla }_{g}V=g^{jk}(g){\mathrm{d}}_{g}V(e_k)e_j(g)=g^{jk}(1)e_k(V)(g){\partial }_1{L}_g{\xi }_j$$

$$e_{i}(V)(g)=\mathrm{d}(W\circ\iota)(g\cdot\xi_{i})=\frac{\partial W}{\partial{x^{jk}}}\big{(}\partial(\iota\circ L_{g})\xi_{i})\big{)}_{jk}=\mathop{\rm tr}\Big{(}\frac{\partial W}{\partial{x}}^{T}\cdot g\cdot\xi_{i}\Big{)},$$

$$v^{i}(t)\xi_{i}=v^{i}(0)\xi_{i}-tg^{ij}(1)\mathop{\rm tr}\Big{(}\frac{\partial W}{\partial{x}}^{T}\cdot g\cdot\xi_{j}\Big{)}\xi_{i}.$$

$$\gamma_{T}(t)=\big{(}q(0)e^{tv_{0}},v_{0}).$$

$$\dot{v}=-[v_{\mathfrak{k}},v_{\mathfrak{p}}]+{\mathrm{ad}}_{v_{\mathfrak{p}}}^{\ast }{v}_{\mathfrak{k}}.$$

$$\big{(}q(t),v(t)\big{)}=\Big{(}q(0)e^{\big{(}v_{\mathfrak{{p}}}(0)-v_{\mathfrak{{k}}}(0)\big{)}t}e^{2v_{\mathfrak{{k}}}(0)t},v_{\mathfrak{{k}}}(0)+\mathop{\rm Ad}\nolimits e^{-2v_{\mathfrak{{k}}}(0)t}(v_{\mathfrak{{p}}}(0))\Big{)}.$$

$$\dot{q}=v^j{e}_j(q)\dot{v}={\mathrm{ad}}_v^{\ast }v-g^{jk}{e}_j(V)(q){\xi }_{k}v(t)\in \mathfrak{p},$$

$$\dot{q}=v^j{e}_j(q)\dot{v}=-g^{jk}{e}_j(V)(q){\xi }_{k}v(0)\in \mathfrak{p},$$

$$(g,p)=(g,p_{i}\theta^{i}_{g})\mapsto(g,p_{i}\theta^{i}_{1})\mapsto\big{(}g,(p_{1},\ldots,p_{n})\big{)}$$

$$\big{(}F^{*}\tilde{p}_{i}\big{)}(g,\alpha_{g})=(\tilde{p}_{i}\circ F)(g,\alpha_{g})=\tilde{p}_{i}\circ(F(g,\alpha_{g}))=\tilde{p}_{i}\circ(g,L^{*}_{g}\alpha_{g})=L^{*}_{g}\alpha_{g}(e_{i}(1))=p_{i}(\alpha_{g})$$

$$F^{\ast }({\pi }^{\ast }{\theta }^i)={\theta }^i\circ \partial \pi \circ \partial F={\theta }^i\circ \partial (\pi \circ F)={\theta }^i\circ \partial \pi ={\pi }^{\ast }{\theta }^i.$$

$$e_i\circ \mathbb{F}\mathcal{L}(h,v_h)$$

$$\displaystyle=v^{j}_{h}\langle e_{j},e_{i}\rangle(h)=\langle e_{j},e_{i}\rangle(h)\theta^{j}_{h}(v_{h})=\big{(}\theta^{j}\pi_{T\mathcal{G}}^{*}\langle e_{j},e_{i}\rangle\big{)}(h,v_{h}),$$

$$\mathbb{F}\mathcal{L}^{*}\big{(}e_{i}\pi^{*}\theta^{i}\big{)}=\big{(}e_{i}\circ\mathbb{F}\mathcal{L}\big{)}\mathbb{F}\mathcal{L}^{*}\pi^{*}\theta^{i}=\big{(}v^{j}\pi_{T\mathcal{G}}^{*}\langle e_{j},e_{i}\rangle\big{)}\pi^{*}_{T\mathcal{G}}\theta^{i}=v^{j}\pi^{*}_{T\mathcal{G}}\big{(}\langle e_{j},e_{i}\rangle\theta^{i}\big{)}$$

$$f^{\ast }E=E\circ f=V\circ \pi \circ f+\frac{1}{2}⟨\cdot ,\cdot ⟩\circ f=V\circ \pi +\frac{1}{2}⟨f(\cdot ),f(\cdot )⟩$$

$$⟨f(g,v^i{e}_1(1)),f(g,v^j{e}_j(1))⟩$$

$$=L_g^{\ast }⟨v^i{e}_i(1),u^j{e}_j(1)⟩$$

$$f^{\ast }E=V\circ \pi +\frac{1}{2}⟨\cdot ,\cdot {⟩}_1.$$

$${\pi }^{\ast }{\theta }^i\wedge \mathrm{d}{v}^j={\pi }^{\ast }{\theta }^i\otimes \mathrm{d}{v}^j-\mathrm{d}{v}^j\otimes {\pi }^{\ast }{\theta }^i$$

$$g_{ij}\big{(}X^{i}\bar{Y}^{j}-Y^{i}\bar{X}^{j}\big{)}+{\mathchoice{\textstyle{1\over 2}}{\scriptstyle{1\over 2}}{\scriptscriptstyle{1\over 2}}{1/2}}v^{r}c_{rjk}\big{(}X^{j}Y^{k}-Y^{j}X^{k}\big{)}=Y^{i}e_{i}(f)+\bar{Y}^{i}\partial_{v^{i}}f$$

Full text

1 Introduction

The HMC algorithm [14] is a method that generate samples from a probability distribution which has a density known up to constant factor with respect to a reference measure [7, 4, 6]. The method was originally extended to compact Lie groups in [18], and Poisson backets, shadow Hamiltonians and higher order integrators were discussed in [19]. The algorithm was designed by constructing a symplectic structure (and thus a mechanics) on $\mathcal{G}\times\mathfrak{{g}}^{*}$, but the relation with the standard mechanics induced by the Liouville form on $T^{*}\mathcal{G}$ was never explained. One of the difficulties that arise when sampling distributions on manifolds using HMC is the requirement to compute the geodesic flow accurately in order to maintain the motion on the manifold [11, 15, 5].

In the context of Lie groups, the geodesic motion was handled in [18, 19, 13, 20] by assuming the Killing form defined a positive-definite Riemannian metric, in which case geodesics are given by 1-parameter subgroups. This holds for $SO(n)$ and $SU(n)$ but is not case for most other Lie groups, such as $\mathop{\rm GL}(n)$ or $SL(n)$.111In [19] the equations of motion are given for semisimple lie algebras, but the derivation uses the Killing form which is not Riemannian in general and as a result cannot be used for HMC (due to the momentum refreshment step).

In this article we formalise the HMC schemes given in [19, 18] and connect them to standard geometric mechanics on Lie groups, which allows us to remove the compactness assumption on the Lie algebra (or equivalently, the assumption that the Killing form is negative definite).

In §2.2 we show the symplectic structure used in [19, 18] is symplectomorphic to the canonical one on $T^{*}\mathcal{G}$. It follows that many results from geometric mechanics that assume standard mechanics, such as symplectic reduction, can be used in HMC [1]. Since HMC is usually implemented on the tangent bundle, we derive the corresponding symplectic structure on $\mathcal{G}\times\mathfrak{{g}}$ induced by a metric. In [19, 18] the kinetic flow was simply given by left-invariant vector fields as a result of the bi-invariance of the Killing-form. This corresponds to solving the reconstruction equation for a constant curve on $\mathfrak{{g}}$, but, as we show in §2.3, the geodesic flow of more general left-invariant metrics must solve both the Euler-Arnold and reconstruction equations. Hence in §3 we consider the Euler-Arnold equation and explain how it can be simplified and solved by choosing inner products with appropriate symmetries. In particular we derive explicit HMC schemes for $\mathop{\rm GL}(n)$ and semi-simple Lie groups. In §2.4 we provide the explicit relation between the leapfrog scheme on Lie groups as used in [5] and the general one on manifolds; while §2.5 discusses more efficient integrators. In §4 we explain how more generally the mechanical system on homogeneous manifolds can be derived from a constrained mechanics defined on Lie groups. For reductive homogeneous manifolds, the constraints can be naturally handled by choosing a Hamiltonian with sufficient symmetries. Hence we recover the HMC scheme proposed in [5] to sample distributions on naturally reductive homogeneous spaces using Lie group mechanics and symplectic reduction. We will see that this is still true for the discretised motion defined by the leapfrog method or force-gradient integrators (i.e., these preserve the symmetries).

2 Hamiltonian Monte Carlo on Lie Groups

2.1 Lie Groups

Let $L_{g}:\mathcal{G}\to\mathcal{G}$ denote the left translation $L_{g}h=gh$ on a Lie group $\mathcal{G}$, and $\mathfrak{{g}}$ (and $\mathfrak{{g}}_{L}$) denote the lie algebra (of left-invariant vector fields) (see §6.10 for more details). The Maurer-Cartan form $\theta\in\Gamma(T^{*}\mathcal{G}\otimes\mathfrak{{g}})$ 222$\Gamma(\mathcal{M})$ denotes the set of smooth sections of a bundle $\mathcal{M}$. If the base space $\mathcal{B}$ is ambiguous we write $\Gamma(\mathcal{B},\mathcal{M})$. is defined by $\theta:g\mapsto\partial_{g}L_{g^{-1}}$ where $\partial$ denotes the tangent map. Let $\xi_{i}$ be a basis of $\mathfrak{{g}}$, $e_{i}\in\mathfrak{{g}}_{L}$ be the induced left-invariant vector fields (so $e_{i}(1)=\xi_{i}$), and $\theta^{i}$ the dual 1-forms, $\theta^{i}(e_{j})=\delta^{i}_{j}$. We will adopt the Einstein notation and sum over repeated indices.

2.2 Symplectic Structures on Lie Groups

In this section we derive the symplectic structures on $\mathcal{G}\times\mathfrak{{g}}^{*}$, $T\mathcal{G}$ and $\mathcal{G}\times\mathfrak{{g}}$ induced by the canonical 1-form on $T^{*}\mathcal{G}$.

We define the momentum as p\equiv\big{(}\theta^{-1}\big{)}^{*}:T^{*}\mathcal{G}\to\mathfrak{{g}}^{*}. In a basis $e_{i}$ of $\mathfrak{{g}}_{L}$ we may expand these as $\theta=e_{i}(1)\otimes\theta^{i}$ and $p=\theta^{i}(1)\otimes e_{i}$ where we use the canonical isomorphism $T_{g}\mathcal{G}\cong T^{**}_{g}\mathcal{G}$ to view $e_{i}$ as a function $T^{*}\mathcal{G}\to\mathbb{R}$, i.e., e_{i}:(g,\alpha_{g})\mapsto\alpha_{g}\big{(}e_{i}(g)\big{)}. We can then define the 1-form $\eta\equiv e_{i}\pi^{*}\theta^{i}\in\Gamma(T^{*}\mathcal{G},T^{*}T^{*}\mathcal{G})$ on $T^{*}\mathcal{G}$ (here $\pi:T^{*}\mathcal{G}\to\mathcal{G}$), and it is easy to check it yields the canonical Liouville form [5]. The Liouville form may also be expressed in a frame-independent manner as p\big{(}\pi^{*}\theta\big{)}, as well as in terms of the right Maurer Cartan form [2]. The symplectic structure is its exterior derivative333Here we use the fact that $\theta$ is a flat gauge-field, and thus satisfy Maurer-Cartan equation 0=\mathrm{d}\theta+{\mathchoice{\textstyle{1\over 2}}{\scriptstyle{1\over 2}}{\scriptscriptstyle{1\over 2}}{1/2}}[\theta\wedge\theta]=e_{k}(1)\otimes\big{(}\mathrm{d}\theta^{k}+{\mathchoice{\textstyle{1\over 2}}{\scriptstyle{1\over 2}}{\scriptscriptstyle{1\over 2}}{1/2}}c^{k}_{ij}\theta^{i}\wedge\theta^{j}\big{)}.

[TABLE]

and a similar expression holds for any symplectic fiber bundle $\pi:(E,\omega)\to\mathcal{M}$ over a parallelisable manifold $\mathcal{M}$ (i.e., manifolds for which $T\mathcal{M}$ is trivial) [2]. The symplectic volume form is then given by the product of the Haar measure $\theta^{1}\wedge\cdots\wedge\theta^{n}$ on $\mathcal{G}$ with the Lebesgue measure on the fibres

[TABLE]

Originally HMC on compact matrix groups was constructed on the trivial bundle $\mathcal{G}\times\mathfrak{{g}}^{*}$, and although the momentum coordinates were left undefined,444Technically in [19] the symplectic structure is defined as $\mathrm{d}p$ for $p\in T_{g}^{*}\mathcal{G}$ which does not define a differential form over $T^{*}\mathcal{G}$ since the components $p_{i}$ are then constant. we can recover the original mechanical system by defining $\tilde{p}_{i}:\mathcal{G}\times\mathfrak{{g}}^{*}\to\mathbb{R}$ s.t. $\tilde{p}_{i}:(g,\alpha_{1})\to\alpha_{1}(\xi_{i})$, from which we can form the symplectic structure $\tilde{\omega}=-\mathrm{d}(\tilde{p}_{i}\pi^{*}\theta^{i})$. Clearly $\omega$ and $\tilde{\omega}$ are symplectomorphic with respect to the canonical isomorphism $T^{*}\mathcal{G}\cong\mathcal{G}\times\mathfrak{{g}}^{*}$, $(g,\alpha_{g})\mapsto(g,L^{*}_{g}\alpha_{g})$ (see §6.1).

In practice however we usually work on the tangent bundle. Consider a Riemannian metric $\langle\cdot,\cdot\rangle$ on $T\mathcal{G}$. The associated musical isomorphism, $\mathbb{F}\mathcal{L}:T\mathcal{G}\to T^{*}\mathcal{G}$, $v\mapsto\langle v,\cdot\rangle$, enable us to pull-back $\omega$ to a symplectic structure on $T\mathcal{G}$.

The symplectic structure $\omega_{\mathcal{L}}$ induced on $T\mathcal{G}$ by the musical isomorphism is the exterior derivative of the 1-form -v^{j}\pi^{*}\big{(}\langle e_{i},e_{j}\rangle\theta^{i}\big{)}\in\Gamma(T\mathcal{G},T^{*}T\mathcal{G}), where $v^{j}:T\mathcal{G}\to\mathbb{R}$ satisfies $(g,u_{g})\mapsto\theta^{i}_{g}(u_{g})$, and now $\pi:T\mathcal{G}\to\mathcal{G}$ is the tangent bundle projection (see §6.2). To take advantage of the symmetries of $\mathcal{G}$, we assume $\langle\cdot,\cdot\rangle$ is left-invariant and set $g_{ij}\equiv\langle\xi_{i},\xi_{j}\rangle_{1}$. Then

[TABLE]

As above, $\omega_{\mathcal{L}}$ can be pulled-back to define an equivalent mechanics on $\mathcal{G}\times\mathfrak{{g}}$ (the coordinate functions $v^{j}$ are now the maps $(g,X)\mapsto\theta^{j}_{1}(X)$). In particular, if the metric is left-invariant the Hamiltonian $H\equiv V\circ\pi+{\mathchoice{\textstyle{1\over 2}}{\scriptstyle{1\over 2}}{\scriptscriptstyle{1\over 2}}{1/2}}\langle\cdot,\cdot\rangle$ on $T\mathcal{G}$ defines an equivalent Hamiltonian $V\circ\pi+{\mathchoice{\textstyle{1\over 2}}{\scriptstyle{1\over 2}}{\scriptscriptstyle{1\over 2}}{1/2}}\langle\cdot,\cdot\rangle_{1}$ on $\mathcal{G}\times\mathfrak{{g}}$ (see §6.3).

2.3 Hamiltonian Fields and Euler-Arnold Equation

We now focus on the mechanical system on $\mathcal{G}\times\mathfrak{{g}}$ induced by a left-invariant metric $T=v^{i}v^{j}g_{ij}:\mathfrak{{g}}\to\mathbb{R}$. For any function $f:\mathcal{G}\times\mathfrak{{g}}\to\mathbb{R}$, its Hamiltonian vector field $\hat{f}\in\Gamma(T\mathcal{G}\oplus T\mathfrak{{g}})$ is defined by $\omega_{\mathcal{L}}(\hat{f},\cdot)=\mathrm{d}f$. We have (see §6.4)

[TABLE]

For the potential energy and kinetic energy induced by a left-invariant metric

[TABLE]

The Poisson brackets associated to $\omega_{\mathcal{L}}$ are defined by $\{r,h\}\equiv\omega_{\mathcal{L}}(\hat{r},\hat{h})$ for any functions $r,h:\mathcal{G}\times\mathfrak{{g}}\to\mathbb{R}$, and describe the rate of change of $r$ along the flow of $h$ since $\{r,h\}=\omega_{\mathcal{L}}(\hat{r},\hat{h})=\mathrm{d}r(\hat{h})=\hat{h}[r]$. This is often written $\dot{r}=\{r,h\}$ where $\dot{r}:\mathcal{G}\times\mathfrak{{g}}\to\mathbb{R}$ is the function giving the rate of change of $r$ in the direction $\hat{h}$, $\dot{r}(u)\equiv\frac{\mathrm{d}r(\Phi^{\hat{h}}_{u}(t))}{\mathrm{d}t}|_{t=0}$ (here $\Phi^{\hat{h}}:\mathbb{R}\times(\mathcal{G}\times\mathfrak{{g}})\to\mathcal{G}\times\mathfrak{{g}}$ denotes the flow of $\hat{h}$). It follows immediately that a function does not vary along its Hamiltonian flow: $\dot{h}=\{h,h\}=0$, and in fact neither does the symplectic structure $\Phi^{\hat{h}}(t)^{*}\omega_{\mathcal{L}}=\omega_{\mathcal{L}}$, which enables the construction of symplectic integrators (see §6.9). It is also common to write the Poisson bracket of a curve $\gamma:[a,b]\to T\mathcal{G}$ with a function $r$ by setting $\dot{\gamma}=\{\gamma,h\}\equiv\hat{h}(\gamma)$. This is a differential equation stating that $\dot{\gamma}$ is the integral curve of $\hat{h}$. In particular, a curve satisfying $\dot{\gamma}=\{\gamma,H\}$ with initial conditions $(g_{0},v_{0})=\gamma(0)$ is known as a solution of Hamilton’s equations.

An integral curve $\gamma_{V}$ of $\hat{V}$ satisfies $\dot{\gamma}_{V}=\hat{V}(\gamma)$. Let \gamma_{V}(t)=\big{(}q(t),v(t)\big{)}. The equation then reads \big{(}\dot{q},\dot{v}\big{)}=\big{(}0,-g^{jk}e_{j}(V)\xi_{k}\big{)}\big{|}_{\big{(}q(t),v(t)\big{)}}, which implies

[TABLE]

For the integral curves of $\hat{T}$, we have \big{(}\dot{q},\dot{v}\big{)}=\big{(}v^{k}e_{k},g^{ik}v^{r}v^{s}c_{rsi}\xi_{k}\big{)}|_{\gamma_{T}(t)}. The mysterious term $g^{ik}v^{r}v^{s}c_{rsi}\xi_{k}$ disappeared in both [19, 18] (and didn’t appear in [20]) since the structure constants of the Killing form are totally anti-symmetric (see §3.0.1). It is in fact nothing else than the $\langle\cdot,\cdot\rangle$-adjoint of the adjoint representation. Recall that the adjoint representation $\mathop{\rm ad}\nolimits:\mathfrak{{g}}\to\mathop{\rm End(\mathfrak{{g}})}$ satisfies $\mathop{\rm ad}\nolimits_{\xi}\equiv\mathop{\rm ad}\nolimits(\xi)=[\xi,\cdot]$. For each $\xi\in\mathfrak{{g}}$ we can define its adjoint $\mathop{\rm ad}\nolimits^{*}_{\xi}$ with respect to the inner product $\langle\mathop{\rm ad}\nolimits^{*}_{\xi}X,Y\rangle=\langle X,\mathop{\rm ad}\nolimits_{\xi}Y\rangle$. It follows that $\mathop{\rm ad}\nolimits^{*}$ is bilinear so

[TABLE]

Moreover \langle\mathop{\rm ad}\nolimits^{*}_{v}v,\xi_{k}\rangle=\langle\big{(}\mathop{\rm ad}\nolimits^{*}_{v}v\big{)}^{s}\xi_{s},\xi_{k}\rangle=\big{(}\mathop{\rm ad}\nolimits^{*}_{v}v\big{)}^{s}g_{sk}. Combining these two equations gives $\mathop{\rm ad}\nolimits^{*}_{v}v=g^{kl}v^{i}v^{j}c_{jik}\xi_{l}$. Hence the integral curves of $\hat{T}$ satisfy

[TABLE]

The equation $\dot{v}=\mathop{\rm ad}\nolimits^{*}_{v}v$ is a first order differential equation on the lie algebra $\mathfrak{{g}}$ known as the Euler-Arnold equation. In order to derive the integral curves of $\hat{T}$ we can proceed in two steps. We first solve the Euler-Arnold equation to give a solution $s:[0,T]\to\mathfrak{{g}}$. The equation for $\dot{q}$ then states that $\dot{q}=v^{k}e_{k}|_{(q(t),s(t))}=s^{k}(t)e_{k}(q(t))=\partial_{1}L_{q(t)}s(t)$, which is a first order ODE for $q$ called the reconstruction equation.

In summary we have shown that the Hamilton’s equation $\dot{\gamma}=\hat{H}(\gamma)$ defined by the symplectic structure on $\mathcal{G}\times\mathfrak{{g}}$ reads

[TABLE]

which can be recognised as the equations arising from the Hamilton-Pontryagin variational principle.

We finally note that the left-invariance of the metric enabled us to obtain two first order differential equations for the geodesic flow, rather than the usual second-order differential equation on a general Riemannian manifold. This is known as the Euler-Poincaré reduction.

2.4 Leapfrog Integrator

On an arbitrary Riemannian manifold $\mathcal{M}$, we can define the leapfrog integrator for the Hamiltonian system $H=V+T$ as a composition of the kinetic $\hat{T}$ and potential $\hat{V}$ fields defined by the canonical symplectic structure on $T\mathcal{M}$ [22]. The potential flow leads to the velocity update

[TABLE]

(here $\nabla$ is the Riemannian gradient) and the kinetic (or geodesic) flow leads to the update

[TABLE]

where $\exp_{q}v$ denotes the Riemann exponential map and $\mathbb{P}_{q}v$ the parallel transport along the geodesic. Let us see the explicit relation between these updates and the ones derived in §2.3. Suppose that $\mathcal{M}=\mathcal{G}$ is a Lie group and $T$ is defined by a left-invariant metric $\langle\cdot,\cdot\rangle$. From §2.2 we know that, under the symplectomorphism $T\mathcal{G}\cong\mathcal{G}\times\mathfrak{{g}}$, the Hamiltonian vector-field induced by the inner product $\langle\cdot,\cdot\rangle_{1}:\mathfrak{{g}}\to\mathbb{R}$ push-forwards to the geodesic flow induced by $T:T\mathcal{G}\to\mathbb{R}$. More precisely if $(q(t),v(t))\in\mathcal{G}\times\mathfrak{{g}}$ is the integral curve of the Hamiltonian vector field of $\langle\cdot,\cdot\rangle_{1}$, then the geodesic flow of the Riemannian metric $T=\langle\cdot,\cdot\rangle$ on $T\mathcal{G}$ is given by \big{(}q(t),\partial_{1}L_{q(t)}v(t)\big{)}=\big{(}q(t),\dot{q}(t)\big{)}.

For the potential term, observe that555Note that if $\alpha\in\Gamma(T^{*}\mathcal{G})$, then $\alpha=e_{i}(\alpha)\theta^{i}=\alpha(e_{i})\theta^{i}=\alpha_{i}\theta^{i}$. Moreover $u=g^{ij}\alpha_{i}e_{j}$ is the vector field associated to $\alpha$ by the metric since $\alpha(v)=\alpha_{i}\theta^{i}(v^{j}e_{j})=\alpha_{i}v^{i}=\alpha_{k}\delta^{k}_{i}v^{i}=\alpha_{k}g^{kr}g_{ri}v^{i}=u^{r}g_{ri}v^{i}=\langle u,v\rangle$.

[TABLE]

(by left-invariance $g^{ij}$ is constant), and $v(t)=\theta^{i}_{g}(v(t))e_{i}(g)=v^{i}(t)\partial_{1}L_{g}\xi_{i}$. Thus under the symplectomorphism $\mathcal{G}\times\mathfrak{{g}}\to T\mathcal{G}$, we have ($\xi_{i}\to\partial_{1}L_{g}\xi_{i})$ and (5) push-forwards to the velocity update (4).

Finally let us see how the velocity update (4) is computed in practice when $\mathcal{G}$ is a matrix group. In that case we need an extension $W$ of $V$ defined on an open subset of $\mathop{\rm GL}(n)$. Then

[TABLE]

where $\cdot$ is matrix multiplication, $\iota:\mathcal{G}\to\mathop{\rm GL}(n)$ is the inclusion (or a more general injective representation), and in the last line we have identified $g\sim\iota(g),\xi_{i}\sim\partial\iota(\xi_{i})$. We arrive at the general velocity update formula for a matrix group

[TABLE]

In lattice gauge theory we typically have W(x)=Re\mathop{\rm tr}\big{(}Ux\big{)} for a constant complex matrix $U$.

2.5 Efficient Integrators

While the leapfrog integrator is the most common choice of integrator, it is not always the optimal choice. For example [27, 26] consider a symplectic integrator of the form666Using the notation $\exp(t\hat{H})\equiv\Phi^{\hat{H}}(t)$ for the Hamiltonian flow \exp\big{(}\lambda\delta t\hat{T}\big{)}\exp\big{(}{\mathchoice{\textstyle{1\over 2}}{\scriptstyle{1\over 2}}{\scriptscriptstyle{1\over 2}}{1/2}}\delta t\hat{V}\big{)}\exp\big{(}(1-2\lambda)\delta t\hat{T}\big{)}\exp\big{(}{\mathchoice{\textstyle{1\over 2}}{\scriptstyle{1\over 2}}{\scriptscriptstyle{1\over 2}}{1/2}}\delta t\hat{V}\big{)}\exp\big{(}\lambda\delta t\hat{T}\big{)} which, despite being computationally more expensive, they show is $50\%$ more efficient than leapfrog. In [19] a “force-gradient” integrator is examined which involves the Hamiltonian flow $\widehat{\{V,\{V,T\}\}}$, and thus contain second-order information (derivatives) about $V$. The question of tuning HMC using Poisson brackets to minimise the cost was considered in [13]; see also [12] for a discussion on Lie group integrators, and [9] for integrators which outperform leapfrog on $\mathbb{R}^{n}$.

3 Geodesics

We now analyse how Euler-Arnold equation $\dot{v}=\mathop{\rm ad}\nolimits^{*}_{v}v$ simplifies when the inner product on $\mathfrak{{g}}$ is $\mathop{\rm ad}\nolimits_{\mathfrak{{k}}}$-invariant for a subalgebra $\mathfrak{{k}}\subset\mathfrak{{g}}$. When $\mathcal{G}$ is connected, $\mathop{\rm ad}\nolimits_{\mathfrak{{k}}}$-invariant inner products correspond to left $\mathcal{G}$-invariant and right $\mathcal{K}$-invariant Riemannian metrics on $\mathcal{G}$. 777i.e., the left and right actions $L_{g}$ and $R_{k}$ are isometries for $g\in\mathcal{G},k\in\mathcal{K}$ (here $\mathcal{K}$ is the Lie group associated to $\mathfrak{{k}}$). When $\mathcal{G}$ is not connected an analogous correspondence exists with $\mathop{\rm Ad}\nolimits_{\mathcal{K}}$-invariant inner products [5], i.e., inner products s.t. $\langle\mathop{\rm Ad}\nolimits_{k}u,\mathop{\rm Ad}\nolimits_{k}v\rangle_{1}=\langle u,v\rangle_{1}$ for $u,v\in\mathfrak{{g}}$, which are automatically $\mathop{\rm ad}\nolimits_{\mathfrak{{k}}}$-invariant. When $\mathfrak{{k}}=\mathfrak{{g}}$ these correspond to bi-invariant Riemannian metric on $\mathcal{G}$, which always exists when $\mathcal{G}$ is compact [25].

3.0.1 $\mathop{\rm Ad}\nolimits\mathcal{G}$-invariant Inner Products

The simplest case arises when $\mathop{\rm ad}\nolimits^{*}_{v}v=0$. This situation occurs for example when the inner product on $\mathfrak{{g}}$ is $\mathop{\rm ad}\nolimits$-invariant, $\langle\mathop{\rm ad}\nolimits_{u}v,w\rangle+\langle v,\mathop{\rm ad}\nolimits_{u}w\rangle=0$, which implies $\mathop{\rm ad}\nolimits^{*}=-\mathop{\rm ad}\nolimits$ and $\mathop{\rm ad}\nolimits^{*}_{v}v=-[v,v]=0$. We can find $\mathop{\rm ad}\nolimits$-invariant inner products on $\mathfrak{{g}}$ when $\mathcal{G}$ has an $\mathop{\rm Ad}\nolimits\mathcal{G}$-invariant inner product, or equivalently if it has a bi-invariant metric. A particular $\mathop{\rm Ad}\nolimits\mathcal{G}$-invariant inner product is available on Lie groups whose lie algebra is compact, which means the Killing-form B(u,v)\equiv\mathop{\rm tr}\big{(}\mathop{\rm ad}\nolimits(u)\mathop{\rm ad}\nolimits(v)\big{)} is negative definite, and so $-B$ defines a positive definite $\mathop{\rm ad}\nolimits$-invariant inner product.

When $\mathop{\rm ad}\nolimits^{*}_{v}v=0$, then $v(t)=v_{0}$, so the reconstruction equation becomes $\dot{q}=v^{k}e_{k}|_{(q(t),v_{0})}=v_{0}^{k}e_{k}(q(t))=\partial_{1}L_{q(t)}v_{0}=\frac{\mathrm{d}}{\mathrm{d}t}L_{q(t)}e^{tv_{0}}|_{t=0}$, which implies $q(t)=q(0)e^{tv_{0}}$ (here we used the fact that $e^{tv_{0}}$ is tangent to $v_{0}$ at $t=0$). Hence we obtain the kinetic flow as in [19]

[TABLE]

3.0.2 Reductive Decomposition

Sometimes it is not possible to find a bi-invariant metric on $\mathcal{G}$. The natural next step is to consider inner products on $\mathfrak{{g}}$ which are $\mathop{\rm ad}\nolimits_{\mathfrak{{k}}}$-invariant, where $\mathfrak{{k}}$ is a subalgebra of $\mathfrak{{g}}$. This means that $\langle\mathop{\rm ad}\nolimits_{u}v,w\rangle+\langle v,\mathop{\rm ad}\nolimits_{u}w\rangle=0$ for all $u\in\mathfrak{{k}}$, and it follows that $\mathop{\rm ad}\nolimits^{*}_{u}=-\mathop{\rm ad}\nolimits_{u}$. Then $\mathcal{K}$ is totally geodesic in $\mathcal{G}$ and its geodesics are given by 1-parameter subgroups (see §6.5). Let $\mathfrak{{p}}\equiv\mathfrak{{k}}^{\perp}$, so $\mathfrak{{g}}=\mathfrak{{k}}\oplus\mathfrak{{p}}$ and $v=v_{\mathfrak{{k}}}+v_{\mathfrak{{p}}}$. The Euler-Arnold equation reads $\dot{v}=\dot{v}_{\mathfrak{{k}}}+\dot{v}_{\mathfrak{{p}}}=\mathop{\rm ad}\nolimits^{*}_{v}v=\mathop{\rm ad}\nolimits^{*}_{v_{\mathfrak{{k}}}}v_{\mathfrak{{k}}}+\mathop{\rm ad}\nolimits^{*}_{v_{\mathfrak{{k}}}}v_{\mathfrak{{p}}}+\mathop{\rm ad}\nolimits^{*}_{v_{\mathfrak{{p}}}}v_{\mathfrak{{k}}}+\mathop{\rm ad}\nolimits^{*}_{v_{\mathfrak{{p}}}}v_{\mathfrak{{p}}}=-[v_{\mathfrak{{k}}},v_{\mathfrak{{p}}}]+\mathop{\rm ad}\nolimits^{*}_{v_{\mathfrak{{p}}}}v_{\mathfrak{{k}}}+\mathop{\rm ad}\nolimits^{*}_{v_{\mathfrak{{p}}}}v_{\mathfrak{{p}}}$. We can simplify this further if the decomposition $\mathfrak{{k}}\oplus\mathfrak{{p}}$ is naturally reductive, i.e., $[\mathfrak{{k}},\mathfrak{{p}}]\subset\mathfrak{{p}}$ and for all $u,v,w\in\mathfrak{{p}}$, $\langle[u,v]_{\mathfrak{{p}}},w\rangle=\langle u,[v,w]_{\mathfrak{{p}}}\rangle$ [3]. This implies $\mathop{\rm ad}\nolimits^{*}_{v_{\mathfrak{{p}}}}v_{\mathfrak{{p}}}=0$ (see §6.6), and we are left with

[TABLE]

Notably the constant curve $v=v_{0}$ with either $v_{0}\in\mathfrak{{k}}$ or $v_{0}\in\mathfrak{{p}}$ is a solution.

3.0.3 Matrix Groups

Consider a subgroup $\mathcal{G}\subset\mathop{\rm GL}(n)$ and set $\mathfrak{{g}}=\mathfrak{{p}}\oplus\mathfrak{{k}}$ where $\mathfrak{{p}}$ is a subspace of the space of symmetric matrices, while $\mathfrak{{k}}=\mathfrak{so}(n)$. We have in mind the cases $\mathop{\rm GL}(n),SL(n)$ and $SO(n)$ with $\mathfrak{{p}}$ the space of symmetric matrices, the space of traceless symmetric matrices, and $\{0\}$ respectively. We equip $\mathfrak{{g}}$ with the standard inner product $\langle A,B\rangle\equiv\mathop{\rm tr}(A^{T}B)$ which is $\mathop{\rm Ad}\nolimits\mathcal{K}$-invariant. Since the product of a symmetric and an antisymmetric matrix is traceless, we have $\langle A,B\rangle=\mathop{\rm tr}(A_{\mathfrak{{p}}}B_{\mathfrak{{p}}})-\mathop{\rm tr}(A_{\mathfrak{{k}}}B_{\mathfrak{{k}}})$.

Note that for $\mathcal{G}=SO(n)$, this is just the negative of the Killing form (up to a positive constant); while on $\mathop{\rm GL}(n)$ and $SL(n)$, the Killing form, respectively given by $2n\mathop{\rm tr}(AB)-2\mathop{\rm tr}(A)\mathop{\rm tr}(B)$ and $2n\mathop{\rm tr}(AB)$, doesn’t define an inner product (it is degenerate and pseudo-Riemannian respectively [23]).

For our inner product we have $\mathop{\rm ad}\nolimits^{*}_{S}A=[S,A]$ for any $S\in\mathfrak{{p}},A\in\mathfrak{{k}}$, hence Euler-Arnold equation simplifies to $\dot{v}=-2[v_{\mathfrak{{k}}},v_{\mathfrak{{p}}}]$ with solution (see §6.7)

[TABLE]

3.0.4 Semi-simple Symmetric Spaces

The equation for the geodesics above also holds for a semi-simple Lie group $\mathcal{G}$ using an inner product based on the Killing form $B$ (note $\mathop{\rm GL}(n)$ is not semi-simple). Let us consider a subgroup $\mathcal{K}$ for which $\mathcal{G}/\mathcal{K}$ is a symmetric space. In that case we can find a Cartan decomposition $\mathfrak{{g}}=\mathfrak{{k}}\oplus\mathfrak{{p}}$. By corollary 5.4.3 [17] the Killing form is negative-definite on $\mathcal{K}$. We say $\mathcal{G}/\mathcal{K}$ is of compact (non-compact) type if $B$ is negative (positive) definite on $\mathfrak{{p}}$. If $\mathcal{G}/\mathcal{K}$ is of compact type the inner product $\langle v,u\rangle\equiv B(v_{\mathfrak{{k}}},u_{\mathfrak{{k}}})+B(v_{\mathfrak{{p}}},u_{\mathfrak{{p}}})$ is $\mathop{\rm Ad}\nolimits\mathcal{G}$-invariant, while if it is of non-compact type $\langle v,u\rangle\equiv-B(v_{\mathfrak{{k}}},u_{\mathfrak{{k}}})+B(v_{\mathfrak{{p}}},u_{\mathfrak{{p}}})$ is $\mathop{\rm Ad}\nolimits\mathcal{K}$-invariant. Then $\mathop{\rm ad}\nolimits^{*}_{v_{\mathfrak{{p}}}}v_{\mathfrak{{k}}}=[v_{\mathfrak{{p}}},v_{\mathfrak{{k}}}]$ (see Lemma 5.5.4 [17]) and we can find the integral curves of $\hat{T}$ as in §3.0.3.

3.0.5 Approximating the Matrix Exponential

In practice there are several ways to approximate the matrix exponential, for example by combining a Padé approximant with a projection from $\mathop{\rm GL}(n)$ to $\mathcal{G}$ [8]. In particular the Cayley transform $z\mapsto(1+{\mathchoice{\textstyle{1\over 2}}{\scriptstyle{1\over 2}}{\scriptscriptstyle{1\over 2}}{1/2}}z)/(1-{\mathchoice{\textstyle{1\over 2}}{\scriptstyle{1\over 2}}{\scriptscriptstyle{1\over 2}}{1/2}}z)=(1-{\mathchoice{\textstyle{1\over 2}}{\scriptstyle{1\over 2}}{\scriptscriptstyle{1\over 2}}{1/2}}z)^{-1}(1+{\mathchoice{\textstyle{1\over 2}}{\scriptstyle{1\over 2}}{\scriptscriptstyle{1\over 2}}{1/2}}z)$, which is the Padé-$(1,1)$ approximant, maps exactly onto $\mathcal{G}$ for quadratic groups888By a quadratic group (sometimes called $J$-orthogonal) we mean a subgroup $\mathcal{G}=\{M\in\mathop{\rm GL}(n):M^{T}JM=J\}$ for some $J\in\mathop{\rm GL}(n)$(in fact diagonal Padé approximants are product of Cayley transforms) [16] .

4 Constrained Mechanics on Homogeneous Manifolds

It was shown in [5] that HMC on naturally reductive manifolds can be implemented using mechanics on Lie groups. Consider a Homogeneous manifold $\mathcal{G}/\mathcal{K}$ (where $\mathcal{K}$ is a closed subgroup of $\mathcal{G}$) and a Hamiltonian $\tilde{H}=\tilde{V}+\tilde{T}$ on T\big{(}\mathcal{G}/\mathcal{K}\big{)}. Let $\pi:\mathcal{G}\to\mathcal{G}/\mathcal{K}$ be the quotient projection. We can lift $\tilde{H}$ to $T\mathcal{G}$ by defining $H\equiv(\partial\pi)^{*}\tilde{H}$. While $V\equiv\pi^{*}\tilde{V}$ is well-defined on $\mathcal{G}$, the lift $T\equiv(\partial\pi)^{*}\tilde{T}$ defines a degenerate kinetic energy on $T\mathcal{G}$, since it maps all vectors in the vertical space $\text{ver}_{g}\equiv\text{ker}(\partial_{g}\pi)$ to zero. The vertical space is described by the action fields $\text{ver}_{g}=\{\xi_{\mathcal{G}}(g):\xi\in\mathfrak{{k}}\}$ where $\xi_{\mathcal{G}}(g)=\frac{d}{dt}e^{t\xi}\cdot g|_{t=0}=\partial_{1}L_{g}(\xi)$.

In order to circumvent the issues related to the degeneracy of $T$, we need to make a choice of horizontal space complementary to $\text{ver}_{g}$, or equivalently choose a connection on the principal bundle $\mathcal{G}\to\mathcal{G}/\mathcal{K}$. This defines a unique way to lift a curve on $\mathcal{G}/\mathcal{K}$ to a curve on $\mathcal{G}$, and a complementary subspace $\mathfrak{{p}}$ to $\mathfrak{{k}}$, i.e., $\mathfrak{{g}}=\mathfrak{{k}}\oplus\mathfrak{{p}}$ (see §6.11). Then, if $(q,v)$ satisfies the constrained equation

[TABLE]

the projected curve $\pi(q(t))$ solves Hamilton’s equations on $\mathcal{G}/\mathcal{K}$ [21]. We note the constraint $v(t)\in\mathfrak{{p}}$ can be rewritten in the form ver\big{(}v(t)\big{)}=0 where $ver:\mathfrak{{g}}\to\mathfrak{{k}}$ is the vertical projection.

If $\mathcal{G}/\mathcal{K}$ is a reductive homogeneous manifold, there is a one-to-one correspondance between $\mathop{\rm Ad}\nolimits_{\mathcal{G}}\mathcal{K}$-invariant inner products on $\mathfrak{{p}}$ and $\mathcal{G}$-invariant metric on $\mathcal{G}/\mathcal{K}$. In particular if we extend the inner product on $\mathfrak{{p}}$ to a non-degenerate quadratic form on $\mathfrak{{g}}$ s.t. $\mathfrak{{p}}=\mathfrak{{k}}^{\perp}$, it defines a pseudo-Riemannian metric on $\mathcal{G}$ that is left $\mathcal{G}$-invariant and right $\mathcal{K}$-invariant, and a kinetic energy $T$ with $T=(\partial\pi)^{*}\tilde{T}$ on $\mathfrak{{p}}$. The resulting Hamiltonian $H=V+T$ is right $\mathcal{K}$-invariant so by Noether’s theorem, if $v(0)\in\mathfrak{{p}}$, then $v(t)\in\mathfrak{{p}}$ for all $t$. This means the constraint $v(t)\in\mathfrak{{p}}$ is naturally handled by the symmetries of $H$ (see §6.8).

If $\mathcal{G}/\mathcal{K}$ is further naturally reductive, since $\mathop{\rm ad}\nolimits^{*}_{v}v=0$ for $v\in\mathfrak{{p}}$ the above system becomes

[TABLE]

and we recover the HMC algorithm on naturally reductive homogeneous spaces proposed in [5].

It is clear that integrators that are built by alternating between the flows of $V$ and $T$ are automatically horizontal (the trajectory stays in $\mathcal{G}\times\mathfrak{{p}}$) since both $V$ and $T$ are $\mathcal{K}$-invariant. This is still true for the force-gradient integrators §2.5. Indeed, from the $\mathop{\rm ad}\nolimits_{\mathfrak{{k}}}$-invariance of the inner product on $\mathfrak{{g}}$, we have $c_{rij}=-c_{jir}$ for any $i\in\mathfrak{{k}}$ (the structure constants always satisfy $c^{i}_{jk}=-c^{i}_{kj}$). Hence if $V$ is invariant under the right action of $\mathcal{K}$ on $\mathcal{G}$, then $\widehat{\{V,\{V,T\}\}}=-2g^{jk}g^{ls}e_{l}(V)e_{j}e_{s}(V)\xi_{k}\in\mathfrak{{p}}$ (see §6.12).

5 Acknowledgements

We thank Prof. Anthony Kennedy and Prof. Mark Girolami for useful discussions. AB was supported by a Roth Scholarship funded by the Department of Mathematics, Imperial College London, by EPSRC Fellowship (EP/J016934/3) and by The Alan Turing Institute under the EPSRC grant [EP/N510129/1].

6 Supplementary Material

6.1 Symplectomorphism $T^{*}\mathcal{G}\cong\mathcal{G}\times\mathfrak{g}^{*}$

The cotangent bundle over $\mathcal{G}$ is isomorphic to a trivial bundle $T^{*}\mathcal{G}\cong\mathcal{G}\times\mathfrak{{g}}^{*}\cong\mathcal{G}\times\mathbb{R}^{n}$. Indeed fix a frame $e_{i}\in\mathfrak{{g}}_{L}$ and co-frame $\theta^{i}\in\mathfrak{{g}}^{*}_{L}$. We have the maps

[TABLE]

Letting $F:T^{*}\mathcal{G}\to\mathcal{G}\times\mathfrak{{g}}^{*}$ be this diffeomorphism. Then

[TABLE]

and

[TABLE]

Hence $T^{*}\mathcal{G}$ with canonical symplectic structure is symplectomorphic to $\mathcal{G}\times\mathfrak{{g}}^{*}$ with the symplectic structure from [19, 18].

6.2 Symplectic Structure on $T\mathcal{G}$

The pull-back of the coordinate functions $e_{i}$ is $e_{i}\circ\mathbb{F}\mathcal{L}:T\mathcal{G}\to\mathbb{R}$, and

[TABLE]

where $\pi_{T\mathcal{G}}:T\mathcal{G}\to\mathcal{G}$ is the canonical projection. Note in the final expression $\theta^{j}$ is viewed as a map $T\mathcal{G}\to\mathbb{R}$ (which we denote by $v^{j}$) rather than an element of $\Gamma(T^{*}\mathcal{G})$. This gives the definition of the coordinate functions $v^{j}$ on $T\mathcal{G}$. On the other hand, $\mathbb{F}\mathcal{L}^{*}\pi^{*}\theta^{i}=(\pi\circ\mathbb{F}\mathcal{L})^{*}\theta^{i}$ and $\pi\circ\mathbb{F}\mathcal{L}:T\mathcal{G}\to\mathcal{G}$ is simply $\pi_{T\mathcal{G}}$.

In summary,

[TABLE]

6.3 Symplectomorphism $T\mathcal{G}\cong\mathcal{G}\times\mathfrak{{g}}$

The isomorphism $T\mathcal{G}\cong\mathcal{G}\times\mathfrak{{g}}$ can be written as $(g,v)\mapsto(g,\theta_{g}(v))=(g,\partial_{g}L_{g^{-1}}v)$, where $\theta$ is the (left) Maurer-Cartan form. We can thus pull-back our Hamiltonian system from $T\mathcal{G}$ to $\mathcal{G}\times\mathfrak{{g}}$: let $f:\mathcal{G}\times\mathfrak{{g}}\to T\mathcal{G}$ be this map, and let $E=V\circ\pi+{\mathchoice{\textstyle{1\over 2}}{\scriptstyle{1\over 2}}{\scriptscriptstyle{1\over 2}}{1/2}}\langle\cdot,\cdot\rangle$ be an energy function on $T\mathcal{G}$. Then

[TABLE]

and

[TABLE]

Thus if the metric is left-invariant

[TABLE]

6.4 Hamiltonian Vector Fields

We proceed similarly to [19] but work on $\mathcal{G}\times\mathfrak{{g}}$ which will enable us to identify the coordinate vector fields $\partial_{v^{i}}$ with $\xi_{i}$. An arbitrary vector field $Y$ on $\mathcal{G}\times\mathfrak{{g}}$ may be expanded as $Y=Y^{i}e_{i}+\bar{Y}^{i}\partial_{v^{i}}$, where $v^{i}:(g,u_{1})\mapsto\theta^{i}_{1}(u_{1})$ are the fibre coordinates. This follows because any vector tangent to $\mathcal{G}\times\mathfrak{{g}}$ is given by the derivative of a curve $\gamma:[a,b]\to\mathcal{G}\times\mathfrak{{g}}$ which may be written as $\gamma(t)=(q(t),u(t))$; then $\dot{\gamma}=(\dot{q},\dot{u})$, where $\dot{q}\in T\mathcal{G}$ and $\dot{u}\in T\mathfrak{{g}}$ can be expanded in terms of the (global trivialising) frames $e_{i}$ and $\partial_{v^{i}}$.

The equation $\omega_{\mathcal{L}}(\hat{f},Y)=\mathrm{d}f(Y)$ must hold for every $Y$ if $\hat{f}$ is the Hamiltonian vector field of $f$. Note

[TABLE]

and set $\hat{f}=X^{i}e_{i}+\bar{X}^{i}\partial_{v^{i}}$ Then $\pi^{*}\theta^{i}\wedge\mathrm{d}v^{j}(\hat{f},Y)=X^{i}\bar{Y}^{j}-Y^{i}\bar{X}^{j}$, moreover $\pi^{*}\theta^{i}\wedge\pi^{*}\theta^{k}(\hat{f},Y)=X^{j}Y^{k}-Y^{j}X^{k}$, and finally $\mathrm{d}f(Y)=Y(f)=Y^{i}e_{i}(f)+\bar{Y}^{i}\partial_{v^{i}}f$. Hence $\omega_{\mathcal{L}}(\hat{f},Y)=\mathrm{d}f(Y)$ is equivalent to

[TABLE]

where $c_{rjk}\equiv g_{rl}c^{l}_{jk}$. As this must hold for every $Y$, equating the coefficients of $Y^{i}$ and $\bar{Y}^{i}$ yields $g_{ij}X^{i}=\partial_{v^{j}}f$ and $g_{ij}\bar{X}^{j}=v^{r}c_{rji}X^{j}-e_{i}(f)$ respectively, where we have used $c_{rjk}=-c_{rkj}$ since the Lie bracket of vector field is antisymmetric: $[Z,W]=-[W,Z]$.

It follows that

[TABLE]

In particular if $T={\mathchoice{\textstyle{1\over 2}}{\scriptstyle{1\over 2}}{\scriptscriptstyle{1\over 2}}{1/2}}g_{ij}v^{i}v^{j}$ and $V$ is independent of fibre coordinates $v^{i}$, then

[TABLE]

where we have used $e_{i}[T]=0$ since e_{i}[T]={\mathchoice{\textstyle{1\over 2}}{\scriptstyle{1\over 2}}{\scriptscriptstyle{1\over 2}}{1/2}}e_{i}\big{(}g_{jk}(g)\big{)}v^{j}v^{k}={\mathchoice{\textstyle{1\over 2}}{\scriptstyle{1\over 2}}{\scriptscriptstyle{1\over 2}}{1/2}}e_{i}\big{(}g_{jk}(1)\big{)}v^{j}v^{k}=0 by left-invariance of the metric.

In fact the coordinate vectors $\partial_{v^{k}}$ are simply the generators, $\partial_{v^{k}}=\xi_{k}$, since both of them have the integral curve $\gamma:t\mapsto t\xi_{k}$. Indeed, if $v=(v^{1},\ldots,v^{n}):\mathfrak{{g}}\to\mathbb{R}^{n}$ denotes the coordinate chart $\dot{\gamma}[f]=\frac{\partial(f\circ v^{-1})}{\partial v^{i}}\frac{\partial v^{i}\circ\gamma}{dt}|_{t=0}=\frac{\partial(f\circ v^{-1})}{\partial v^{k}}=\frac{\partial}{\partial v^{k}}(f)$.

Note that for a general $f$, $e_{j}(f)$ is a function on $T\mathcal{G}$, however since $V$ does not depend on the fibre coordinates, $e_{j}(V)$ is just a function of $\mathcal{G}$.

6.5 Totally Geodesic Subgroups

Recall we call a subgroup $\mathcal{K}\subset\mathcal{G}$ totally geodesic, if the geodesics of $\mathcal{K}$ are geodesics of $\mathcal{G}$.

Suppose we have an $\mathop{\rm ad}\nolimits_{\mathfrak{{k}}}$-invariant inner product on $\mathfrak{{g}}$. Let $\mathfrak{{p}}\equiv\mathfrak{{k}}^{\perp}$, then for any $u\in\mathfrak{{k}},v\in\mathfrak{{p}}$, $\langle u,\mathop{\rm ad}\nolimits_{u}v\rangle=-\langle\mathop{\rm ad}\nolimits_{u}u,v\rangle=-\langle[u,u],v\rangle=0$, so Theorem 2 [24] is satisfied. In fact $\mathfrak{{k}}$ is “easy totally geodesic” in $\mathfrak{{g}}$, i.e., $\mathop{\rm ad}\nolimits_{\mathfrak{{k}}}\mathfrak{{p}}\subset\mathfrak{{p}}$, since for any $u,w\in\mathfrak{{k}}$, $v\in\mathfrak{{p}}$ we find $\langle\mathop{\rm ad}\nolimits_{u}v,w\rangle=\langle v,\mathop{\rm ad}\nolimits^{*}_{u}w\rangle=-\langle v,[w,u]\rangle=0$ since $[w,u]\in\mathfrak{{k}}$ and $v\in\mathfrak{{p}}$, hence $\mathop{\rm ad}\nolimits_{u}v\in\mathfrak{{p}}$ (see def 4.1 [24]).

Finally we note that the restriction of the inner product to $\mathfrak{{k}}$ is $\mathop{\rm ad}\nolimits$-invariant and thus the geodesics are given by the 1-parameter subgroups [24].

6.6 Naturally Reductive Geodesics

Recall that we have an $\mathop{\rm ad}\nolimits_{\mathfrak{{k}}}$-invariant inner product, i.e. $\langle\mathop{\rm ad}\nolimits_{u}v,w\rangle+\langle v,\mathop{\rm ad}\nolimits_{u}w\rangle=0$ for all $u\in\mathfrak{{k}}$, and a naturally reductive decomposition implies for $u,v,w\in\mathfrak{{p}}$, $\langle[u,v]_{\mathfrak{{p}}},w\rangle=\langle u,[v,w]_{\mathfrak{{p}}}\rangle$, where $\mathfrak{{p}}=\mathfrak{{k}}^{\perp}$.

We want to show that $\mathop{\rm ad}\nolimits^{*}_{v}v=0$ for any $v\in\mathfrak{{p}}$. Let $z\in\mathfrak{{g}}$. Then

[TABLE]

where we have used

[TABLE]

by orthogonality and the fact $[v,v]=0$.

Moreover $\langle v,\mathop{\rm ad}\nolimits_{z_{\mathfrak{{k}}}}v\rangle=-\langle\mathop{\rm ad}\nolimits_{z_{\mathfrak{{k}}}}v,v\rangle=-\langle v,\mathop{\rm ad}\nolimits_{z_{\mathfrak{{k}}}}v\rangle$, so $\langle v,\mathop{\rm ad}\nolimits_{z_{\mathfrak{{k}}}}v\rangle=0$. In fact this only requires the inner product on $\mathfrak{{g}}$ to be $\mathop{\rm ad}\nolimits_{\mathfrak{{k}}}\mathfrak{{p}}$-invariant, that is $\langle\mathop{\rm ad}\nolimits_{u}v,w\rangle+\langle v,\mathop{\rm ad}\nolimits_{u}w\rangle=0$ for all $u\in\mathfrak{{k}},v,w\in\mathfrak{{p}}$.

Finally let us consider the term $\mathop{\rm ad}\nolimits^{*}_{v_{\mathfrak{{p}}}}v_{\mathfrak{{k}}}$. Note let $S_{i},A_{j}$ basis of $\mathfrak{{p}},\mathfrak{{k}}$ respectively, and $T_{k}$ arbitrary basis element of $\mathfrak{{g}}$. Then $\langle\mathop{\rm ad}\nolimits^{*}_{S_{i}}A_{j},T_{k}\rangle=\langle A_{j},[S_{i},T_{k}]\rangle$. This is zero if $T_{k}\in\mathfrak{{k}}$, since $[\mathfrak{{p}},\mathfrak{{k}}]\subset\mathfrak{{p}}$ and $\mathfrak{{p}}=\mathfrak{{k}}^{\perp}$. Moreover $\langle\mathop{\rm ad}\nolimits^{*}_{S_{i}}A_{j},T_{k}\rangle=\langle A_{j},c^{r}_{ik}T_{r}\rangle=c^{r}_{ik}g_{jr}=c_{jik}$ where $j\in\mathfrak{{k}}$, $i\in\mathfrak{{p}}$. On the other hand $\langle[S_{i},A_{j}],T_{k}\rangle$ is again zero if $T_{k}\in\mathfrak{{k}}$; and $\langle[S_{i},A_{j}],T_{k}\rangle=c^{r}_{ij}g_{rk}=c_{kij}=-c_{kji}=c_{ijk}$ for $i\in\mathfrak{{p}}$, $j\in\mathfrak{{k}}$, where in the last step we used $\mathop{\rm ad}\nolimits_{\mathfrak{{k}}}$-invariance. It doesn’t seem to be possible in general to relate $c_{jik}$ with $c_{ijk}$.

6.7 Geodesics Matrix Groups

Let us first verify that $\mathop{\rm ad}\nolimits^{*}_{S}A=[S,A]$, for $S\in\mathfrak{{p}}$, $A\in\mathfrak{{k}}$. Let $T\in\mathfrak{{g}}$. Then

[TABLE]

and since $T$ is arbitrary the result holds. Here we have used the property of the trace that $\mathop{\rm tr}(ABC)=\mathop{\rm tr}(BCA)$ and the fact that $[S,A]$ is symmetric, [S,A]^{T}=\big{(}SA-AS)^{T}=-AS+SA=[S,A]. It follows that the Euler-Arnold equation becomes $\dot{v}=\dot{v}_{\mathfrak{{k}}}+\dot{v}_{\mathfrak{{p}}}=-2[v_{\mathfrak{{k}}},v_{\mathfrak{{p}}}]$. Recall $\mathop{\rm ad}\nolimits_{u}=\partial_{1}\mathop{\rm Ad}\nolimits(u)=\frac{d}{dt}\mathop{\rm Ad}\nolimits(e^{tu})|_{t=0}$ for any $u\in\mathfrak{{g}}$.

For any $u,w\in\mathfrak{{g}}$ the curve $t\mapsto\mathop{\rm Ad}\nolimits(e^{tu})(w)$ is tangent to $[u,\mathop{\rm Ad}\nolimits(e^{tu})(w)]$ at $t$, since (for matrix groups) $\mathop{\rm Ad}\nolimits(e^{tu})w=e^{tu}we^{-tu}$ and differentiating gives $ue^{tu}we^{-tu}-e^{tu}we^{-tu}u=u\cdot\mathop{\rm Ad}\nolimits(e^{tu})w-\mathop{\rm Ad}\nolimits(e^{tu})w\cdot u=[u,\mathop{\rm Ad}\nolimits(e^{tu})w]$; in particular at $t=0$ this is $[u,w]$. Moreover note that $\mathop{\rm Ad}\nolimits\mathcal{K}(\mathfrak{{p}})\subset\mathfrak{{p}}$ (i.e., the decomposition of $\mathfrak{{g}}$ is reductive). Hence if $v(0)=v_{\mathfrak{{k}}}(0)+v_{\mathfrak{{p}}}(0)$ and using $\mathop{\rm Ad}\nolimits(1)(u)=u$, we find that the curve

[TABLE]

satisfies the Euler-Arnold equation (i.e. $v_{\mathfrak{{k}}}(t)=v_{\mathfrak{{k}}}(0)$ and $v_{\mathfrak{{p}}}(t)=\mathop{\rm Ad}\nolimits e^{-2v_{\mathfrak{{k}}}(0)t}(v_{\mathfrak{{p}}}(0))$).

The reconstruction equation then implies [10]

[TABLE]

Note we can also write

[TABLE]

6.8 Symmetries of Hamiltonian

For $v\in\mathfrak{{g}}$, the Nother’s current $J$ associated to the lifted right action of $\mathcal{K}$ on $T\mathcal{G}$ is $J\circ\mathbb{F}\mathcal{L}(v)=\langle v,\cdot\rangle_{1}\in\mathfrak{{k}}^{*}$ (see cororally 4.2.11 [1]). Since $H$ is invariant under the action, by theorem 4.2.2 we know that $J$ is conserved along the Hamiltonian flow. In particular if $v(0)\in\mathfrak{{p}}=\mathfrak{{k}}^{\perp}$, i.e., $J\circ\mathbb{F}\mathcal{L}(v_{0})=0$, then $J\circ\mathbb{F}\mathcal{L}(v(t))\in\mathfrak{{p}}$ for all time.

6.9 Symplectic Integrators and Shadow Hamiltonian

That the flow of a Hamiltonian vector field preserve the symplectic structure is proposition 3.3.4 [1].

In order to simulate Hamilton’s equations we need to build flows that approximate the Hamiltonian flow $\dot{\gamma}=\{\gamma,H\}=\hat{H}(\gamma)$. Given an initial point $\gamma(0)\in\mathcal{G}\times\mathfrak{{g}}$, the solution to Hamilton’s equations satisfy \gamma(t)=\Phi^{\hat{H}}_{\gamma(0)}(t)\equiv\Phi^{\hat{H}}(t)\big{(}\gamma(0)\big{)}. We say an approximate flow $\Psi$ is said to be a symplectic integrator if it preserves the symplectic structure $\Psi^{*}\omega_{\mathcal{L}}=\omega_{\mathcal{L}}$ (technically we only need integrators that preserve the volume $\omega_{\mathcal{L}}^{n}$). Since Hamiltonian vector fields preserve the symplectic structure, a simple way to build integrators of $H$ consists of alternating between the flows of simpler Hamiltonians. For example the leapfrog integrator is defined as $\Psi_{\gamma(0)}(h)\equiv\Phi^{\hat{V}}(h/2)\circ\Phi^{\hat{T}}(h)\circ\Phi^{\hat{V}}_{\gamma(0)}(h/2)$. This integrator has the important property to also be reversible.

6.10 Lie Algebra

The tangent space $T_{1}\mathcal{G}$ may be identified with the vector space of left-invariant vector fields $\mathfrak{{g}}_{L}$, where we associate to each $\xi\in T_{1}\mathcal{G}$ the left-invariant vector field $X:g\mapsto\partial_{1}L_{g}\xi$, where left-invariance means $(L_{g})_{*}X=X$ and $(L_{g})_{*}$ is the push-forward. The lie algebra $\mathfrak{{g}}$ is $T_{1}\mathcal{G}$ equipped with the Lie bracket induced by the commutator of vector fields and the isomorphism $T_{1}\mathcal{G}\cong\mathfrak{{g}}_{L}$.

6.11 Horizontal Space

Recall that a principal bundle connection on $\mathcal{G}\to\mathcal{G}/\mathcal{K}$ is simply a smooth choice of horizontal spaces $H_{g}\subset T_{g}\mathcal{G}$, s.t. $T\mathcal{G}=ver_{g}\oplus H_{g}$ and \partial_{g}R_{k}\big{(}H_{g}\mathcal{G}\big{)}=H_{gk}\mathcal{G}. In particular we set $\mathfrak{{p}}\equiv H_{1}$.

If $v_{1},v_{2}\in\mathfrak{{p}}$ and $\partial_{1}\pi(v_{1})=\partial_{1}\pi(v_{2})$, then $\partial_{1}\pi(v_{1}-v_{2})=0$ so $v_{1}-v_{2}\in\mathfrak{{k}}$ and thus $v_{1}-v_{2}=0$, since $\mathfrak{{k}}$ and $\mathfrak{{p}}$ are complementary subspaces. Hence $v_{1}=v_{2}$.

6.12 Force Integrators

Using (3), and the fact $\{r,h\}=\omega_{\mathcal{L}}(\hat{r},\hat{h})=\hat{h}[r]$, we find

[TABLE]

where $c_{r}^{lk}\equiv g^{ik}g^{lj}c_{rji}$. In particular

[TABLE]

and \{V,\{V,T\}\}=g^{jk}\xi_{k}\big{(}v^{r}e_{r}(V)\big{)}e_{j}(V)=g^{jk}e_{k}(V)e_{j}(V). From (3) we find

[TABLE]

where we use Leibniz Rule.

To see that $-2g^{jk}g^{ls}e_{l}(V)e_{j}e_{s}(V)\xi_{k}\in\mathfrak{{p}}$, we first observe that $e_{j}e_{s}(V)=[e_{j},e_{s}](V)-e_{s}e_{j}(V)$. Moreover we can restrict to $s\in\mathfrak{{p}}$ since otherwise $e_{s}(V)=0$ (note here the basis fields $(e_{i})$ are generated from an orthogonal decomposition $\mathfrak{{g}}=\mathfrak{{k}}\oplus\mathfrak{{p}}$). Then $e_{j}(V)=0$ for all $j\in\mathfrak{{k}}$, and using the properties of the stucture constant we have $[e_{j},e_{s}]=c^{r}_{js}e_{r}$. We have $r\in\mathfrak{{p}}$ (again we have $e_{r}(V)=0$), and note

[TABLE]

Hence, since $\mathfrak{{p}}=\mathfrak{{k}}^{\perp}$, $k\in\mathfrak{{k}}$ implies $j\in\mathfrak{{k}}$, implies $c_{ajs}=-c_{sja}$ and since $g^{ls}g^{ra}e_{l}(V)e_{r}(V)$ is symmetric in $a,s$ the term sums to [math]. We have thus shown that $\widehat{\{V,\{V,T\}\}}\in\mathfrak{{p}}$.