RKMK Numerical Integration on Product Manifolds

Lie Group × Euclidean Factors

Abstract

Many flight-dynamics and GNC state vectors naturally live on product manifolds, e.g., attitude on \(\SO(3)\) (or unit quaternions) together with translational states in \(\R^n\). Standard Runge–Kutta methods assume a vector space and can break manifold constraints (e.g. orthonormality, unit norm). Runge–Kutta–Munthe-Kaas (RKMK) methods lift the integration to the Lie algebra and update the group by an exponential (or a retraction), achieving high order while preserving the manifold structure.

These notes explain RKMK in a product-manifold setting, emphasizing the typical case \(\SO(3)\times \R^n\), and how to consistently couple the Lie-group and Euclidean factors.

Why product manifolds show up in GNC?

A common state layout is

\[x = (R,\, p,\, v,\, \omega,\, b,\, \dots) \in \SO(3)\times \R^3\times \R^3\times \R^3\times \R^m \times \cdots.\]

The configuration (attitude) is not a vector: \(R\in \SO(3)\) satisfies \(R^\top R = I\) and \(\det(R)=1\). If we apply a vanilla Runge–Kutta scheme directly to a matrix parameterization (or integrate a quaternion and re-normalize), we risk:

drift from the manifold (loss of orthonormality / unit norm),
loss of accuracy order due to ad-hoc projection/normalization,
inconsistency between continuous-time model and discrete-time propagation used by an EKF (Extended Kalman Filter).

RKMK addresses this by:

performing Runge–Kutta stages in the Lie algebra (a vector space),
mapping algebra increments back to the group via \(\exp\) (or a retraction),
integrating Euclidean factors normally (addition).

Product manifolds and tangent spaces

Definition

A product manifold is \(M = M_1 \times \cdots \times M_k\) where each \(M_i\) is a manifold. A point is \(x=(x_1,\dots,x_k)\) with \(x_i\in M_i\).

Tangent space of a product

The tangent space factorizes:

\[T_{(x_1,\dots,x_k)}M \;\cong\; T_{x_1}M_1 \times \cdots \times T_{x_k}M_k.\]

So a vector field \(\dot x = F(x)\) on \(M\) splits into components

\[\dot x_i = F_i(x_1,\dots,x_k), \qquad i=1,\dots,k,\]

but the components can still be coupled through the shared dependence on all factors.

When factors are Lie groups

If \(G\) is a Lie group, then each tangent space \(T_g G\) can be identified with the Lie algebra \(\Lie{g} = T_e G\) via left translation:

\[T_g G \ni \dot g \quad\leftrightarrow\quad \xi = g^{-1}\dot g \in \Lie{g}.\]

For matrix Lie groups (like \(\SO(3)\)), \(g^{-1} = g^\top\) and \(\xi\) is a matrix in the algebra.

For a direct-product Lie group \(G = G_1 \times G_2\), the Lie algebra is the direct sum

\[\Lie{g} = \Lie{g}_1 \oplus \Lie{g}_2,\]

and \(\exp_G\) acts component-wise:

\[\exp_G(\xi_1,\xi_2) = (\exp_{G_1}(\xi_1),\exp_{G_2}(\xi_2)).\]

In practice, most GNC states look like a Lie group factor \(G\) times Euclidean \(\R^n\) (which is itself a Lie group under addition).

Dynamics on \(G\times \R^n\)

Let \(M = G \times \R^n\), with state \(x=(g,r)\). A broad class of coupled dynamics is:

\[\dot g = g\, \Omega(g,r), \qquad \dot r = f(g,r),\]

where \(\Omega(g,r)\in \Lie{g}\) is the left-trivialized “angular” rate (for \(G=\SO(3)\) this corresponds to body angular velocity), and \(f:\,G\times \R^n\to \R^n\) is the Euclidean part.

Key idea: \(\dot g\) lives on the group; \(\dot r\) lives in a vector space. RKMK integrates them together without leaving \(G\).

Retractions and the exponential map

Exponential on the group

For a Lie group \(G\), the exponential map

\[\exp:\ \Lie{g} \to G\]

maps an algebra element (a tangent vector at identity) to a group element.

Given a base point \(g_n\in G\), an increment \(\eta\in \Lie{g}\) updates the state as

\[g_{n+1} = g_n \exp(\eta).\]

This is a retraction based at \(g_n\).

Euclidean factor as a trivial Lie group

On \(\R^n\), the exponential is just addition:

\[r_{n+1} = r_n + \Delta r.\]

This matches the product-manifold view: we retract on \(G\) by \(\exp\), and on \(\R^n\) by \(+\).

From classical RK to RKMK

Classical Runge–Kutta (RK) methods are defined for ODEs in a vector space:

\[\dot y = \phi(y),\quad y\in \R^m.\]

RKMK adapts an RK tableau \((a_{ij}, b_i, c_i)\) to Lie-group evolution by:

representing stage states as \(g_n \exp(\eta_i)\) for algebra stage variables \(\eta_i\),
evaluating the group vector field in left-trivialized form \(\Omega(\cdot)\),
correcting for the nonlinear relationship between \(\eta\) and \(g\) using \(\dexp^{-1}\).

The \(\dexp\) operator

Define \(\eta(t)\in \Lie{g}\) by writing the group trajectory as

\[g(t) = g_n \exp(\eta(t)).\]

Differentiate:

\[\dot g(t) = g_n \frac{d}{dt}\exp(\eta(t)).\]

Left-trivialize at \(g(t)\):

\[g(t)^{-1}\dot g(t) = \Omega(g(t), r(t)).\]

A standard identity gives

\[g(t)^{-1}\dot g(t) = \dexp_{\eta(t)}(\dot \eta(t)),\]

so

\[\dot \eta(t) = \dexp_{\eta(t)}^{-1}\!\big(\Omega(g_n\exp(\eta(t)), r(t))\big).\]

This is the core: it converts the group ODE into an algebra ODE (a vector space), where we can apply RK.

RKMK on \(G\times \R^n\): coupled stage equations

Let step size be \(h\) and tableau \((a_{ij},b_i,c_i)\) have \(s\) stages. We propagate from \((g_n,r_n)\) to \((g_{n+1},r_{n+1})\).

Stage reconstruction on the product manifold

We maintain:

\[\eta_i \in \Lie{g},\qquad \rho_i \in \R^n,\]

and define stage states

\[g_i = g_n \exp(\eta_i),\qquad r_i = r_n + \rho_i.\]

Stage derivatives

Compute:

\[K_i = \dexp_{\eta_i}^{-1}\!\left(\Omega(g_i,r_i)\right)\in \Lie{g}, \qquad k_i = f(g_i,r_i)\in \R^n.\]

Stage accumulation

For an explicit R-K tableau, we set

\[\eta_i = h \sum_{j=1}^{i-1} a_{ij} K_j,\qquad \rho_i = h \sum_{j=1}^{i-1} a_{ij} k_j.\]

For implicit tableaus, these become coupled non-linear equations within the stages.

Step update

Finally,

\[\eta_{n+1} = h\sum_{i=1}^s b_i K_i,\qquad r_{n+1} = r_n + h\sum_{i=1}^s b_i k_i,\qquad g_{n+1} = g_n \exp(\eta_{n+1}).\]

Interpretation on a product manifold:

The group factor is advanced by a single exponential at the end (plus exponentials inside stages).
The Euclidean factor is advanced by the same R-K weights, but with ordinary addition.
Coupling enters only through the evaluations \(\Omega(g_i,r_i)\) and \(f(g_i,r_i)\).

Concrete special case: \(\SO(3)\times \R^n\)

SO(3) algebra and the hat operator

Identify \(\Lie{so}(3)\cong \R^3\) via the hat map:

\[\begin{split}\phi \in \R^3 \quad\mapsto\quad [\phi]_\times \in \Lie{so}(3), \qquad [\phi]_\times = \begin{bmatrix} 0 & -\phi_3 & \phi_2\\ \phi_3 & 0 & -\phi_1\\ -\phi_2 & \phi_1 & 0 \end{bmatrix}.\end{split}\]

Write \(\eta \equiv \phi\) (vector form) and \(\exp(\eta) \equiv \Exp([\phi]_\times)\).

Exponential map on SO(3)

Let \(\theta=\norm{\phi}\), \(A=[\phi]_\times\). Then Rodrigues’ formula:

\[\Exp(A) = I + \frac{\sin\theta}{\theta}A + \frac{1-\cos\theta}{\theta^2}A^2,\]

with the usual smooth limits as \(\theta\to 0\) (use series to avoid division by small \(\theta\)).

Left Jacobian and \(\dexp^{-1}\) on SO(3)

For \(\SO(3)\), \(\dexp\) corresponds to the left Jacobian \(J(\phi)\):

\[\dexp_{\phi}(u) \equiv J(\phi)\,u,\qquad u\in\R^3.\]

A standard closed form is

\[J(\phi) = I - \frac{1-\cos\theta}{\theta^2}[\phi]_\times + \frac{\theta-\sin\theta}{\theta^3}[\phi]_\times^2.\]

Hence

\[\dexp^{-1}_{\phi}(v) \equiv J(\phi)^{-1} v.\]

A commonly used closed form for \(J(\phi)^{-1}\) (for \(\theta\not\approx 0\)) is

\[J(\phi)^{-1} = I + \frac{1}{2}[\phi]_\times + \left(\frac{1}{\theta^2} - \frac{1+\cos\theta}{2\theta\sin\theta}\right)[\phi]_\times^2.\]

For small \(\theta\), use the series expansion

\[J(\phi)^{-1} \approx I + \frac{1}{2}[\phi]_\times + \frac{1}{12}[\phi]_\times^2 + \mathcal{O}(\theta^3).\]

This is often the most AD-friendly route: branch on \(\theta^2\) and use series near zero.

Typical coupled state example

Let \((R, z)\in \SO(3)\times \R^n\) and define

\[\dot R = R[\omega(R,z)]_\times,\qquad \dot z = f(R,z).\]

Then in RKMK you compute stages:

\[K_i = J(\phi_i)^{-1}\,\omega(R_i,z_i),\]

with \(R_i = R_n \Exp([\phi_i]_\times)\), \(z_i = z_n + \rho_i\), and the same R-K accumulation rules.

Note \(\R^n\) is also a Lie group

\((\R^n,+)\) is a Lie group with Lie algebra \(\R^n\), exponential \(\exp(v)=v\) (identity), and group action by addition. So in principle you can treat the whole state as a Lie group element in a direct product \(G\times \R^n\).

In practice, you still keep the Euclidean part in vector form because:

its exponential and \(\dexp\) are trivial,
you avoid unnecessary machinery,
it keeps Jacobians/simple linear algebra straightforward.

RKMK on product manifolds is exactly this: nontrivial Lie factors get \(\exp\) and \(\dexp^{-1}\); Euclidean factors reduce to standard R-K.

Explicit example: RK4-RKMK on \(G\times \R^n\)

Using classical RK4 tableau, stages are explicit. A compact description:

Given \((g_n,r_n)\) and step \(h\).

\[\eta_1 = 0,\quad \rho_1=0, \qquad (g_1,r_1) = (g_n,r_n),\]

\[K_1 = \dexp_{\eta_1}^{-1}\!\left(\Omega(g_1,r_1)\right),\quad k_1 = f(g_1,r_1),\]

\[\eta_2 = \tfrac{h}{2}K_1,\quad \rho_2=\tfrac{h}{2}k_1,\qquad (g_2,r_2) = (g_n\exp(\eta_2),\, r_n+\rho_2),\]

\[K_2 = \dexp_{\eta_2}^{-1}\!\left(\Omega(g_2,r_2)\right),\quad k_2 = f(g_2,r_2),\]

\[\eta_3 = \tfrac{h}{2}K_2,\quad \rho_3=\tfrac{h}{2}k_2,\qquad (g_3,r_3) = (g_n\exp(\eta_3),\, r_n+\rho_3),\]

\[K_3 = \dexp_{\eta_3}^{-1}\!\left(\Omega(g_3,r_3)\right),\quad k_3 = f(g_3,r_3),\]

\[\eta_4 = hK_3,\quad \rho_4=h k_3,\qquad (g_4,r_4) = (g_n\exp(\eta_4),\, r_n+\rho_4),\]

\[K_4 = \dexp_{\eta_4}^{-1}\!\left(\Omega(g_4,r_4)\right),\quad k_4 = f(g_4,r_4).\]

Then update:

\[\eta_{n+1} = \frac{h}{6}(K_1+2K_2+2K_3+K_4),\qquad r_{n+1} = r_n + \frac{h}{6}(k_1+2k_2+2k_3+k_4),\]

\[g_{n+1} = g_n \exp(\eta_{n+1}).\]

Implicit RKMK on product manifolds (stiff dynamics)

When dynamics are stiff (e.g. strong aerodynamic damping, tight actuator models, flexible modes, or coupled constraints), implicit R-K methods (Radau IIA, Gauss–Legendre) are attractive.

In implicit RKMK, the stage equations

\[\eta_i = h\sum_{j=1}^s a_{ij}K_j,\qquad \rho_i = h\sum_{j=1}^s a_{ij}k_j\]

are solved simultaneously with

\[K_i = \dexp_{\eta_i}^{-1}\!\left(\Omega(g_i,r_i)\right),\qquad k_i = f(g_i,r_i).\]

This forms a nonlinear system in \(\{\eta_i,\rho_i\}_{i=1}^s\).

Newton structure on a product manifold

Your unknown vector is a stacked Euclidean vector consisting of:

Lie algebra coordinates for each \(\eta_i\) (e.g. \(\R^3\) per \(\SO(3)\) factor),
Euclidean increments \(\rho_i\in\R^n\).

The nonlinearities enter through \(\exp(\eta_i)\) and \(\dexp^{-1}_{\eta_i}\) (plus whatever nonlinearities are in \(\Omega\) and \(f\)). This is exactly where algorithmic differentiation (e.g. CppAD) is helpful.

Quaternion normalization vs. Lie-group integration

If you represent attitude by a quaternion \(q\in S^3\) and integrate \(\dot q = \tfrac12 q\otimes \omega\) with RK, you often re-normalize: \(q \leftarrow q/\norm{q}\).

This is a projection back to \(S^3\). It can work well, but:

the normalization is an extra nonlinear operation not accounted for in R-K order conditions,
it can reduce the global order (especially for higher-order schemes) or complicate consistent linearization,
the EKF propagation Jacobians become sensitive to the chosen “projection” step.

RKMK instead updates by a group exponential (or a retraction) and stays on the manifold by construction:

\[R_{n+1} = R_n\Exp([\phi]_\times),\quad\text{or}\quad q_{n+1} = q_n \otimes \exp_{S^3}(\tfrac12\phi).\]

That tends to be cleaner for both accuracy and filter consistency.

Summary for \(G\times\R^n\) implementations

Write group dynamics in left-trivialized form: \(\dot g = g\,\Omega(g,r)\).
Choose a base R-K tableau (explicit or implicit).
In each stage, reconstruct \((g_i,r_i)\) via \((g_n\exp(\eta_i),\, r_n+\rho_i)\).
Evaluate \(\Omega(g_i,r_i)\) and \(f(g_i,r_i)\).
Convert the group rate to the algebraic rate with \(\dexp^{-1}_{\eta_i}\).
Accumulate stages using R-K coefficients.
Update \(g\) with a final \(\exp(\cdot)\) and \(r\) with a final weighted sum.
For \(\SO(3)\), implement stable \(\Exp\), \(J(\phi)\), \(J(\phi)^{-1}\) with small-angle series.

Appendix

Useful SO(3) identities

For \(\phi\in\R^3\) and \(A=[\phi]_\times\):

\[A^2 = \phi\phi^\top - \theta^2 I,\qquad \theta=\norm{\phi}.\]

These help implement Rodrigues formula and Jacobians efficiently.