  We derive a family of high-order, structure-preserving approximations of the
Riemannian exponential map on several matrix manifolds, including the group of
unitary matrices, the Grassmannian manifold, and the Stiefel manifold. Our
derivation is inspired by the observation that if $\Omega$ is a skew-Hermitian
matrix and $t$ is a sufficiently small scalar, then there exists a polynomial
of degree $n$ in $t\Omega$ (namely, a Bessel polynomial) whose polar
decomposition delivers an approximation of $e^{t\Omega}$ with error
$O(t^{2n+1})$. We prove this fact and then leverage it to derive high-order
approximations of the Riemannian exponential map on the Grassmannian and
Stiefel manifolds. Along the way, we derive related results concerning the
supercloseness of the geometric and arithmetic means of unitary matrices.

Numerical Analysis