Geometric Neural Ordinary Differential Equations: From Manifolds to Lie Groups

1. Introduction

Ordinary differential equations (ODEs) are ubiquitious in the engineering sciences, from modeling and control of simple physical systems like pendulums and mass-spring dampers, or more complicated robotic arms and drones, to the description of high-dimensional spatial discretizations of distributed systems, such as fluid flows, chemical reactions or quantum oscillators. Neural ordinary differential equations (neural ODEs) [1,2] are ODEs parameterized by neural networks. Given a state x, and parameters $\theta$ representing weights and biases of a neural network, a neural ODE reads:

$$\dot{x}=f_{\theta }(x,t),x\left(0\right)=x_0.$$

(1)

First introduced by [1] as the continuum limit of recurrent neural networks, the number of applications of neural ordinary differential equations quickly exploded beyond simple classification tasks: learning highly nonlinear dynamics of multi-physical systems from sparse data [3,4,5], optimal control of nonlinear systems [6], medical imaging [7] and real-time handling of irregular time-series[8], to name but a few. Discontinuous state-transitions and dynamics [9,10], time-dependent parameters [11], augmented neural ODEs [12] and physics preserving formulations [13,14] present further extensions that increase the expressivity of neural ODEs.

However, these methods are typically phrased for states $x{\in }^n$. For many physical systems of interest, such as robot arms, humanoid robots and drones, the state lives on differentiable manifolds and Lie groups [15,16]. More generally, the manifold hypothesis in machine learning raises the expectation that many high-dimensional data-sets evolve on intrinsically lower-dimensional, albeit more complicated manifolds [17]. Neural ODEs on manifolds [18,19] presented significant steps to address this gap, presenting first optimization methods for neural ODEs on manifolds. Yet, the general tools and approaches available on ${\mathbb{R}}^n$, such as including running costs, augmented states, time-dependent parameters, control-inputs or discontinuous state-transitions, are rarely addressed in a manifold context. Similar issues persist in a Lie group context, were Neural ODEs on Lie groups [20,21] were formalized.

Our goal is to extend more of the methods for neural ODEs on ${\mathbb{R}}^n$ to arbitrary manifolds, and in particular Lie groups. To this end, we present a systematic approach to the design of neural ODEs on manifolds and Lie groups. Specifically, our contributions are:

• Systematic derivation of neural ODEs on manifolds and Lie groups, highlighting differences and equivalence of various approaches - for an overview, see also Table 1;
• Summarizing the state of the art on manifold and Lie group neural ODEs, by formalizing the notion of extrinsic and intrinsic neural ODEs;
• A tutorial-like introduction, to assist the reader in implementing various neural ODE methods on manifolds and Lie groups, presenting coordinate expressions alongside geometric notation.

The remainder of the article is organized as follows. A brief state-of-the-art on neural ODEs concludes this introduction. Section 2 provides background on differentiable manifolds, Lie groups, and the coordinate-free adjoint method. Section 3 describes neural ODEs on manifolds, and derives parameter updates via the adjoint method for various common architectures and cost-functions, including time-dependent parameters, augmented neural ODEs, running costs, and intermediate cost-terms. Section 4 describes neural ODEs on matrix Lie groups, explaining merits of treating Lie groups separately from general differentiable manifolds. Both Section 3 and Section 4 also classify methods into extrinsic and intrinsic approaches. We conclude with a discussion in Section 5, highlighting advantages, disadvantages, challenges and promises of the presented material. The Appendix includes background on Hamiltonian systems, which appear when transforming the adjoint method into a form that is unique to Lie groups.

1.1. Literature review

For a general introduction to neural ODEs, see [25]. Neural ODEs on ${\mathbb{R}}^n$ with fixed parameters were first introduced by [1], and parameter optimization via the adjoint method allowed for intermittent and final cost terms on each trajectory. The generalized adjoint method [2] also allows for running cost terms. Memory-efficient checkpointing is introduced in [26] to address stability issues of adjoint methods. Augmented neural ODEs [12] introduced augmented state-spaces to allow neural ODEs to express arbitrary diffeomorphisms. Time-varying parameters were introduced by [11], with similar benefits to augmented neural ODEs. Neural ODEs with discrete transitions were formulated in [9,10], with [9] also learning event-triggered transitions common in engineering application. Neural controlled differential equations (CDEs) were introduced in [27] for handling irregular time-series, and parameter updates reapply the adjoint method [1]. Neural stochastic differential equations (SDEs) were introduced in [28], relying on a stochastic variant of the adjoint method for the parameter update. The previously mentioned literature phrases dynamics of neural ODEs on ${\mathbb{R}}^n$.

Neural ODEs on manifolds were first introduced by [22], including an adjoint method on manifolds for final cost terms, and application to continuous normalizing flows on Riemannian manifolds, but embedding manifolds into ${\mathbb{R}}^n$. Neural ODEs on Riemannian manifolds are expressed in local exponential charts in [18], avoiding embedding into ${\mathbb{R}}^n$, and considering final cost terms in the optimization. Charts for unknown, nontrivial latent manifolds, and dynamics in local charts, are learned from high-dimensional data in [29], including also discretized solutions to partial differential equations. Parameterized equivariant Neural ODEs on manifolds are constructed in [23], also commenting on state augmentation to express arbitrary (equivariant) flows on manifolds.

Neural ODEs on Lie groups were first introduced in [30] on the Lie group $SE\left(3\right)$, to learn the port-Hamiltonian dynamics of a drone from experiment, expressing group elements on an embedding ${\mathbb{R}}^{12}$, and the approach was formalized to port-Hamiltonian systems on arbitrary matrix Lie groups in [20], embedding $m\times m$ matrices in ${\mathbb{R}}^{m^2}$.

Neural ODEs on $SE\left(3\right)$ were phrased in local exponential charts in [24] to optimize a controller for a rigid body, using a chart-based adjoint method in local expential charts. An alternative, Lie algebra based adjoint method on general Lie groups was introduced in [21], foregoing Lie group specific numerical issues of applying the adjoint method in local charts.

1.2. Notation

For a complete introduction to differential geometry see e.g, [31], and for Lie group theory see [32].

Calligraphic letters $\mathcal{M},\mathcal{N},\dots$ denote smooth manifolds. For conceptual clarity, the reader may think of these manifolds as embedded in a high dimensional ${\mathcal{R}}^N$, e.g., $\mathcal{M}\subset {\mathcal{R}}^N$. The set $C^{\infty }(\mathcal{M},\mathcal{N})$ contains smooth functions between $\mathcal{M}$ and $\mathcal{N}$, and we define $C^{\infty }\left(\mathcal{M}\right):=C^{\infty }(\mathcal{M},\mathcal{R})$:.

The tangent space at $x\in \mathcal{M}$ is $T_x\mathcal{M}$ and the cotangent space is $T_x^{*}\mathcal{M}$. The tangent bundle of $\mathcal{M}$ is $T\mathcal{M}$, and the cotangent bundle of $\mathcal{M}$ is $T^{*}\mathcal{M}$. Then $\left(\mathcal{M}\right)$ denotes the set of vector fields over $\mathcal{M}$, and ${\mathsf{\Omega }}^k\left(\mathcal{M}\right)$ denotes the set of k forms, where ${\mathsf{\Omega }}^1\left(\mathcal{M}\right)$ are co-vector fields, and ${\mathsf{\Omega }}^0\left(\mathcal{M}\right)=C^{\infty }\left(\mathcal{M}\right)$ are smooth functions $V:\mathcal{M}\to \mathcal{R}$. The exterior derivative is denoted as $:{\mathsf{\Omega }}^k\left(\mathcal{M}\right)\to {\mathsf{\Omega }}^{k+1}\left(\mathcal{M}\right)$. For functions $V\in C^{\infty }(\mathcal{M}\times \mathcal{N},\mathbb{R})$, with $x\in \mathcal{M}$, $y\in \mathcal{N}$, we denote by $\left({‌}_{x}V\right)\left(y\right)\in T_x^{*}\mathcal{M}$ the partial differential at $x\in \mathcal{M}$. Curves $x:\mathcal{R}\to \mathcal{M}$ are denoted as $x\left(t\right)$, and their tangent vectors are denoted as $\dot{x}\in T_{x\left(t\right)}\mathcal{M}$.

A Lie group is denoted by G, its elements by $g,h$. The group identity is $e\in G$, and I denotes the identity matrix. The Lie algebra of G is , and its dual is ${‌}^{*}$. Letters $\tilde{A},\tilde{B}$ denote vectors in the Lie algebra, while letters $A,B$ denote vectors in ${\mathcal{R}}^n$.

In coordinate expressions lower indices are covariant and upper indices are contravariant components of tensors. For example, for a $(0,2)$-tensor M, the components $M_{ij}$ are covariant, and for non-degenerate M, the components of its inverse $M^{-1}$ are $M^{ij}$, which are contravariant. We use the Einstein summation convention $a_i{b}^i:={\sum }_i{a}_i{b}^i$, i.e., the product of variables with repeated lower and upper indices implies a sum.

Denoting W a topological space, D the Borel $\sigma$-algebra and $\mathbb{P}:D\to [0,1]$ a probability measure, then the tuple $(W,D,\mathbb{P})$ denotes a probability space. Given a vector space L and a random variable $C:X\to L$, then the expectation of C w.r.t. $\mathbb{P}$ is ${\mathbb{E}}_{w\sim \mathbb{P}}\left(C\right):={\int }_{W}C\left(w\right)\mathbb{P}\left(w\right)$.

2. Background

2.1. Smooth Manifolds

Given an n-dimensional manifold$\mathcal{M}$, with $U\subset \mathcal{M}$ an open set and $:U\to {\mathcal{R}}^n$ a homeomorphism, we call $(U,)$ a chart and we denote the coordinates of $x\in U$ as

$$(q^1,\dots {,}^{q^n}):=\left(x\right),x\in U\subset \mathcal{M}.$$

(2)

Smooth manifolds admit charts $\left(U_1{,}_1\right)$ and $\left(U_2{,}_2\right)$ with smooth transition maps${‌}_{21}{=}_2{\circ }_1^{-1}$ defined on the intersection $U_1\bigcap U_2$, and a collection $\mathcal{A}$ of charts $(U,)$ with smooth transition maps is called a (smooth) atlas. A vector field $f\in \left(\mathcal{M}\right)$ associates to any point $x\in \mathcal{M}$ a vector in $T_x\mathcal{M}$. It defines a dynamic system

$$\dot{x}=f\left(x\right);x\left(0\right)=x_0,$$

(3)

and we denote the solution of (3) by the flow operator

$${\mathsf{\Psi }}_f^t:\mathcal{M}\to \mathcal{M};{\mathsf{\Psi }}_f^t\left(x_0\right):=x\left(t\right).$$

(4)

For a real valued function $V\in C^{\infty }\left(\mathcal{M}\right)$, its differential is the covector field

$$V\in {\mathsf{\Omega }}^1\left(\mathcal{M}\right);V=\left[V\right]{{‌}^i}^i.$$

(5)

Given additionally a smooth manifold $\mathcal{N}$ and a smooth map $\phi :\mathcal{N}\to \mathcal{M}$, with $(U,)$ and $(\overline{U},\overline{‌})$ appropriate charts of $\mathcal{M}$ and $\mathcal{N}$, respectively. Then the pullback of V via $\phi$ is

$${\phi }^{*}V\in {\mathsf{\Omega }}^1\left(\mathcal{N}\right);{\phi }^{*}V:=(V\circ \phi )=\left[{\phi }^j\right]{\overline{‌}}^i\left[V\right]{‌}^j{\overline{‌}}^i.$$

(6)

With a Riemannian metric M (i.e., a symmetric, non-degenerate (0,2) tensor field) on $\mathcal{M}$, the gradient of V is a uniquely defined vector field $\nabla V\in \left(\mathcal{M}\right)$ given by

$$\nabla V:=M^{-1}V=M^{ij}\left[V\right]{‌}^j\left[\right]{‌}^i.$$

(7)

When $\mathcal{M}={\mathcal{R}}^n$, we assume that M is the Euclidean metric, and pick coordinates such that the components of the gradient and differential are the same. Finally, we define the Lie derivative of 1-forms, which differentiates $\omega \in {\mathsf{\Omega }}^1\left(\mathcal{M}\right)$ along a vector field $f\in \left(\mathcal{M}\right)$ and returns ${\mathcal{L}}_f\omega \in {\mathsf{\Omega }}^1\left(\mathcal{M}\right)$:

$${\mathcal{L}}_f\omega :={\displaystyle \frac{‌}{t}}{\left({{\mathsf{\Psi }}_f^t}^{*}\omega \right)}_{t=0}={\omega }_j{\left({\displaystyle \frac{\partial }{{\partial }^i}}{f}^j\right)}^i+\left({\displaystyle \frac{\partial }{{\partial }^j}}{\omega }_i\right)f^{j}i.$$

(8)

2.2. Lie groups

Lie groups are smooth manifolds with a compatible group structure. We consider real matrix lie groups $G\subseteq GL(m,\mathcal{R})$, i.e., subgroups of the general linear group

$$GL(m,\mathcal{R}):=\{g\in {\mathcal{R}}^{m\times m}|\det \left(g\right)\ne 0\},$$

(9)

For $g,h\in G$ the left and right translations by h are, respectively, the matrix multiplications

$$\begin{array}{c}\hfill L_h\left(g\right):=hg,\end{array}$$

(10)

$$\begin{array}{c}\hfill R_h\left(g\right):=gh.\end{array}$$

(11)

The Lie algebra of G is the vector space $\subseteq \mathfrak{gl}(m,\mathcal{R})$, with $\mathfrak{gl}(m,\mathcal{R})={\mathcal{R}}^{m\times m}$ the Lie algebra of $GL(m,\mathcal{R})$.

Define a basis $E:=\{{\tilde{E}}_1,\dots ,{\tilde{E}}_n\}$ with ${\tilde{E}}_i\in \subset {\mathcal{R}}^{m\times m}$, and Define the (invertible, linear) map $\Lambda :{\mathcal{R}}^n\to$ as1

$$\Lambda :{\mathcal{R}}^n\to ;(A^1,\dots ,A^n)\mapsto \sum _i{A}^i{\tilde{E}}_i.$$

(12)

The dual of is denoted ${‌}^{*}$, and given the map $\Lambda$ we call ${\Lambda }^{*}{:}^{*}\to {\mathcal{R}}^n$ its dual. For $\tilde{A},\tilde{B}\in$ the small adjoint ${\mathrm{ad}}_{\tilde{A}}\left(\tilde{B}\right)$ is a bilinear map, and the large Adjoint is the linear map ${\mathrm{Ad}}_g\left(\tilde{A}\right)$ is a linear map

$$\begin{array}{cc}\hfill \mathrm{ad}& :\times \to ;{\mathrm{ad}}_{\tilde{A}}\left(\tilde{B}\right)=\tilde{A}\tilde{B}-\tilde{B}\tilde{A},\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill \mathrm{Ad}& :G\times \to ;{\mathrm{Ad}}_g\left(\tilde{A}\right)=g\tilde{A}{g}^{-1}.\hfill \end{array}$$

(14)

In the remainder of the article, we exclusively use the adjoint representation ${\mathrm{ad}}_A:{\mathcal{R}}^n\to {\mathcal{R}}^n$, written without a tilde in the subscript A, and Adjoint representation ${\mathrm{Ad}}_g:{\mathcal{R}}^n\times {\mathcal{R}}^n$ which are obtained as

$$\begin{array}{cc}\hfill {\mathrm{ad}}_A& :={\Lambda }^{-1}\left({\mathrm{ad}}_{\Lambda \left(A\right)}\Lambda (·)\right),\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill {‌}_g& :={\Lambda }^{-1}\left({\mathrm{Ad}}_g\Lambda (·)\right)\hfill \end{array}$$

(16)

The exponential map $exp:\to G$ is a local diffeomorphism given by the matrix exponential [Chapter 3.7 [32]:

$$exp\left(\tilde{A}\right):=\sum _{n=0}^{\infty }{\displaystyle \frac{1}{n!}}{\tilde{A}}^n.$$

(17)

Its inverse $log:U_{log}\to$ is a given by the matrix logarithm, and it is well-defined on a subset $U_{log}\subset G$ [Chapter 2.3 [32]:

$$log\left(g\right)=\sum _{n=1}^{\infty }{(-1)}^{n+1}{\displaystyle \frac{{(g-I)}^n}{n}}.$$

(18)

Often, these infinite sums in (17) and (18) can further be reduced to a finite sums in m terms, by use of the Cayley-Hamilton Theorem [34]. A chart $\left(U_h{,}_h\right)$ on G that assigns zero coordinates to $h\in G$ can be defined using (18) and (12)

$$\begin{array}{cc}\hfill U_h& =\left\{hg\right|g\in U_{log}\},\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill {‌}_h& :U_h\to {\mathcal{R}}^n;g\mapsto {\Lambda }^{-1}log\left(h^{-1}g\right),\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill {‌}_h^{-1}& :{\mathcal{R}}^n\to G;q\mapsto hexp\left(\Lambda \left(q\right)\right).\hfill \end{array}$$

(21)

The chart $\left(U_h{,}_h\right)$ is called an exponential chart, and a collection $\mathcal{A}$ of exponential charts $\left(U_h{,}_h\right)$ that cover the manifold is called an exponential atlas.

The differential of a function $V\in C^{\infty }(G,\mathcal{R})$ is the co-vector field $V\in {\mathsf{\Omega }}^1\left(G\right)$ (see also Equation (5)). For any given $g\in G$ we further transform the co-vector $V\left(g\right)\in T_g^{*}G$ to a left-trivialized differential, which collects the components of the gradient expressed in ${‌}^{*}$:

$$\begin{array}{c}\hfill {‌}_g^{L}V:={\Lambda }^{*}{L}_g^{*}V\left(g\right)=\left[\right]‌V{\left(g\left(I+\Lambda (‌)\right)\right)}_{|=0}\in {\mathcal{R}}^n.\end{array}$$

(22)

For a derivation of this coordinate expression, see [Sec. 3 [21].

2.3. Gradient over a flow

We will be interested in computing the gradient of functions with respect to the initial state of a flow. The adjoint sensitivity equations are a set of differential equations that achieve this. In the following, we show a derivation of the adjoint sensitivity on manifolds [21]. Given a function $:\mathcal{M}\to \mathcal{R}$, a vector field $f\in \left(\mathcal{M}\right)$, the associated flow ${\mathsf{\Psi }}_f^t:\mathcal{M}\to \mathcal{M}$, and a final time $T\in \mathcal{R}$, the goal of the adjoint sensitvity method on manifolds is to compute the gradient

$$\left(C\circ {\mathsf{\Psi }}_f^T\right)\left(x_0\right).$$

In the adjoint method we define a co-state $\lambda \left(t\right)=(C\circ {\mathsf{\Psi }}^{T-t})\left(x\left(t\right)\right)\in T_{x\left(t\right)}^{*}\mathcal{M}$, which represents the differential of $C\left(x\left(T\right)\right)$ with respect to $x\left(t\right)$. The adjoint sensitivity method describes its dynamics, which are integrated backwards in time from the known final condition $\lambda \left(T\right)=C\left(x\left(T\right)\right)$, see also Figure 1. The adjoint sensitivity method is stated in Theorem 1.

Theorem 1 

(Adjoint sensitivity on manifolds). The gradient of a function $C\circ {\mathsf{\Psi }}_f^T$ is

$$\left(C\circ {\mathsf{\Psi }}_f^T\right)\left(x_0\right)=\lambda \left(0\right),$$

(23)

where $\lambda \left(t\right)\in T_{x\left(t\right)}^{*}\mathcal{M}$ is the co-state. In a local chart $(U,)$ of $\mathcal{M}$ with induced coordinates on $T^{*}U$, $x\left(t\right)$ and $\lambda \left(t\right)$ satisfy the dynamics

$$\begin{array}{cc}\hfill {\dot{‌}}^j& =f^j\left(\right),\left(0\right)=\left(x_0\right),\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill {\dot{\lambda }}_i& =-{\lambda }_j{\displaystyle \frac{\partial }{{\partial }^i}}{f}^j\left(\right),{\lambda }_i\left(T\right)=\left[C\right]{‌}^i\left(x\left(T\right)\right)\hfill \end{array}$$

(25)

Proof. 

Define the co-state $\lambda \left(t\right)\in T_{x\left(t\right)}^{*}\mathcal{M}$ as

$$\begin{array}{cc}\hfill \lambda \left(t\right)& :={\left({\mathsf{\Psi }}_f^{T-t}\right)}^{*}C\left(x\left(T\right)\right).\hfill \end{array}$$

(26)

Then Equation (23) is recovered by application of Equation (6):

$$\lambda \left(0\right)={\left({\mathsf{\Psi }}_f^T\right)}^{*}C\left(x\left(T\right)\right)=(C\circ {\mathsf{\Psi }}_f^T)\left(x_0\right),$$

(27)

A derivation of the dynamics governing $\lambda \left(t\right)$ constitutes the remainder of this proof. By definition of $\lambda \left(t\right)$ and the Lie derivative (8), we have that ${\mathcal{L}}_f\lambda \left(t\right)=0$:

$$\begin{array}{cc}\hfill {\mathcal{L}}_f\lambda \left(t\right)& ={\displaystyle \frac{d}{ds}}{\left({\left({\mathsf{\Psi }}_f^s\right)}^{*}\lambda (t+s)\right)}_{s=0}\hfill \\ \hfill & ={\displaystyle \frac{d}{ds}}\lambda \left(t\right)=0.\hfill \end{array}$$

(28)

If we further treat $\lambda$ as a 1-form $\lambda \in {\mathsf{\Omega }}^1\left(\mathcal{M}\right)$ (denoted as $\lambda$, by an abuse of notation), we obtain:

$$\begin{array}{cc}\hfill {\mathcal{L}}_f\lambda =& {\lambda }_j{\left({\displaystyle \frac{\partial }{{\partial }^i}}{f}^j\right)}^i+\left({\displaystyle \frac{\partial }{{\partial }^j}}{\lambda }_i\right)f^{j}i=0.\hfill \end{array}$$

The components satisfy the partial differential equation

$${\lambda }_j\left[\right]{‌}^i{f}^j+f^j\left[\right]{‌}^j{\lambda }_i=0.$$

(29)

Impose that $\lambda \left(t\right)=\lambda \left({\mathsf{\Psi }}_f^t\left(x_0\right)\right)$ (this defines the 1-form $\lambda$ along $x\left(t\right)$), then

$${\dot{\lambda }}_i=\left[{\lambda }_i\right]{‌}^j{\dot{‌}}^j=\left[{\lambda }_i\right]{‌}^j{f}^j.$$

(30)

Combining Equations (29) and (30) leads to Equation ():

$${\dot{\lambda }}_i=-{\lambda }_j{\displaystyle \frac{\partial }{{\partial }^i}}{f}^j.$$

(31)

Expanding the final condition $\lambda \left(T\right)=C\left(x\left(T\right)\right)$ in local coordinates (see Equation (5)) gives

$$\lambda \left(T\right)=\left[C\right]{‌}^i{\left(x\left(T\right)\right)}^i={\lambda }_i{\left(T\right)}^i\iff {\lambda }_i\left(T\right)=\left[C\right]{‌}^i\left(x\left(T\right)\right).$$

(32)

□

A fact that will become useful in Section 4, is that the equations (24) and () have a Hamiltonian form. Define the control Hamiltonian $H_c:T^{*}\mathcal{M}\to \mathbb{R}$ as

$$H_c(x,\lambda )=\lambda \left(f(x,t)\right)={\lambda }_i\left(f^i(,t)\right).$$

(33)

Then Equation (24) and Equation (), respectively, of Theorem 1 follow as the Hamiltonian equations on $T^{*}\mathcal{M}$:

$$\begin{array}{cc}\hfill {\dot{‌}}^j& ={\displaystyle \frac{\partial H_c}{\partial {\lambda }_j}}=f^j(,t),\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill {\dot{\lambda }}_i& =-{\displaystyle \frac{\partial H_c}{{\partial }^i}}=-{\lambda }_j{\displaystyle \frac{\partial }{{\partial }^i}}{f}^j(,t).\hfill \end{array}$$

(35)

For background on Hamilton’s equations, see also Appendix A.1.

3. Neural ODEs on Manifolds

A neural ODE on a manifold is an NN-parameterized vector field in $\left(\mathcal{M}\right)$ – or including time dependence, it is an NN-parameterized vector field in $(\mathcal{M}\times \mathcal{R})$, with t in the $\mathcal{R}$ slot and $\dot{t}=1$. Given parameters $\theta \in {\mathcal{R}}^{n_{\theta }}$, we denote this parameterized vector field as $f_{\theta }(x,t):=f(x,t,\theta )$. This results in the dynamic system

$$\dot{x}=f_{\theta }(x,t),x\left(0\right)=x_0,$$

(36)

The key idea of neural ODEs is to tackle various flow approximation tasks, by optimising the parameters with respect to a to-be-specified optimization problem. Denote a finite time horizon T and intermittent times $T_1,T_2,\dots <T$. Denote a general trajectory cost by

$$C_{f_{\theta }}^T(x_0,\theta )=F(\theta ,{\mathsf{\Psi }}_{f_{\theta }}^{T_0}\left(x_0\right),{\mathsf{\Psi }}_{f_{\theta }}^{T_1}\left(x_0\right),\dots ,{\mathsf{\Psi }}_{f_{\theta }}^T\left(x_0\right))+{\int }_0^{T}r({\mathsf{\Psi }}_{f_{\theta }}^s\left(x_0\right),s)s,$$

(37)

with intermittent and final cost term F, and running cost r. Indicating a probability space $(\mathcal{M},D,\mathbb{P})$, we define the total cost as

$$J\left(\theta \right):={\mathbb{E}}_{x_0\sim \mathbb{P}}{C}_{f_{\theta }}^T(x_0,\theta ).$$

(38)

The minimization problem takes the form

$$\begin{array}{c}\hfill \underset{\theta }{min}J\left(\theta \right).\end{array}$$

(39)

Note that (39) is not subject to any dynamic constraint - the flow already appears explicitly in the cost $C_{f_{\theta }}^T$.

Normally, the optimization problem is solved by means of a stochastic gradient descent optimization algorithm [35]. In this, a batch of N initial conditions $x_i$ is sampled from the probability distribution corresponding to the probability measure $\mathbb{P}$. Writing $C_i=C_{f_{\theta }}^T(x_i,\theta )$, the parameter gradient $\left[\right]\theta J\left(\theta \right)$ is approximated as

$${\displaystyle \frac{\partial }{\partial \theta }}J\left(\theta \right)={\mathbb{E}}_{x_0\sim \mathbb{P}}{\displaystyle \frac{\partial }{\partial \theta }}{C}_{f_{\theta }}^T\left(x_0\right)\approx {\displaystyle \frac{1}{N}}\sum _{i=0}^N{\displaystyle \frac{\partial }{\partial \theta }}{C}_i.$$

(40)

In this section, we show how to optimize the parameters $\theta$ for various choices of neural ODEs and cost functions, with (37) the most general case of a cost, and highlight similarities in the various derivations. In the following, the gradient ${\displaystyle \frac{\partial }{\partial \theta }}{C}_i$ is computed via the adjoint method on manifolds, for various scenarios. The advantage of the adjoint method over e.g., automatic differentiation of $C_i$ / back-propagation through an ODE solver, is that it has a constant memory efficiency with respect to the network depth T.

3.1. Constant parameters, running and final cost

Here we consider neural ODEs of the form (36), with constant parameters $\theta$, and cost functions of the form

$$C_{f_{\theta }}^T(x_0,\theta )=F({\mathsf{\Psi }}_{f_{\theta }}^T\left(x_0\right),\theta )+{\int }_0^{T}r({\mathsf{\Psi }}_{f_{\theta }}^s\left(x_0\right),\theta ,s)s,$$

(41)

with a final cost term F and a running cost term r. This generalizes [1,2] to manifolds. Compared to existing manifold methods for neural ODES [18,19], the running cost is new.

The parameter gradient’s components $\left[\right]\theta C_{f_{\theta }}^T\left((x_0,t_0),\theta \right)\in {\mathcal{R}}^{n_{\theta }}$ are then computed by Theorem 2 (see also [21]):

Theorem 2 

(Generalized Adjoint Method on Manifolds). Given the dynamics (36) and the cost (41), the parameter gradient’s components $\left[\right]\theta C_{f_{\theta }}^T\left((x_0,t_0),\theta \right)\in {\mathcal{R}}^{n_{\theta }}$ are computed by

$$\left[\right]\theta C_{f_{\theta }}^T\left((x_0,t_0),\theta \right)=\left(\left[F\right]\theta \right)(x\left(T\right),\theta )+{\int }_0^T\left[\right]\theta \left({\lambda }_j{f}_{\theta }^j\left(\left(s\right)\right)+r(\left(s\right),\theta ,s)\right)s.$$

(42)

where the state $x\left(s\right)\in \mathcal{M}$ and co-state $\lambda \left(s\right)\in T_{x\left(s\right)}^{*}\mathcal{M}$ satisfy, in a local chart $(U,)$ with $\left(t\right)=\left(x\right(t\left)\right)$, $\lambda \left(t\right)={\lambda }_i{\left(t\right)}^i$

$$\begin{array}{cc}\hfill {\dot{‌}}^j& =f_{\theta }^j(,t),\left(0\right)=\left(x_0\right),t\left(0\right)=t_0,\hfill \end{array}$$

(43)

$$\begin{array}{cc}\hfill {\dot{\lambda }}_i& =-{\lambda }_j{\displaystyle \frac{\partial }{{\partial }^i}}{f}_{\theta }^j(,t)-{\displaystyle \frac{\partial r}{{\partial }^i}},{\lambda }_i\left(T\right)=\left[F\right]{‌}^i\left(x\left(T\right),\theta \right).\hfill \end{array}$$

(44)

Proof. 

Define the augmented state space as ${\mathcal{M}}^{\prime }=\mathcal{M}\times {\mathcal{R}}^{n_{\theta }}\times \mathcal{R}\times \mathcal{R}$, to include the original state $x\in \mathcal{M}$, parameters $\theta \in {\mathcal{R}}^{n_{\theta }}$, accumulated running cost $L\in \mathcal{R}$ and time $t\in \mathcal{R}$ in the augmented state $x^{\prime }:=(x,\theta ,L,t)\in {\mathcal{M}}^{\prime }$. In addition, define the augmented dynamics $f_{\mathrm{aug}}\in \left({\mathcal{M}}^{\prime }\right)$ as

$${\dot{x}}^{\prime }=f_{\mathrm{aug}}\left(x^{\prime }\right)=\left(\begin{array}{c}{f}_{\theta }(x,t)\\ 0\\ r(x,\theta ,t)\\ 1\end{array}\right),x^{\prime }\left(0\right)=x_0^{\prime }:=\left(\begin{array}{c}{x}_0\\ \theta \\ 0\\ t_0\end{array}\right).$$

(45)

This is an autonomous system with final state $x^{\prime }\left(T\right)=(x\left(T\right),\theta ,{\int }_0^{T}r(x,\theta ,s)s,T)$. Next, define the cost $C_{\mathrm{aug}}:{\mathcal{M}}^{\prime }\to \mathcal{R}$ on the augmented space:

$$C_{\mathrm{aug}}\left(x^{\prime }\right)=F(x,\theta )+L.$$

(46)

Then Equation (41) can be rewritten as the evaluation of a terminal cost $C_{\mathrm{aug}}\left(x^{\prime }\left(T\right)\right)$:

$$C_{f_{\theta }}^T\left(x_0\right)=(C_{\mathrm{aug}}\circ {\mathsf{\Psi }}_{f_{\mathrm{aug}}}^T)\left(x_0^{\prime }\right).$$

(47)

By Theorem 1, the gradient $\left(C_{\mathrm{aug}}\circ {\mathsf{\Psi }}_{f_{\mathrm{aug}}}^T\right)$ is given by

$$(C_{\mathrm{aug}}\circ {\mathsf{\Psi }}_{f_{\mathrm{aug}}}^T)\left(x_0^{\prime }\right)=\lambda \left(0\right),$$

(48)

and by Equation (), the componentents of $\lambda \left(s\right)$ satisfy

$$\begin{array}{c}\hfill {\dot{\lambda }}_i=-{\lambda }_j{\displaystyle \frac{\partial }{\partial {{‌}^i}^{\prime }}}{f}_{\mathrm{aug}}^j,{\lambda }_i\left(T\right)=\left[\right]{{‌}^i}^{\prime }{C}_{\mathrm{aug}}\left(x^{\prime }\left(T\right)\right)\end{array}$$

(49)

Split the co-state into ${\lambda }_{,}{\lambda }_{\theta },{\lambda }_L,{\lambda }_t$, then their components’ dynamics are:

$$\begin{array}{cc}\hfill {\dot{\lambda }}_{,i}& =-\left[\right]{‌}^i\left({\lambda }_{,j}{f}_{\theta }^j(,t)+{\lambda }_{L}r(,\theta ,t)\right),{\lambda }_{(}T)=\left[F\right]‌(x\left(T\right),\theta ),\hfill \end{array}$$

(50)

$$\begin{array}{cc}\hfill {\dot{\lambda }}_{\theta ,i}& =-\left[\right]{\theta }^i\left({\lambda }_{,j}{f}_{\theta }^j(,t)+{\lambda }_{L}r(,\theta ,t)\right),{\lambda }_{\theta }\left(T\right)=\left[F\right]\theta (x\left(T\right),\theta )\hfill \end{array}$$

(51)

$$\begin{array}{cc}\hfill {\dot{\lambda }}_L& =0,{\lambda }_L\left(T\right)=\left[\right]L{C}_{\mathrm{aug}}(x\left(T\right),\theta )=1\hfill \end{array}$$

(52)

$$\begin{array}{cc}\hfill {\dot{\lambda }}_t& =-\left[\right]t\left({\lambda }_{,j}{f}_{\theta }^j(,t)+{\lambda }_{L}r(,\theta ,t)\right),{\lambda }_t\left(T\right)=\left[\right]t{C}_{\mathrm{aug}}(x\left(T\right),\theta )=0.\hfill \end{array}$$

(53)

The component ${\lambda }_L=1$ is constant, so equation (50) coincides with (). Integrating () from $s=0$ to $s=T$ recovers Equation (42). ${\lambda }_t$ does not appear in any of the other equations, such that Equation () may be ignored. □

In summary, the above proof dependeds on identifying a suitable augmented manifold ${\mathcal{M}}^{\prime }$, with the goal that augmented dynamics $f_{\mathrm{aug}}\in \left({\mathcal{M}}^{\prime }\right)$ are autonomous, and that the cost-function $C_{\mathrm{aug}}:{\mathcal{M}}^{\prime }\to \mathcal{R}$ on the augmented manifold rephrases the cost (41) as a final cost $C_{\mathrm{aug}}\left(x\left(T\right)\right)$, which allows to apply Theorem 1, and express any gradients of the original cost. In later sections (Sec. Section 3.2), this process will be the main technical tool for generalizations of Theorem (2). The next sections describe common special cases of (36) and Theorem 2.

3.1.1. Vanilla neural ODEs and extrinsic neural ODEs on manifolds

The case of neural ODEs on ${\mathcal{R}}^n$ (e.g., [1,2]) is obtained by setting $\mathcal{M}={\mathcal{R}}^n$. In this case, one global chart can be used to represent all quantities.

This overlaps with extrinsic neural ODEs on manifolds (described, for instance, in [22]), which optimize the neural ODE on an embedding space ${\mathcal{R}}^N$. We denote this embedding as $\iota :\mathcal{M}\to {\mathcal{R}}^N$, let $x\in \mathcal{M}$ and $y\in {\mathcal{R}}^N$. Optimizing the neural ODE on ${\mathcal{R}}^N$ requires extending the dynamics $f_{\theta }(x,t)\in \left(\mathcal{M}\right)$ to a vector field $f_{\theta }^{\uparrow }(y,t)\in \left({\mathcal{R}}^N\right)$, such that

$${\iota }_{*}{f}_{\theta }(x,t)=f_{\theta }^{\uparrow }(\iota \left(x\right),t)$$

(54)

The dynamics $f_{\theta }^{\uparrow }(y,t)$ are then used in Theorem 2, and also the costate lives in $T^{*}{\mathcal{R}}^N$.

As shown in [22], the resulting parameter gradients are equivalent to those resulting from an application in local charts, as long as it can be guaranteed that the integral curves of $f^{\uparrow }(y,t)\in X\left({\mathcal{R}}^N\right)$ remain within $\iota \left(\mathcal{M}\right)\subset {\mathcal{R}}^N$, i.e., are geometrically preserving.

A strong upside of an extrinsic formulation is that existing neural ODE packages (e.g., [36]) can be applied directly. Possible downsides of extrinsic neural ODEs are that finding $f^{\uparrow }(y,t)\in X\left({\mathcal{R}}^N\right)$ may not be immediate, a geometrically preserving integration has to be guaranteed separately, and that N is larger than the intrinsic dimension $n=dim\mathcal{M}$, leading to computational overhead.

3.1.2. Intrinsic neural ODEs on manifolds

The intrinsic case of neural ODEs on manifolds [18] is described by integrating the dynamics in local charts. Given a chart-transition from a chart $\left(U_1{,}_1\right)$ to a chart $\left(U_2{,}_2\right)$ chart transitions of the state and costate components (${‌}_1^i,{\lambda }_{1,i}$ and ${‌}_2^i,{\lambda }_{2,i}$, respectively) are given by

$$\begin{array}{cc}\hfill {‌}_2^i& {=}_2^i{\circ }_1^{-1}{(}_1),\hfill \end{array}$$

(55)

$$\begin{array}{cc}\hfill {\lambda }_{i,2}& =A_i^j{\lambda }_{j,1}.\hfill \end{array}$$

(56)

with $A_i^j={\partial }_i\left({‌}_1^j{\circ }_2^{-1}\right)$. The advantage of intrinsic neural ODEs on manifolds over extrinsic neural ODEs on manifolds is that the dimension of the resulting equations is as low as possible, for the given manifold. A disadvantage lies in having to determine charts, and chart-switching procedures. In available state-of-the-art packages for neural ODEs, these are phrased as discontinuous dynamics with state transitions ${‌}_1^j{\circ }_2^{-1}:{\mathcal{R}}^n\to {\mathcal{R}}^n$. For details on chart-switching methods see [18,21,29].

3.2. Extensions

The proof of Theorem 2 depended on identifying a suitable augmented manifold ${\mathcal{M}}^{\prime }$, autonomous augmented dynamics $f_{\mathrm{aug}}\in \left({\mathcal{M}}^{\prime }\right)$ and augmented cost-function $C_{\mathrm{aug}}:{\mathcal{M}}^{\prime }\to \mathcal{R}$ that rephrases the cost (41) as a final cost $C_{\mathrm{aug}}\left(x\left(T\right)\right)$, to apply Theorem 1. This approach generalizes to various other scenarios, including different cost-terms, augmented neural ODEs on manifolds and time-dependent parameters, presented in the following.

3.2.1. Nonlinear and intermittent cost terms

We consider here the case of neural ODEs on manifolds of the form (36) with cost (37). For the final and intermittent cost term $F_{\theta }:\mathcal{M}\times \mathcal{M}\times \dots \times \mathcal{M}\to \mathcal{R}$ we denote by ${‌}_k{F}_{\theta }\in T_x^{*}\mathcal{M}$ the differential w.r.t. the k-th slot, and denote $\theta$ as a subscript to avoid confusion. The components of ${‌}_{k}F$ will be denoted $\left[F\right]{‌}^{k}i$. In this case, the parameter gradient is determined by repeated application of Theorem 2:

Theorem 3 

(Generalized Adjoint Method on Manifolds). Given the dynamics (36) and the cost (37), the parameter gradient’s components $\left[\right]\theta C_{f_{\theta }}^T\left((x_0,t_0),\theta \right)\in {\mathcal{R}}^{n_{\theta }}$ are computed by

$$\begin{array}{cc}\hfill \left[\right]\theta C_{f_{\theta }}^T\left((x_0,t_0),\theta \right)=& \left(\left[F\right]\theta \right)(\theta ,x\left(T_1\right),x\left(T_2\right),\dots ,x\left(T\right))\hfill \\ \hfill & +{\int }_0^T\left[\right]\theta \left({\lambda }_j{f}_{\theta }^j\left(\left(s\right)\right)+r(\left(s\right),\theta ,s)\right)s.\hfill \end{array}$$

(57)

where the state $x\left(s\right)\in \mathcal{M}$ satisfies (43) and the co-state $\lambda \left(s\right)\in T_{x\left(s\right)}^{*}\mathcal{M}$ satisfies dynamics with discrete updates at times $T_1,\dots ,T_{N-1}$ given by

$$\begin{array}{cc}\hfill {\dot{\lambda }}_{,i}& =\left[\right]{‌}^i\left({\lambda }_{,j}{f}_{\theta }^j(,t)+r{(}_N,\theta ,t)\right);{\lambda }_{,i}\left(T\right)=\left[F_{\theta }\right]{‌}^{N}i(x\left(T_1\right),\dots ,x\left(T\right))\hfill \end{array}$$

(58)

$$\begin{array}{cc}\hfill {\lambda }_i\left(T_{k,-}\right)& ={\lambda }_i\left(T_{k,+}\right)+\left[F_{\theta }\right]{‌}^{k}i\left(x\left(T_1\right),\dots ,x\left(T\right)\right),\hfill \end{array}$$

(59)

with $T_{k,-}$ the instance after a discrete update at time $T_k$ (recall that co-state dynamics are integrated backwards, so $T_{k,-}<T_k<T_{k,+}$), and $T_{k,+}$ the instance before.

Proof. 

We introduce an augmented manifold ${\mathcal{M}}^{\prime }=\mathcal{M}\times \dots \times \mathcal{M}\times {\mathcal{R}}^{n_{\theta }}\times \mathcal{R}\times \mathcal{R}$, to include N copies of the original state $x\in \mathcal{M}$, parameters $\theta \in {\mathcal{R}}^{n_{\theta }}$, accumulated running cost $L\in \mathcal{R}$ and time $t\in \mathcal{R}$ in the augmented state $x^{\prime }:=(x_1,\dots ,x_N,\theta ,L,t)\in {\mathcal{M}}^{\prime }$. Let

$${\varrho }_{T_i}\left(t\right)=\left\{\begin{array}{c}1t\le T_i\hfill \\ 0t>T_i\hfill \end{array}\right.,$$

(60)

and define the augmented dynamics $f_{\mathrm{aug}}\in \left({\mathcal{M}}^{\prime }\right)$ as

$${\dot{x}}^{\prime }=f_{\mathrm{aug}}\left(x^{\prime }\right)=\left(\begin{array}{c}{\varrho }_{T_1}\left(t\right)f_{\theta }(x_1,t)\\ ⋮\\ {\varrho }_{T_{N-1}}\left(t\right)f_{\theta }(x_{N-1},t)\\ f_{\theta }(x_N,t)\\ 0\\ r(x_N,\theta ,t)\\ 1\end{array}\right),x^{\prime }\left(0\right)=x_0^{\prime }:=\left(\begin{array}{c}{x}_0\\ ⋮\\ x_0\\ x_0\\ \theta \\ 0\\ t_0\end{array}\right).$$

(61)

This is an autonomous system with final state

$$x^{\prime }\left(T\right)=(x\left(T_1\right),\dots ,x\left(T_{N-1}\right),x\left(T\right),\theta ,{\int }_0^{T}r(x,\theta ,s)s,T).$$

(62)

Next, define the cost $C_{\mathrm{aug}}:{\mathcal{M}}^{\prime }\to \mathcal{R}$ on the augmented space:

$$C_{\mathrm{aug}}\left(x^{\prime }\right)=F_{\theta }(x_1,\dots ,x_N)+L.$$

(63)

Then Equation (37) can be rewritten as the evaluation of a terminal cost $C_{\mathrm{aug}}\left(x^{\prime }\left(T\right)\right)$:

$$C_{f_{\theta }}^T\left(x_0\right)=(C_{\mathrm{aug}}\circ {\mathsf{\Psi }}_{f_{\mathrm{aug}}}^T)\left(x_0^{\prime }\right).$$

(64)

Apply Equation (), and split the co-state into ${\lambda }_{{‌}_1},\dots ,{\lambda }_{{‌}_N},{\lambda }_{\theta },{\lambda }_L,{\lambda }_t$, then their components’ dynamics are:

$$\begin{array}{cc}\hfill {\dot{\lambda }}_{{‌}_1,i}& =-\left[\right]{‌}^i\left({\lambda }_{{‌}_1,j}{\varrho }_{T_1}\left(t\right)f_{\theta }^j{(}_1,t)\right),{\lambda }_{{‌}_1}\left(T\right)=\left[F_{\theta }\right]{‌}^1(x\left(T_1\right),\dots ,x\left(T\right)),\hfill \\ \hfill ⋮\end{array}$$

(65)

$$\begin{array}{cc}\hfill {\dot{\lambda }}_{{‌}_N,i}& =-\left[\right]{‌}^{N}i\left({\lambda }_{{‌}_N,j}{f}_{\theta }^j{(}_N,t)+{\lambda }_{L}r{(}_N,\theta ,t)\right),{\lambda }_{{‌}_N}\left(T\right)=\left[F_{\theta }\right]{‌}^N(x\left(T_1\right),\dots ,x\left(T\right)),\hfill \end{array}$$

(66)

$$\begin{array}{cc}\hfill {\dot{\lambda }}_{\theta ,i}& =-\left[\right]{\theta }^i\left({\lambda }_{{‌}_1,j}{\varrho }_{T_1}\left(t\right)f_{\theta }^j{(}_1,t)+\dots +{\lambda }_{{‌}_N,j}{f}_{\theta }^j{(}_N,t)+{\lambda }_{L}r(,\theta ,t)\right),\hfill \\ \hfill {\lambda }_{\theta }\left(T\right)& =\left[F_{\theta }\right]\theta (x\left(T_1\right),\dots ,x\left(T\right)).\hfill \end{array}$$

(67)

We excluded the dynamics of ${\lambda }_t$, which does not appear in any of the other equations, and the constant ${\lambda }_L=1$. Finally, define the cumulative co-state

$${\lambda }_{‌}={\varrho }_{T_1}\left(t\right){\lambda }_{{‌}_1}+\dots +{\lambda }_{{‌}_N}.$$

(68)

Its dynamics at $t\in [0,T]T_1,·,T_{N-1}$ are given by the sum of (65) to (), letting ${=}_N$:

$$\begin{array}{cc}\hfill {\dot{\lambda }}_{,i}& ={\dot{\lambda }}_{{‌}_1,i}+\dots +{\dot{\lambda }}_{{‌}_N,i}\hfill \end{array}$$

(69)

$$\begin{array}{cc}\hfill & =\left[\right]{‌}^i\left({\lambda }_{,j}{f}_{\theta }^j(,t)+r{(}_N,\theta ,t)\right)\hfill \end{array}$$

(70)

$$\begin{array}{cc}\hfill {\lambda }_{‌}\left(T\right)& =\left[F_{\theta }\right]{‌}^N(x\left(T_1\right),\dots ,x\left(T\right)),\hfill \end{array}$$

(71)

with discrete jumps (58) accounting for the final conditions of ${\lambda }_{{‌}_1},\dots ,{\lambda }_{{‌}_N}$, and the dynamics () can be rewritten as

$${\dot{\lambda }}_{\theta ,i}=\left[\right]{\theta }^i\left({\lambda }_{,j}{f}_{\theta }^j(,t)+r(,\theta ,t)\right);{\lambda }_{\theta }\left(T\right)=\left[F_{\theta }\right]\theta (x\left(T_1\right),\dots ,x\left(T\right)).$$

(72)

Integrating this from $s=0$ to $s=T$ recovers Equation (57). □

Cost terms of this form are interesting for optimiziation of e.g., periodic orbits [37] or trajectories on manifolds, where conditions at multiple checkpoints ${\mathsf{\Psi }}_{f_{\theta }}^{T_i}\left(x_0\right)$ may appear in the cost.

3.2.2. Augmented neural ODEs on manifolds and time-dependent parameters

With state $x\in \mathcal{M}$, augmented state $\alpha \in \mathcal{N}$ (not to be confused with $x^{\prime }\in {\mathcal{M}}^{\prime }$), and parameterized ${\phi }_{\theta }:\mathcal{M}\to \mathcal{N}$, augmented neural ODEs on manifolds are neural ODEs on the manifold $\mathcal{M}\times \mathcal{N}$ of the form

$$\left(\begin{array}{c}\dot{x}\\ \dot{\alpha }\end{array}\right)=\left(\begin{array}{c}{f}_{\theta }(x,\alpha )\\ g_{\theta }(x,\alpha )\end{array}\right);\left(\begin{array}{c}x\left(0\right)\\ \alpha \left(0\right)\end{array}\right)=\left(\begin{array}{c}{x}_0\\ {\phi }_{\theta }\left(x_0\right)\end{array}\right).$$

(73)

Time t is not included explicitly in these dynamics, since it can be included in $\alpha$. This case also includes the scenario of time-dependent parameters $\overline{\theta }\left(t\right)$ as part of $\alpha$. As the trajectory cost, we take a final cost

$$C_{f_{\theta },g_{\theta }}^T(x_0,\theta )=F({\mathsf{\Psi }}_{f_{\theta },g_{\theta }}^T\left(x_0,{\phi }_{\theta }\left(x_0\right)\right),\theta ).$$

(74)

Theorem 4 

(Adjoint Method for Augmented Neural ODEs on Manifolds). Given the dynamics (73) and the cost (74), the parameter gradient’s components $\left[\right]\theta C_{f_{\theta }}^T\left((x_0,t_0),\theta \right)\in {\mathcal{R}}^{n_{\theta }}$ are computed by

$$\begin{array}{cc}\hfill \left[\right]\theta C_{f_{\theta },g_{\theta }}^T\left((x_0,\phi \left(x_0\right)),\theta \right)=& \left(\left[F\right]\theta \right)(x\left(T\right),\alpha \left(T\right),\theta )+\left[{\phi }^j\right]\theta {\lambda }_{\alpha ,j}\left(0\right)\hfill \\ \hfill & +{\int }_0^T\left[\right]\theta \left({\lambda }_{x,j}{f}_{\theta }^j\left(\left(s\right)\right)+{\lambda }_{\alpha ,j}{g}_{\theta }^j\left(\left(s\right)\right)\right)s.\hfill \end{array}$$

(75)

where the states $x\left(s\right)\in \mathcal{M},\alpha \left(s\right)\in \mathcal{N}$ satisfy (73) and co-states ${\lambda }_x\left(s\right)\in T_{x\left(s\right)}^{*}\mathcal{M},{\lambda }_{\alpha }\left(s\right)\in T_{\alpha \left(s\right)}^{*}\mathcal{N},$ satisfy, in a local chart $(U,)$ on $\mathcal{M}$ and $\overline{U},\overline{‌}$ on $\mathcal{N}$:

$$\begin{array}{cc}\hfill {\dot{\lambda }}_{x,i}& =-{\displaystyle \frac{\partial }{{\partial }^i}}\left({\lambda }_{x,j}{f}_{\theta }^j(,\overline{‌},t)+{\lambda }_{\alpha ,j}{g}_{\theta }^j(,\overline{‌},t)\right),{\lambda }_{x,i}\left(T\right)=\left[F\right]{‌}^i\left(x\left(T\right),\alpha \left(T\right),\theta \right),\hfill \end{array}$$

(76)

$$\begin{array}{cc}\hfill {\dot{\lambda }}_{\alpha ,i}& =-{\displaystyle \frac{\partial }{\partial {\overline{‌}}^i}}\left({\lambda }_{x,j}{f}_{\theta }^j(,\overline{‌},t)+{\lambda }_{\alpha ,j}{g}_{\theta }^j(,\overline{‌},t)\right),{\lambda }_{\alpha ,i}\left(T\right)=\left[F\right]{\overline{‌}}^i\left(x\left(T\right),\alpha \left(T\right),\theta \right).\hfill \end{array}$$

(77)

Proof. 

Define the augmented state space as ${\mathcal{M}}^{\prime }=\mathcal{M}\times \mathcal{N}\times {\mathcal{R}}^{n_{\theta }}$, to include the states $x\in \mathcal{M},\alpha \left(s\right)\in \mathcal{N}$ and parameters $\theta \in {\mathcal{R}}^{n_{\theta }}$ in the augmented state $x^{\prime }:=(x,\alpha ,\theta )\in {\mathcal{M}}^{\prime }$. In addition, define the augmented dynamics $f_{\mathrm{aug}}\in \left({\mathcal{M}}^{\prime }\right)$ as

$${\dot{x}}^{\prime }=f_{\mathrm{aug}}\left(x^{\prime }\right)=\left(\begin{array}{c}{f}_{\theta }(x,\alpha )\\ g_{\theta }(x,\alpha )\\ 0\end{array}\right),x^{\prime }\left(0\right)=x_0^{\prime }:=\left(\begin{array}{c}{x}_0\\ {\phi }_{\theta }\left(x_0\right)\\ \theta \end{array}\right).$$

(78)

This is an autonomous system with final state $x^{\prime }\left(T\right)=(x\left(T\right),\alpha \left(T\right),\theta )$. Next, define the cost $C_{\mathrm{aug}}:{\mathcal{M}}^{\prime }\to \mathcal{R}$ on the augmented space:

$$C_{\mathrm{aug}}\left(x^{\prime }\right)=F(x,\alpha ,\theta ).$$

(79)

Then Equation (41) can be rewritten as the evaluation of a terminal cost $C_{\mathrm{aug}}\left(x^{\prime }\left(T\right)\right)$. The gradient $\left(C_{\mathrm{aug}}\circ {\mathsf{\Psi }}_{f_{\mathrm{aug}}}^T\right)$ is given by an application of Equation (). Split the co-state into ${\lambda }_x,{\lambda }_{\alpha },{\lambda }_{\theta }$, then their components’ dynamics are:

$$\begin{array}{cc}\hfill {\dot{\lambda }}_{x,i}& =-{\displaystyle \frac{\partial }{{\partial }^i}}\left({\lambda }_{x,j}{f}_{\theta }^j(,\overline{‌},t)+{\lambda }_{\alpha ,j}{g}_{\theta }^j(,\overline{‌},t)\right),{\lambda }_x\left(T\right)=\left[F\right]‌(x\left(T\right),\alpha \left(T\right),\theta ),\hfill \end{array}$$

(80)

$$\begin{array}{cc}\hfill {\dot{\lambda }}_{\alpha ,i}& =-{\displaystyle \frac{\partial }{\partial {\overline{‌}}^i}}\left({\lambda }_{x,j}{f}_{\theta }^j(,\overline{‌},t)+{\lambda }_{\alpha ,j}{g}_{\theta }^j(,\overline{‌},t)\right),{\lambda }_{\alpha ,i}\left(T\right)=\left[F\right]{\overline{‌}}^i\left(x\left(T\right),\alpha \left(T\right),\theta \right).\hfill \end{array}$$

(81)

$$\begin{array}{cc}\hfill {\dot{\lambda }}_{\theta ,i}& =-\left[\right]{\theta }^i({\lambda }_{x,j}{f}_{\theta }^j(,\overline{‌},t)+{\lambda }_{\alpha ,j}{g}_{\theta }^j(,\overline{‌},t)),{\lambda }_{\theta }\left(T\right)=\left[F\right]\theta (x\left(T\right),\alpha \left(T\right),\theta )\hfill \end{array}$$

(82)

Since $\alpha \left(0\right)={\phi }_{\theta }\left(x_0\right)$ also depends on $\theta$, the total gradient of the cost w.r.t. $\theta$ is given by

$$\left[\right]{\theta }^i{C}_{f_{\theta }}^T\left((x_0,{\phi }_{\theta }\left(x_0\right)),\theta \right)={\lambda }_{\theta ,i}\left(0\right)+\left[{\phi }^j\right]{\theta }^i{\lambda }_{\alpha ,j}\left(0\right).$$

(83)

This recovers Equation (75). □

A further degenerate application of Theorem 4 is obtained by removing x. Then both dynamics $g_{\theta }\left(\alpha \right)$ and initial condition $\alpha \left(0\right)={\phi }_{\theta }\left(0\right)$ are parameterized by $\theta$, allowing joint optimization of parameters and initial condition. This is interesting for joint optimization and numerical continuation, e.g. [37].

4. Neural ODEs on Lie Groups

Just as a neural ODE on a manifold is an NN-parameterized vector field in $\left(\mathcal{M}\right)$ (or, including time, $(\mathcal{M}\times \mathcal{R})$), a neural ODE on a Lie group can be seen as a parameterized vector field in $\left(G\right)$ (or $(G\times \mathcal{R})$, respectively). Similarly to Equation (36), this results in a dynamic system

$$\dot{g}=f_{\theta }(g,t),g\left(0\right)=g_0,$$

(84)

Yet, Lie groups offers more structure than manifolds: the Lie algebra provides a canonical space to represent tangent vectors, and its dual ${‌}^{*}$ provides a canonical space to represent the co-state. Similarly, canonical (exponential) charts offer structure for integrating dynamic systems [38]. Frequently, dynamics on a Lie group induced dynamics on a manifold $\mathcal{M}$: by means of an action

$$\Phi :G\times \mathcal{M}\to \mathcal{M};(g,x)\mapsto \Phi (g,x),$$

(85)

evolutions $g\left(t\right)$ induce evolutions $x\left(t\right)=\Phi (g\left(t\right),x_0)$ on $\mathcal{M}$. This makes neural ODEs on Lie groups interesting in their own right.

In this section, we describe optimizing (39) for the cost

$$C_{f_{\theta }}^T(g_0,\theta )=F({\mathsf{\Psi }}_{f_{\theta }}^T\left(g_0\right),\theta )+{\int }_0^{T}r({\mathsf{\Psi }}_{f_{\theta }}^s\left(g_0\right),\theta ,s)s,$$

(86)

with a final cost term F and a running cost term r. We highlight the extrinsic approach, and two intrinsic approaches, where one of the latter is peculiar to Lie groups.

4.1. Extrinsic neural ODEs on Lie groups

The extrinsic formulation of neural ODEs on Lie groups was first introduced by [20], and applies ideas of [19] (see also Section 3.1.1). Given $G\subset GL(m,\mathcal{R})$, this formulation treats the dynamic system (84) as a dynamic system on ${\mathcal{R}}^{m^2}$. Denote $\mathrm{vec}:{\mathcal{R}}^{m\times m}\to {\mathcal{R}}^{m^2}$ an invertible map that stacks the components of an input matrix matrix into a component vector2 and let ${\mathrm{proj}}_G:{\mathcal{R}}^{m\times m}\to G$ be a projection onto $G\subset {\mathcal{R}}^{m\times m}$. Further denote $A_y={\mathrm{vec}}^{-1}\left(y\right)$ and $g_y={\mathrm{proj}}_G\left(A_y\right)$. A lift $f_{\theta }^{\uparrow }(y,t)$ can then be defined as

$$f_{\theta }^{\uparrow }(y,t)=\mathrm{vec}\left(A_y{g}_y^{-1}f(g_y,\theta ,t)\right).$$

(87)

As was the case for extrinsic neural ODEs on manifolds, the cost-gradient resulting from this optimization is well-defined and equivalent to any intrinsically defined procedure. However, dimension $m^2$ of the vectorization can be significantly larger than the intrinsic dimension of the Lie group.

4.2. Intrinsic neural ODEs on Lie groups

Theorem 2 directly applies to optimization of neural ODEs on Lie groups, given the local exponential charts (19), () on G. This does not make full use of the available structure on Lie groups. Frequently, dynamical systems are of a left-invariant form (88) or a right-invariant form ()

$$\begin{array}{cc}\hfill \dot{g}& =g\Lambda \left({\rho }_{\theta }^L(g,t)\right),\hfill \end{array}$$

(88)

$$\begin{array}{cc}\hfill \dot{g}& =\Lambda \left({\rho }_{\theta }^R(g,t)\right)g.\hfill \end{array}$$

(89)

Denoting $K\left(\right):T_{‌}{\mathcal{R}}^n\to {\mathcal{R}}^n$ the derivative of the exponential map (see [21] for details). Then the chart-representatives $f_{\theta }^i$ in a local exponential chart $\left(U_h{,}_h\right)$ are

$$\begin{array}{cc}\hfill f_{\theta }^{L,i}(,t)& ={\left(K^{-1}\right)}_j^i\left(\right){\rho }^{L,j}{(}_h^{-1}\left(\right)),\hfill \end{array}$$

(90)

$$\begin{array}{cc}\hfill f_{\theta }^{R,i}(,t)& ={\left(K^{-1}\right)}_j^i\left(\right){\mathrm{Ad}}_{{‌}_h^{-1}\left(\right)}{\rho }^{R,j}{(}_h^{-1}\left(\right)).\hfill \end{array}$$

(91)

Application of Theorem 2 then requires computing $\left[\right]{‌}^j{f}_{\theta }^{L,i}(,t)$ or $\left[\right]{‌}^j{f}_{\theta }^{R,i}(,t)$. But this leads to significant computational overhead due differentiation of the terms ${\left(K^{-1}\right)}_j^i\left(\right)$ (see [21]). Instead of applying Theorem 2, i.e., expressing dynamics in local charts, the dynamics can also be expressed at the Lie algebra . Theorem has a Hamiltonian form, which can be directly transformed into Hamiltonian equations on a Lie group (see also Appendix A.1). Applying this reasoning to Theorem 2, we arrive at the following form, which foregoes differentiating ${\left(K^{-1}\right)}_j^i\left(x\right)$:

Theorem 5 

(Left Generalized Adjoint Method on Matrix Lie Groups). Given are the dynamics (88) and the cost (86), or the dynamics () with ${\rho }_{\theta }^L(g,t)={\mathrm{Ad}}_{g^{-1}}{\rho }_{\theta }^R(g,t)$. Then the parameter gradient $\left[\right]\theta C_{f_{\theta }}^T\left(g_0\right)$ of the cost is given by the integral equation

$$\left[\right]\theta C_{f_{\theta }}^T\left(g_0\right)=\left[F\right]\theta (g\left(T\right),\theta )+{\int }_0^T\left[\right]\theta \left({\lambda }_g^{\top }{\rho }_{\theta }^L(g,s)+r(g,\theta ,s)\right)s,$$

(92)

where the state $g\left(t\right)\in G$ and co-state ${\lambda }_g\left(t\right)\in {\mathbb{R}}^n$ are the solutions of the system of equations

$$\begin{array}{cc}\hfill \dot{g}& =f_{\theta }(g,t),g\left(0\right)=g_0,\hfill \end{array}$$

(93)

$$\begin{array}{cc}\hfill {\dot{\lambda }}_g& ={-}_g^L\left({\lambda }_g^{\top }{\rho }_{\theta }^L(g,s)+r(g,\theta ,s)\right){+}_{{\rho }_{\theta }^L(g,t)}^{\top }{\lambda }_g,{\lambda }_g\left(T\right){=}_g^{L}F(g\left(T\right),\theta ).\hfill \end{array}$$

(94)

Proof. 

This is proven in two steps. First, define the time-and-parameter-dependent control-Hamiltonian $H_c:T^{*}\mathcal{M}\times {\mathcal{R}}^{n_{\theta }}\times \mathbb{R}\to \mathbb{R}$ as

$$H_c(x,\lambda ,\theta ,t)=\lambda \left(f_{\theta }(x,t)\right)+r(x,\theta ,t)={\lambda }_i\left(f_{\theta }^i(,t)\right)+r(,\theta ,t).$$

(95)

The equations for the state- and co-state dynamics (43) and (), respectively, of Theorem 2 follow as the Hamiltonian equations on $T^{*}\mathcal{M}$:

$$\begin{array}{cc}\hfill {\dot{‌}}^j& ={\displaystyle \frac{\partial H_c}{\partial {\lambda }_j}}=f_{\theta }^j(,t),\hfill \end{array}$$

(96)

$$\begin{array}{cc}\hfill {\dot{\lambda }}_i& =-{\displaystyle \frac{\partial H_c}{{\partial }^i}}=-{\lambda }_j{\displaystyle \frac{\partial }{{\partial }^i}}{f}_{\theta }^j(,t)-{\displaystyle \frac{\partial r}{{\partial }^i}}.\hfill \end{array}$$

(97)

And the integral equation (42) reads

$$\left[\right]\theta C_{f_{\theta }}^T\left((x_0,t_0),\theta \right)=\left[F\right]\theta (x\left(T\right),\theta )+{\int }_0^T\left[H_c\right]\theta t.$$

(98)

Second, rewrite the control Hamiltonian (95) on a Lie group G, i.e., $H_c:T^{*}G\times \times {\mathcal{R}}^{n_{\theta }}\times \mathbb{R}\to \mathbb{R}$. By substituting ${\lambda }_g\left(t\right)={\Lambda }^{*}{L}_g^{*}\lambda \left(t\right)$ (see also Equation (A6)), this induces ${\mathcal{H}}_c:G\times {\mathfrak{g}}^{*}\times {\mathcal{R}}^{n_{\theta }}\times \mathbb{R}\to \mathbb{R}$

$${\mathcal{H}}_c(g,{\lambda }_g,\theta ,t)={\lambda }_g^{\top }{\rho }_{\theta }^L(g,t)+r(g,\theta ,t).$$

(99)

Finally Hamilton’s equations (96), () are rewritten in their form on a matrix Lie group by means of (A7), (), which recovers equations (93) and ():

$$\begin{array}{cc}\hfill \dot{g}& =g\Lambda \left(\left[{\mathcal{H}}_c\right]{\lambda }_g\right),\hfill \end{array}$$

(100)

$$\begin{array}{cc}\hfill {\dot{\lambda }}_g& ={-}_g^L{\mathcal{H}}_c+{\mathrm{ad}}_{\left[{\mathcal{H}}_c\right]{\lambda }_g}^{\top }{\lambda }_g,\hfill \end{array}$$

(101)

To find the final condition for ${\lambda }_g$, use that ${\lambda }_g\left(t\right)={\Lambda }^{*}{L}_g^{*}\lambda \left(t\right)$:

$${\lambda }_g\left(T\right)={\Lambda }^{*}{L}_g^{*}\lambda \left(T\right)={\Lambda }^{*}{L}_g^{*}F(g\left(T\right),\theta )=d_g^{L}F(g\left(T\right),\theta ).$$

(102)

□

Similar equations also hold on abstract (non-matrix) Lie groups, see [21]. Compared to the extrinsic method of Section 4.1, Theorem 5 has the advantage that the dimension of the co-state ${\lambda }_g$ is as low as possible. Compared to the chart-based approach on Lie groups, Theorem 5 foregoes differentiating through the terms $K_j^i\left(\right)$, avoiding overhead. Compared to a chart-based approach on manifolds, also the choice of charts is canonical on Lie groups. Although the Lie-group approach foregoes many of the pitfalls of intrinsic neural ODEs on manifolds, implementation in existing neural ODE packages is currently cumbersome: the adjoint-sensitivity equations () have a non-standard form, requiring an adapted dynamics of the co-state $\lambda$, but these equations are rarely intended for modification, in existing packages. Packages for geometry-preserving integrators on Lie groups, such as [38] are also not readily available for arbitrary Lie groups.

4.3. Extensions

The proof of Theorem 5 relied on finding a control-Hamiltonian formulation for Theorem 2. This approach generalizes to methods in Section 3.2, which rely on the use of Theorem 1. That is because Theorem 1 itself has a Hamiltonian form ([19,21]).

5. Discussion

We discuss advantages and disadvantages of the main flavors of the presented formulations for manifold neural ODEs, expanding on the previous sections. We focus on extrinsic (embedding dynamics in ${\mathcal{R}}^N$) and intrinsic (integrating in local charts) formulations. Summarizing prior comments:

• the extrinsic formulation is readily implemented if the low-dimensional manifold $\mathcal{M}$ and an embedding into ${\mathcal{R}}^N$ are known. This comes at the possible cost of geometric inexactness, and a higher dimension of the co-state and sensitivity equations
• the co-state in the intrinsic formulation has a generally lower dimension, which reduces the dimension of the sensitivity equations. The chart-based formulation also guarantees geometrically exact integration of dynamics. This comes at the mild cost of having to define local charts and chart-transitions.

This dimensionality reduction is unlikely to have a high impact when the manifold $\mathcal{M}$ is known and low dimensional, e.g., for the sphere $\mathcal{M}=S^2$ or similar manifolds. However, when applying the manifold hypothesis to high-dimensional data, there might be non-trivial latent manifolds for which the embedding is not immediate, and where the latent manifold is of a much lower dimension than the embedding data-manifold. Then the intrinsic method becomes difficult to avoid. If geometric exactness of the integration is desired, local charts need to be defined also for the extrinsic approach, in which case the intrinsic approach may offer further advantages.

In order to derive neural ODEs on Lie groups, three approaches were possible: the extrinsic and intrinsic formulations on manifolds directly carry over to matrix Lie groups, embedding $G\subset GL(m,\mathcal{R})$ in ${\mathcal{R}}^{m^2}$ or using local exponential charts, respectively. A third option a novel intrinsic method for neural ODEs on matrix Lie groups, which made full use of the Lie group structure by phrasing dynamics on (as is more common on Lie groups), and the co-state on ${‌}^{*}$, avoiding difficulties of the chart-based formalism in differentiating extra terms.

Summarizing prior comments on advantages and disadvantages of these flavors:

• the extrinsic formulation on matrix Lie groups can come at much higher cost than that on manifolds, since the intrinsic dimension of G can be much lower than $m^2$, and a higher dimension of the co-state and sensitivity equations. Geometrically exact integration procedures are more readily available for matrix Lie groups, integrating $\dot{g}$ in local exponential charts.
• the chart-based formulation on matrix Lie groups struggles when are not naturally phrased in local charts. This is common, dynamics are often more naturally phrased on . This was alleviated by an algebra-based formulation on matrix Lie groups. Both are intrinsic approaches, that feature a co-state dynamics that are as low as posssible. However, the algebra-based approach still misses readily available software implementation.

The authors believe that the algebra-based formulation is more convenient, in principle, and consider software implementations of the algebra-based approach as possible future work.

In summary, we presented a unified, geometric approach to extend various methods for neural ODEs on ${\mathcal{R}}^N$ to neural ODEs on manifolds and Lie groups. Optimization of neural ODEs on manifolds was based on the adjoint method on manifolds. Given a novel cost-function C and neural ODE architecture f, the strategy to present the results in a unified fashion was to identify a suitable augmented manifold ${\mathcal{M}}_{\mathrm{aug}}$, augmented dynamics $f_{\mathrm{aug}}\in \left({\mathcal{M}}_{\mathrm{aug}}\right)$, and cost $C_{\mathrm{aug}}:{\mathcal{M}}_{\mathrm{aug}}\to \mathcal{R}$ such that the original cost-function can be rephrased as $C=C_{\mathrm{aug}}\circ {\mathsf{\Psi }}_{f_{\mathrm{aug}}}^T$. To further derive optimization of intrinsic neural ODEs on Lie groups was based on finding a Hamiltonian formulation of the adjoint method on manifolds, and to subsequently transformed them into the Hamiltonian equations on a matrix Lie group.

Appendix A.1. Hamiltonian dynamics on Lie groups

We briefly review Hamiltonian systems on manifolds and matrix Lie groups (see also [21]).

Given a manifold $\mathcal{Q}$ with coordinate maps ${‌}^i:\mathcal{Q}\to \mathbb{R}$ and $p_i$ in the basis ${‌}^i$ on $T_q^{*}\mathcal{Q}$, we define the symplectic form $\omega \in {\mathsf{\Omega }}^2\left(T^{*}\mathcal{M}\right)$ as

$$\omega =p_i{\wedge }^i.$$

(A1)

Let $Y\in \left(T^{*}\mathcal{Q}\right)$, then a Hamiltonian $H\in C^{\infty }(T^{*}\mathcal{Q},\mathbb{R})$ implicitly defines a unique vector field $X_H\in \left(T^{*}\mathcal{Q}\right)$ by

$$H\left(Y\right)=\omega (X_H,Y).$$

(A2)

In coordinates, $X_H$ has the components

$$\begin{array}{cc}\hfill {\dot{‌}}^i& =\left[H\right]p_i,\hfill \end{array}$$

(A3)

$$\begin{array}{cc}\hfill {\dot{p}}_i& =-\left[H\right]{‌}^i.\hfill \end{array}$$

(A4)

On a Lie group G, the group structure allows the identification $T^{*}G\equiv G\times {\mathfrak{g}}^{*}\equiv G\times {\mathcal{R}}^n$. E.g., using the pull back ${L_g}^{*}:T_g^{*}G\to {\mathfrak{g}}^{*}$ of left-translation map $L_g:G\to G$, and ${\Lambda }^{*}:\to {\mathcal{R}}^n$, to define $P_g\in {\mathcal{R}}^n$ as

$$P_g={\Lambda }^{*}{L}_g^{*}P.$$

(A5)

Then the left Hamiltonian ${\mathcal{H}}^L:G\times {\mathfrak{g}}^{*}\to \mathbb{R}$ is defined in terms of $H:T^{*}G\to \mathfrak{g}$ as

$${\mathcal{H}}^L(g,P_g)=H\left(g,P)\right).$$

(A6)

For a matrix Lie group the left Hamiltonian equations read:

$$\begin{array}{cc}\hfill \dot{g}& =g\Lambda \left(\left[{\mathcal{H}}^L\right]P\right),\hfill \end{array}$$

(A7)

$$\begin{array}{cc}\hfill \dot{P}& ={-}_g^L{\mathcal{H}}^L+{\mathrm{ad}}_{\left[{\mathcal{H}}^L\right]P}^{\top }P,\hfill \end{array}$$

(A8)

with $\Lambda :{\mathbb{R}}^n\to \mathfrak{g}$ as in (12) and ${‌}_g^L\mathcal{H}\in {\mathbb{R}}^n$ as in (22).