# The Fisher Geometry and Geodesics of the Multivariate Normals, without
Differential Geometry

Brodie A. J. Lawson111Centre for Data Science, Queensland University of
Technology, Australia 222School of Mathematical Sciences, Queensland
University of Technology, Australia 333ARC Centre of Excellence for Plant
Success in Nature and Agriculture, Australia Kevin Burrage††footnotemark:
††footnotemark:  444Department of Computer Science, University of Oxford,
United Kingdom (Visiting Professorship) Kerrie Mengersen††footnotemark:
††footnotemark:  Rodrigo Weber dos Santos555Graduate Program in Computational
Modeling, Universidade Federal de Juiz de Fora, Brazil

###### Abstract

Choosing the Fisher information as the metric tensor for a Riemannian manifold
provides a powerful yet fundamental way to understand statistical distribution
families. Distances along this manifold become a compelling measure of
statistical distance, and paths of shorter distance improve sampling
techniques that leverage a sequence of distributions in their operation.
Unfortunately, even for a distribution as generally tractable as the
multivariate normal distribution, this information geometry proves unwieldy
enough that closed-form solutions for shortest-distance paths or their lengths
remain unavailable outside of limited special cases. In this review we present
for general statisticians the most practical aspects of the Fisher geometry
for this fundamental distribution family. Rather than a differential geometric
treatment, we use an intuitive understanding of the covariance-induced
curvature of this manifold to unify the special cases with known closed-form
solution and review approximate solutions for the general case. We also use
the multivariate normal information geometry to better understand the paths or
distances commonly used in statistics (annealing, Wasserstein). Given the
unavailability of a general solution, we also discuss the methods used for
numerically obtaining geodesics in the space of multivariate normals,
identifying remaining challenges and suggesting methodological improvements.

Keywords: Information geometry, Information theory, Statistical distance

##  1 Introduction

The concept of statistical distance, some measure of discrepancy between two
distributions, is ubiquitous in statistics. Classical tests such as the
Cramér–von Mises criterion or Anderson–Darling test use a distance between
(cumulative) density functions to identify whether given data is expected to
have come from a given distribution [1]. The expectation maximisation
algorithm, and variational approaches to Bayesian inference, operate by
minimising an information theoretic discrepancy measure — the Kullback–Leibler
(KL) divergence [2]. The Wasserstein distance (also known as the
Kantrovich–Rubenstein metric or earth-mover’s distance) describes in a sense
the amount of work required to physically reshape one distribution into
another, a concept that brings with it both appealing mathematics [3], and
demonstrated potential in statistics and machine learning [4, 5, 6]. These are
only a few examples of the many statistical discrepancy measures that see use
(for example see also [7, 8]).

Spoiled for choice in such a manner, it is often desirable to select a
discrepancy that is mathematically specified by a distance metric, such that
we may understand our measure of statistical distance through analogs to the
familiar Euclidean distance (satisfying key properties such as symmetry,
positivity, and the triangle inequality). If we also choose this metric such
that our notion of distance is tied to statistical information, then lengths
of paths through distribution space become measures of the information
required to traverse that path, and shortest-distance paths become “paths of
least statistical resistance”. This is of practical concern, as can be seen
within sampling approaches using distribution sequences such as sequential
Monte Carlo [9] or parallel tempering [10], as well as normalising constant
(ratio) estimators such as annealed importance sampling [11, 12].

The Fisher information matrix serves as just such a (Riemannian) metric. Use
of this metric began with Rao in 1945, who wile presenting the well-known
Cramér–Rao bound in statistical estimation also discussed the idea of members
of a distribution family sitting on what is now known as a “statistical
manifold” [13]. The geometry of this Riemannian manifold is warped by the
information already present in the distribution, in a fashion similar to the
way Riemannian manifolds in astrophysics describe the gravitational warping of
space-time [14, Ch. 10]. Especially after their modernisation by Shun’ichi
[15], these ideas of information geometry have often concentrated on the
theoretical properties of the manifolds associated with statistical models or
distribution families [16], but have also seen use in generating new
statistical algorithms [17], explaining convergence behaviour [18] and making
learning algorithms parameterisation-invariant [19]. As a means of quantifying
statistical distance, information geometry has a natural application to
hypothesis testing [20], and is also intrinsically linked to the question of
identifiability of parameters in models of complex systems [21, 22].

For a given family of distributions, information geometry provides an
information-aware measure of distance termed the Fisher–Rao distance (or
sometimes the Rao distance) that behaves like a distance, unlike say the KL
divergence [8]. However, owing to the way Riemannian geometry defines distance
in an infinitesimal fashion, finding the Fisher–Rao distance between two
points on the statistical manifold requires first finding the shortest-
distance path along the manifold between those two points. Known as geodesics,
such paths generalise straight lines to a non-Euclidean geometry, but are
challenging to generate. This, together with the general requirement of
familiarity with differential geometry, means that information geometry
remains somewhat unwieldy.

In this work, we provide an introduction to, and practical understanding of,
information geometry for general statisticians. We do so by avoiding as many
differential geometric concepts and as much of the corresponding notation as
possible, and by concentrating on a fundamental distribution family, the
multivariate normals (MVNs). What may at first appear to be a simplistic
choice reveals a rich topic area of ongoing research that we review, with
applications for example to segmentation of diffusion tensor imaging data [23]
and wireless signal processing [24, 25]. We provide an intuitive way of
understanding this statistical manifold’s geometry, and the form taken by its
geodesic paths. This approach allows us to unify the existing cases for which
a closed-form expression for the geodesic is known, and informs approximations
for use when a closed-form solution is not available. We also discuss how
other commonly-used paths through distribution space (annealing/geometric,
Wasserstein-optimal transport) compare to these information-optimal paths in
the space of MVNs, providing a deeper understanding of their function and
expected performance.

Given the need to numerically find the geodesic in the majority of cases where
a closed-form solution is unavailable, we also thoroughly discuss the shooting
methods used for this purpose. As part of our review, we link recent
developments in approximately-geodesic curves to shooting method performance,
by demonstrating the potential of the former as an initialisation. We conclude
with a discussion of how the concepts we discuss extend beyond the MVN family,
and the unsolved problem of a closed form for the geodesic and Fisher–Rao
distance between arbitrary MVNs.

##  2 The Fisher Geometry of the Normal Distribution

###  2.1 The Riemannian Metric

Consider a family of distributions parameterised by $\boldsymbol{\theta}$. For
the MVN family we consider here, the typical choice of parameters is the mean
vector, $\boldsymbol{\mu}$, and covariance matrix, $\boldsymbol{\Sigma}$. As
such, the parameter space for $d$-variate MVNs could be taken as
$\boldsymbol{\theta}\in\mathbb{R}^{d+d(d+1)/2}$, with the vector
$\boldsymbol{\theta}$ housing the $d$ components of $\boldsymbol{\mu}$ and the
$d(d+1)/2$ unique elements of the symmetric $\boldsymbol{\Sigma}$. Owing to
the requirement that $\boldsymbol{\Sigma}$ be a valid (positive definite)
covariance matrix, the parameter space for the family of $d$-variate normals
in fact occupies some subspace of $\mathbb{R}^{d+d(d+1)/2}$.

Asked for the distance between two points in the parameter space for
$\boldsymbol{\theta}$, we could say calculate the Euclidean norm of the
difference in their associated vectors,
$\|\boldsymbol{\theta}_{2}-\boldsymbol{\theta}_{1}\|_{2}$. This would
completely ignore the statistical information represented within
$\boldsymbol{\theta}$, through its constituent components $\boldsymbol{\mu}$
and $\boldsymbol{\Sigma}$ and what they imply about the MVN-distributed random
variables. Information geometry rectifies this, by treating the parameter
space instead as a curved manifold that explicitly acknowledges this
statistical context.

As stated, our goal here is to avoid a differential geometric treatment of the
statistical manifold associated with MVNs, which has already been
comprehensively provided [20, 26, 27]. However, in describing the (Riemannian)
geometry of this manifold, the metric tensor is an inescapable concept.
Commonly denoted $g_{ij}$, the metric tensor defines the warping of space by
defining the distance $ds$ that is traversed when the co-ordinates are
infinitesimally incremented by some adjustment $d\boldsymbol{\theta}$,

$ds^{2}=\sum_{i}\sum_{j}g_{ij}\,d\theta_{i}d\theta_{j}=d\boldsymbol{\theta}^{T}\mathbf{G}d\boldsymbol{\theta},$
(1)

with $\theta_{i}$ denoting the $i$-th co-ordinate of $\boldsymbol{\theta}$,
and $\mathbf{G}$ the matrix with elements $g_{ij}$. When $\mathbf{G}$ is the
identity matrix, we have
$ds^{2}=d\boldsymbol{\theta}^{T}d\boldsymbol{\theta}=\|d\boldsymbol{\theta}\|_{2}^{2}$
and so recover the familiar Euclidean distance for a non-warped space. On the
other hand, where $\mathbf{G}$ is itself a function of the parameters, as it
is for the MVN statistical manifold, the distance associated with
$d\boldsymbol{\theta}$ depends also on the current position on this manifold.

Information geometry typically imposes this curvature on the parameters of a
family of distributions by choosing $\mathbf{G}$ to be the Fisher information
matrix associated with that distribution family [15]. For a family with
density functions $f(\mathbf{x};\boldsymbol{\theta})$, the Fisher information
matrix ${\cal I}$ has elements defined by

${\cal I}_{ij}=E_{\mathbf{X}}\left[\frac{\partial}{\partial\theta_{i}}\log
f(\mathbf{x};\boldsymbol{\theta})\,\frac{\partial}{\partial\theta_{j}}\log
f(\mathbf{x};\boldsymbol{\theta})\right].$

For $d$-variate MVNs with log-density

$\log f(\mathbf{x};\boldsymbol{\mu},\boldsymbol{\Sigma})=-\frac{d}{2}\log
2\pi-\frac{1}{2}\log\det\boldsymbol{\Sigma}-\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu})^{T}\boldsymbol{\Sigma}^{-1}(\mathbf{x}-\boldsymbol{\mu}),$

we may show that the distance metric takes the form [28]

$ds^{2}=d\boldsymbol{\mu}^{T}\boldsymbol{\Sigma}^{-1}d\boldsymbol{\mu}+\frac{1}{2}\operatorname{tr}\Bigl{[}(\boldsymbol{\Sigma}^{-1}d\boldsymbol{\Sigma})^{2}\Bigr{]}=\left\|\boldsymbol{\Sigma}^{-1/2}d\boldsymbol{\mu}\right\|^{2}+\frac{1}{2}\left\|\boldsymbol{\Sigma}^{-1/2}d\boldsymbol{\Sigma}\boldsymbol{\Sigma}^{-1/2}\right\|_{F}^{2},$
(2)

where $\|\cdot\|_{F}$ denotes the Frobenius norm. Equation (2) completely
defines the geometry of the Fisher manifold for MVNs, but only for two points
infinitesimally close together (separated by $\boldsymbol{\mu}$ and
$\boldsymbol{\Sigma}$), and it is not immediately obvious what minimum
distance paths between finitely-separated points should look like. Before we
discuss such paths, however, it is beneficial to gain a high-level
understanding of how the warping of space defined by equation (2) acts.

The first thing we may observe is that for any point on the manifold
$\boldsymbol{\theta}=(\boldsymbol{\mu},\boldsymbol{\Sigma})$, the distance is
affected only by $\boldsymbol{\Sigma}$ (which appears only as the precision,
$\boldsymbol{\Sigma}^{-1}$). This makes good statistical sense, in that
different values of $\boldsymbol{\mu}$ do not inherently correspond to
differing amounts of information, only to information pointing to different
locations in the sample space for the associated random variable. The
precision matrix, however, describes how well $\boldsymbol{\mu}$ is specified,
and thus has direct effect on how the distribution will respond to additional
information. That is, the greater the precision, the more information that is
required to obtain a shift in the distribution’s parameters. This information
cost is composed of wholly separate contributions due to the movement in
$\boldsymbol{\mu}$-spaec and the movement in $\boldsymbol{\Sigma}$-space,
owing to the block-diagonal structure of the Fisher information matrix in this
case [29].

The second piece of intuition comes from the norm-based definition for
$ds^{2}$ in equation (2). Considering the property for an MVN-distributed
random variable $\mathbf{x}\sim{\cal N}(\boldsymbol{\mu},\boldsymbol{\Sigma})$
that the affine-transformed variable
$\mathbf{y}=\mathbf{P}\mathbf{x}+\mathbf{r}$ is distributed
$\mathbf{y}\sim{\cal
N}(\mathbf{P}\boldsymbol{\mu}+\mathbf{r},\mathbf{P}\boldsymbol{\Sigma}\mathbf{P}^{T})$,
we describe affine transformations of an MVN’s parameters as those of form

$\boldsymbol{\mu}\mapsto\mathbf{P}\boldsymbol{\mu}+\mathbf{r},\qquad\qquad\boldsymbol{\Sigma}\mapsto\mathbf{P}\boldsymbol{\Sigma}\mathbf{P}^{T},\qquad\qquad\det(\mathbf{P})>0.$
(3)

We then see through the second expression in equation (2) that the elements of
$d\boldsymbol{\theta}=(d\boldsymbol{\mu},d\boldsymbol{\Sigma})$ do in fact add
as a sum of squares reminiscent of Euclidean distance, but only after the
effect of the current precision is removed by applying the affine
transformation $\mathbf{P}=\boldsymbol{\Sigma}^{-1/2}$,
$\mathbf{r}=\mathbf{0}$ that achieves an identity covariance (known as a
sphering or whitening transformation). Indeed, the well-known Mahalanobis
distance [30] may be thought of as the Euclidean distance after application of
a whitening transformation, and the first term in (2) is precisely the
(squared) Mahalanobis distance between two means separated by
$d\boldsymbol{\mu}$. This crystallises the property of the MVN statistical
manifold that shifts in the mean in strongly-specified directions cost more
information than shifts in weakly-specified directions.

It can also be shown via direct substitution that the metric tensor for MVNs
(2) is invariant to further affine transformations of the form (3).
Interpreted through the lens of the norm-based definition, this invariance may
be intuitively understood as pure translations of the mean (through
$\mathbf{r}$) not affecting $d\boldsymbol{\mu}$, and the sphering
transformation applied in (3) removing any scaling effects imparted by
$\mathbf{P}$ (rotation may still occur). This invariance of the metric is an
important property, as for a given pair of points, we may first transform
these points so that one lies at the origin,
$(\boldsymbol{\mu},\boldsymbol{\Sigma})=(\mathbf{0},\mathbf{I})$, without
affecting the lengths or intrinsic shapes of these paths [31] between them.

###  2.2 Geodesics in Multivariate Normal Distribution Space

The distance metric (2) directly specifies infinitesimal distances on the
statistical manifold associated with MVNs. To define finite distances on this
manifold, these infinitesimal increments must be integrated along the path of
minimal distance that connects them — the geodesic. In Euclidean space, the
geodesics are straight lines, and thus the Euclidean metric functions for both
infinitesimal and finite distances. For MVNs, however, the geodesics are
nontrivial and unavailable in closed form. When found, however, their length
gives the Fisher–Rao distance between a pair of MVNs, a statistically-informed
measure of their separation. For very close members of a distribution family
(MVN or otherwise), half of the square of the Fisher–Rao distance approximates
the KL divergence [32]. For more distant distributions, the two become
distinct.

Before describing the geodesics for the space of MVNs mathematically, we
present in Fig. 1 some example geodesics in the space of bivariate and
trivariate normal distributions that illustrate how the metric (2) influences
their shapes. As discussed above, only the covariance matrix
$\boldsymbol{\Sigma}$ contributes to the warping of space, and so we may
interpret $\boldsymbol{\mu}$-space as functioning in a Euclidean fashion
outside of this “external” influence. As such, where we discuss “straight” or
“curved” paths in $\boldsymbol{\mu}$-space we use these terms in the Euclidean
sense (in the Riemannian sense, the geodesics function as the analog of
straight lines by definition). Following others [23, 33], we also plot
bivariate and trivariate MVNs and their geodesics in this sense, as points and
paths through $\boldsymbol{\mu}$-space with covariance indicated by
ellipses/ellipsoids with axis lengths given by (a scaling of) the eigenvalues
of $\boldsymbol{\Sigma}$, and orientations given by its eigenvectors.

Geodesic paths are seen to curve in $\boldsymbol{\mu}$-space in order to
strike the balance between moving towards the target $\boldsymbol{\mu}$ value
in the Euclidean sense, and travelling along higher-variance directions (which
as per equation (2) incur less distance). In general, the geodesics initially
increase the variance in the direction of travel, easing their movement
through $\boldsymbol{\mu}$-space before restoring the covariance to that of
the curve’s endpoint. The endpoint’s covariance can also bias the direction in
$\boldsymbol{\mu}$-space from which the geodesic approaches, providing a
motivation to leave the starting point from a direction that does not
immediately point at the target. Notably however, when a principal direction
of variance (as defined by the covariance matrix’s eigenvectors) is aligned
with the Euclidean connection between the two points in
$\boldsymbol{\mu}$-space, the geodesic travels along this path. We return to
this point in more detail in Section 2.3.

Figure 1: Example geodesics starting from the origin for bivariate (left) and
trivariate (right) MVNs. Drawn paths indicate movement through
$\boldsymbol{\mu}$-space, while positions in $\boldsymbol{\Sigma}$-space at
regular intervals along this path are indicated by (scaled)
ellipses/ellipsoids. Geodesic paths curve through $\boldsymbol{\mu}$-space to
optimise movement in higher-variance directions and often increase variance in
the direction of travel.

Geodesics on Riemannian manifolds are found by solving the geodesic equations,
which are ordinary differential equations (ODEs) defined wholly by the metric
tensor. For the space of MVN distributions, these equations take the neat form
[31]

$\displaystyle\frac{d^{2}\boldsymbol{\mu}}{dt^{2}}$
$\displaystyle=\frac{d\boldsymbol{\Sigma}}{dt}\,\boldsymbol{\Sigma}^{-1}\,\frac{d\boldsymbol{\mu}}{dt}$
(4) $\displaystyle\frac{d^{2}\boldsymbol{\Sigma}}{dt^{2}}$
$\displaystyle=\frac{d\boldsymbol{\Sigma}}{dt}\boldsymbol{\Sigma}^{-1}\frac{d\boldsymbol{\Sigma}}{dt}-\left(\frac{d\boldsymbol{\mu}}{dt}\right)\left(\frac{d\boldsymbol{\mu}}{dt}\right)^{T}\\!\\!.$

In order to simplify the geodesic equations, it is beneficial to consider
geodesics originating from the origin,
$(\boldsymbol{\mu},\boldsymbol{\Sigma})=(\mathbf{0},\mathbf{I})$. Making use
of the freedom to transform according to (3) without changing the metric, we
may choose

$\mathbf{P}=\boldsymbol{\Sigma}^{1/2}\mathbf{R},\qquad\mathbf{r}=\boldsymbol{\mu}$
(5)

to map the origin, $(\mathbf{0},\mathbf{I})$, to any given point
$(\boldsymbol{\mu},\boldsymbol{\Sigma})$. Here, $\mathbf{R}$ is an orthogonal
matrix $\mathbf{R}\mathbf{R}^{T}=\mathbf{I}$ that reflects the freedom to
choose $\mathbf{P}$ as any decomposition
$\mathbf{P}\mathbf{P}^{T}=\boldsymbol{\Sigma}$. This is equivalent to the
freedom to choose from a variety of transformation matrices that all carry out
the process of data whitening [34].

Considering (5), we see that the geodesic between a starting point,
$(\boldsymbol{\mu}_{0},\boldsymbol{\Sigma}_{0})$, and a target point,
$(\boldsymbol{\mu}_{1},\boldsymbol{\Sigma}_{1})$, can instead be treated as a
connection between the origin and a new transformed target,

$(\boldsymbol{\mu}_{t},\boldsymbol{\Sigma}_{t})=\Bigl{(}\mathbf{R}^{T}\boldsymbol{\Sigma}_{0}^{-1/2}(\boldsymbol{\mu}_{1}-\boldsymbol{\mu}_{0}),\,\,\mathbf{R}^{T}\boldsymbol{\Sigma}_{0}^{-1/2}\boldsymbol{\Sigma}_{1}\boldsymbol{\Sigma}_{0}^{-1/2}\mathbf{R}\Bigr{)}.$
(6)

The desired geodesic can then be obtained by applying the transform (5) to
this connection according to (3). For a geodesic beginning at the origin, the
geodesic equations (4) can be integrated once to become

$\displaystyle\frac{d\boldsymbol{\mu}}{dt}$
$\displaystyle=\boldsymbol{\Sigma}\mathbf{x}$
$\displaystyle\mathbf{x}=\left.\frac{d\boldsymbol{\mu}}{dt}\right|_{t=0}$ (7)
$\displaystyle\frac{d\boldsymbol{\Sigma}}{dt}$
$\displaystyle=\boldsymbol{\Sigma}(\mathbf{B}-\mathbf{x}\boldsymbol{\mu}^{T})$
$\displaystyle\mathbf{B}=\left.\frac{d\boldsymbol{\Sigma}}{dt}\right|_{t=0}.$

Several authors [20, 31, 35] have shown that these equations permit a closed-
form solution when expressed in terms of the canonical co-ordinates
$(\boldsymbol{\delta},\boldsymbol{\Delta})=(\boldsymbol{\Sigma}^{-1}\boldsymbol{\mu},\boldsymbol{\Sigma}^{-1})$,

$\displaystyle\boldsymbol{\delta}(t)$
$\displaystyle=-\mathbf{B}\Bigl{(}\cosh(t\mathbf{G})-\mathbf{I}\Bigr{)}(\mathbf{G}^{-})^{2}\mathbf{x}+\sinh(t\mathbf{G})\mathbf{G}^{-}\mathbf{x},$
(8) $\displaystyle\boldsymbol{\Delta}(t)$
$\displaystyle=\mathbf{I}+\frac{1}{2}\Bigl{(}\cosh(t\mathbf{G})-\mathbf{I}\Bigr{)}+\frac{1}{2}\mathbf{B}\Bigl{(}\cosh(t\mathbf{G})-\mathbf{I}\Bigr{)}(\mathbf{G}^{-})^{2}\mathbf{B}-\frac{1}{2}\sinh(t\mathbf{G})\mathbf{G}^{-}\mathbf{B}-\frac{1}{2}\mathbf{B}\sinh(t\mathbf{G})\mathbf{G}^{-}.$

Here $\mathbf{G}^{2}=\mathbf{B}^{2}+2\mathbf{x}\mathbf{x}^{T}$, and
$\mathbf{G}^{-}$ denotes the generalised inverse of $\mathbf{G}$, satisfying
$\mathbf{B}\mathbf{G}\mathbf{G}^{-}=\mathbf{G}\mathbf{G}^{-}\mathbf{B}=\mathbf{B}$
and $\mathbf{G}\mathbf{G}^{-}\mathbf{x}=\mathbf{G}^{-}\mathbf{G}\mathbf{x}$
(see [31] for details). Notably, owing to the forms of the (globally
convergent) Taylor series for the hyperbolic sine and cosine functions, it can
be seen that equations (8) calls for only even powers of $\mathbf{G}$, and so
there is no degeneracy associated with knowing only $\mathbf{G}^{2}$, and not
$\mathbf{G}$.

Although equations (8) can be evaluated directly, in our experience points
along the geodesic path are best numerically evaluated using the form of the
solution presented in [35]. In this form, $\boldsymbol{\delta}(t)$ and
$\boldsymbol{\Delta}(t)$ are found by reading off a portion of the result of
exponentiating a matrix of augmented dimension,

$\begin{pmatrix}\boldsymbol{\Delta}&\boldsymbol{\delta}&\cdot\,\,\\\
\boldsymbol{\delta}^{T}&\cdot&\cdot\,\,\\\
\cdot&\cdot&\cdot\,\,\end{pmatrix}=\exp\left(t\begin{pmatrix}-\mathbf{B}&\mathbf{x}&\mathbf{0}\\\
\mathbf{x}^{T}&0&-\mathbf{x}^{T}\\\
\mathbf{0}&-\mathbf{x}&\mathbf{B}\end{pmatrix}\right).$ (9)

Here elements of the matrix marked with $\cdot$ are unused.

By observing that the velocity along geodesic curves is constant, the length
of these curves over the segment $t\in[0,1]$ (the Fisher–Rao distance between
the distributions at their start and endpoints) can be calculated by
evaluating the velocity at the origin according to the metric (2), giving [33]

$d_{F}=\sqrt{\frac{1}{2}\operatorname{tr}(\mathbf{B}^{2})+\mathbf{x}^{T}\mathbf{x}}=\sqrt{\frac{1}{2}\operatorname{tr}(\mathbf{G}^{2})}\,.$
(10)

Given a geodesic’s initial velocity $\mathbf{v}=(\mathbf{x},\mathbf{B})$, we
thereby have closed-form expressions for both the path the geodesic takes (8),
and the Fisher–Rao distance to the point at $t=1$ (10). However, as yet there
is no known analytic solution to the problem of choosing $\mathbf{v}$ so as to
produce a geodesic passing through a given arbitrary target point. Outside of
the special cases discussed subsequently, geodesics between arbitrary points
in the MVN space can be found using a specially-tailored shooting method [23],
as we discuss in Section 4. If only the Fisher–Rao distance between a pair of
MVNs is required, closed-form bounds may suffice [27, 36].

###  2.3 Special Cases with Known Solution

####  2.3.1 Regions in $\boldsymbol{\mu}$-space Admitting Closed-Form Solution

For several special cases, the velocity of the geodesic connecting a pair of
MVNs is known, as summarised by Pinele et al. [33] and also Nielsen [27] using
the lens of totally geodesic submanifolds. Also relevant is the work of
Herntier and Peter [28], who consider solutions for closest points on
submanifolds (and geodesic paths to them), given some starting point on the
full statistical manifold for bivariate MVNs. Here, we present the cases with
known closed-form solution without any further invocation of differential
geometry.

As equation (8) gives the form of geodesics that begin at the origin, we
discuss the cases with known solution in the context of having already been
affine transformed so that they begin at the origin. That is, rather than the
problem of finding the geodesic connecting two MVNs
$(\boldsymbol{\mu}_{0},\boldsymbol{\Sigma}_{0})$ and
$(\boldsymbol{\mu}_{1},\boldsymbol{\Sigma}_{1})$ we consider the analogous
problem of finding the geodesic connecting the pair $(\mathbf{0},\mathbf{I})$
and $(\boldsymbol{\mu}_{t},\boldsymbol{\Sigma}_{t})$, with the latter given by
equation (6). The desired initial velocity in this transformed space is
$\mathbf{v}=(\mathbf{x},\mathbf{B})$, and when it is obtained, the velocity of
the geodesic connecting the actual targets is given by undoing the transform,

$\left.\dot{\boldsymbol{\mu}}\right|_{t=0}=\boldsymbol{\Sigma}_{0}^{1/2}\mathbf{R}\mathbf{x},\qquad\qquad\left.\dot{\boldsymbol{\Sigma}}\right|_{t=0}=\boldsymbol{\Sigma}_{0}^{1/2}\mathbf{R}\mathbf{B}\mathbf{R}^{T}\boldsymbol{\Sigma}_{0}^{1/2}.$

Owing to the invariance of the metric to these affine transforms, the
Fisher–Rao distance between the original pair and the transformed pair remains
the same.

We provide the cases with known solution in Table 1, and the associated
Fisher–Rao distances in Table 2. Although these special cases may appear
distinct, they are all in fact examples of a single type of solution. To
demonstrate this, we observe that the choice of $\mathbf{R}$ in equation (6)
is free, and so we may choose $\mathbf{R}$ so that it diagonalises
$\boldsymbol{\Sigma}_{0}^{-1/2}\boldsymbol{\Sigma}_{1}\boldsymbol{\Sigma}_{0}^{-1/2}$
and the target covariance becomes diagonal,
$\boldsymbol{\Sigma}_{t}=\mathbf{D}$. Given this, applicability of the “axis-
aligned targets” case in Tables 1 and 2 depends only on its second condition,
that the precision-weighted target mean,
$\boldsymbol{\delta}_{t}=\boldsymbol{\Sigma}^{-1}\boldsymbol{\mu}_{t}$, is
aligned with a principal axis after applying this choice of transformation,

$\boldsymbol{\Sigma}_{t}^{-1}\boldsymbol{\mu}_{t}=\delta_{t}\mathbf{e}_{k}.$

Here $\delta_{t}=\|\boldsymbol{\delta}_{t}\|_{2}$ is a scalar and
$\mathbf{e}_{k}$ an elementary basis vector with $k$-th element unity and all
other elements zero. This condition is satisfied only if
$\boldsymbol{\mu}_{t}$ is aligned with an eigenvector of
$\boldsymbol{\Sigma}_{t}$.

Case | $\mathbf{B}$ | $\mathbf{x}$
---|---|---
Univariate, | $\dfrac{\alpha}{\gamma}g$ | $\dfrac{2g\delta_{t}\Delta_{t}}{\gamma}$
$\mu_{t}=\delta_{t}/\Delta_{t},\,\,\Sigma_{t}=1/\Delta_{t}$
Equal means, $\boldsymbol{\mu}_{t}=0$ | $\log\left(\boldsymbol{\Sigma}_{t}\right)$ | $\mathbf{0}$
Equal variances, $\boldsymbol{\Sigma}_{t}=\mathbf{I}$ | $\dfrac{\operatorname{acosh}\left(1+\|\boldsymbol{\mu}_{t}\|^{2}+\frac{1}{8}\|\boldsymbol{\mu}_{t}\|^{4}\right)}{\|\boldsymbol{\mu}_{t}\|\sqrt{\|\boldsymbol{\mu}_{t}\|^{2}+8}}\,\,\boldsymbol{\mu}_{t}\boldsymbol{\mu}_{t}^{T}$ | $\dfrac{2\operatorname{acosh}\left(1+\|\boldsymbol{\mu}_{t}\|^{2}+\frac{1}{8}\|\boldsymbol{\mu}_{t}\|^{4}\right)}{\|\boldsymbol{\mu}_{t}\|\sqrt{\|\boldsymbol{\mu}_{t}\|^{2}+8}}\,\,\boldsymbol{\mu}_{t}$
Axis-aligned targets, | $\Bigl{(}\dfrac{\alpha}{\gamma}g\,\mathbf{I}-\log(\mathbf{D})\Bigr{)}\mathbf{e}_{k}\mathbf{e}_{k}^{T}+\log(\mathbf{D})$ | $\dfrac{2g\delta_{t}\Delta_{t}}{\gamma}\,\,\mathbf{e}_{k}$
$\boldsymbol{\Sigma}_{t}^{-1}\boldsymbol{\mu}_{t}=\delta_{t}\mathbf{e}_{k},\,\,\boldsymbol{\Sigma}_{t}=\mathbf{D}$
Constants:
$\Delta_{t}=1/D_{kk},\quad\alpha=\delta_{t}^{2}-2\Delta_{t}^{2}+2\Delta_{t},\quad\gamma=\sqrt{\alpha^{2}+8\delta_{t}^{2}\Delta_{t}^{2}}$,
$g=\operatorname{acosh}\left(1+\gamma^{2}/8\Delta_{t}^{3}\right)$
Table 1: Special cases for which the geodesic connecting the origin to a target point, $(\boldsymbol{\mu}_{t},\boldsymbol{\Sigma}_{t})$, is known. The target point $(\boldsymbol{\mu}_{t},\boldsymbol{\Sigma}_{t})$ associated with a given start and end point is provided by equation (6). The condition that the target mean is aligned with an elementary basis vector, $\boldsymbol{\Sigma}^{-1}_{t}\boldsymbol{\mu}_{t}=\delta_{t}\mathbf{e}_{k}$ defines the index $k$, and for the univariate case $k=1$. $\mathbf{D}$ denotes a diagonal matrix. Case | $d_{F}$
---|---
Univariate, | $\dfrac{1}{\sqrt{2}}\operatorname{acosh}\left(1+\dfrac{(\delta_{t}^{2}-2\Delta_{t}^{2}+2\Delta_{t})^{2}+8\delta_{t}^{2}\Delta_{t}^{2}}{8\Delta_{t}^{3}}\right)$
$\mu_{t}=\delta_{t}/\Delta_{t},\,\,\Sigma_{t}=1/\Delta_{t}$
Equal means, $\boldsymbol{\mu}_{t}=0$ | $\sqrt{\dfrac{1}{2}\sum_{i=1}(\log\lambda_{i})^{2}},\qquad\mbox{$\lambda_{i}$ the eigenvalues of $\boldsymbol{\Sigma}_{t}$}$
Equal variances, $\boldsymbol{\Sigma}_{t}=\mathbf{I}$ | $\dfrac{1}{\sqrt{2}}\operatorname{acosh}\left(1+\|\boldsymbol{\mu}_{t}\|^{2}+\frac{1}{8}\|\boldsymbol{\mu}_{t}\|^{4}\right)$
Axis-aligned targets, | $\sqrt{\frac{1}{2}\operatorname{acosh}^{2}\left(1+\dfrac{(\delta_{t}^{2}-2D_{kk}^{2}+2D_{kk})^{2}+8\delta_{t}^{2}D_{kk}^{2}}{8D_{kk}^{3}}\right)+\frac{1}{2}\sum\limits_{j\neq k}(\log D_{jj})^{2}}$
$\boldsymbol{\Sigma}_{t}^{-1}\boldsymbol{\mu}_{t}=\delta_{t}\mathbf{e}_{k},\,\,\boldsymbol{\Sigma}_{t}=\mathbf{D}$
Table 2: Fisher–Rao distances $d_{F}$ associated with the special cases given
in Table 1. The target point $(\boldsymbol{\mu}_{t},\boldsymbol{\Sigma}_{t})$
associated with a given start and end point is provided by equation (6). The
condition that the target mean is aligned with an elementary basis vector,
$\boldsymbol{\Sigma}^{-1}_{t}\boldsymbol{\mu}_{t}=\delta_{t}\mathbf{e}_{k}$
defines the index $k$, and for the univariate case $k=1$. $\mathbf{D}$ denotes
a diagonal matrix.

Such a condition is satisfied by all of the simpler cases in Tables 1 and 2.
For the univariate case, all 1-vectors are aligned by definition and the
condition is trivially satisfied. For the case of equal means,
$\boldsymbol{\mu}_{t}=0$ and the zero-length vector can be considered aligned
with any eigenvector of the target covariance. For the case of equal
covariance matrices, $\boldsymbol{\Sigma}_{t}=\mathbf{I}$ and choice of
$\mathbf{R}$ is not restricted to be the diagonaliser. As such, we can instead
choose $\mathbf{R}$ that rotates $\boldsymbol{\mu}_{t}$ to align with a co-
ordinate axis and hence satisfy the condition.

More precisely, the requirement for the use of this solution is that
$\boldsymbol{\mu}_{t}$ lies within the eigenspace associated with one of the
eigenvalues of $\boldsymbol{\Sigma}_{t}$. Given that the geometric and
algebraic multiplicities of all covariance matrix eigenvalues are equal, each
eigenvalue repeated $n$ times thus corresponds to a (hyper)plane of dimension
$n$ in which $\boldsymbol{\mu}_{t}$ (after rotating to diagonalise
$\boldsymbol{\Sigma}_{t}$) may lie for a closed-form solution to be available.
Fig. 2 presents a graphical demonstration of this concept, showing how in the
space of trivariate MVNs, an eigenvalue of multiplicity two and another of
multiplicity one correspond to a plane and a line, respectively, of target
means $\boldsymbol{\mu}_{t}$ permitting closed-form geodesics.

Figure 2: Spaces of target means $\mu_{1}$ reachable by “straight line”
geodesics, as imposed by the eigenstructure of
$\boldsymbol{\Sigma}_{t}=\boldsymbol{\Sigma}_{0}^{-1/2}\boldsymbol{\Sigma}_{1}\boldsymbol{\Sigma}_{0}^{-1/2}$.
Pictured are a set of example geodesics for trivariate MVNs, each starting
from the same point $(\boldsymbol{\mu}_{0},\boldsymbol{\Sigma}_{0})$ and
sharing the same $\boldsymbol{\Sigma}_{1}$. Here, $\boldsymbol{\Sigma}_{t}$
has a repeated eigenvalue (multiplicity 2) and produces a plane of
$\boldsymbol{\mu}_{1}$ values for which geodesic paths (blue) travel straight
lines through $\boldsymbol{\mu}$ space with closed form (Table 1). Geodesics
to points with $\boldsymbol{\mu}_{1}$ outside of the plane (red) only take
this form if they align with the remaining eigenvector of
$\boldsymbol{\Sigma}_{t}$ (dotted line).

The closed-form geodesics travel straight lines in $\boldsymbol{\mu}$ space.
Mathematically, these correspond to the cases where $\mathbf{B}$ and
$\mathbf{G}$ are mutually diagonalisable, allowing for each dimension in
equations (8) to be decoupled via a similarity transform (see supplementary
material). In these cases, the covariance matrix eigenvectors remain fixed in
specific directions along the geodesic, and geodesics travel along straight
lines oriented in one of these directions. This property is key to a geometric
understanding of these solutions, and helps understand the geodesics of the
MVN statistical manifold more generally.

As the metric (2) dictates, increasing the variance in a direction lowers the
information cost of moving in that direction. Aligning the covariance
eigenvectors with the direction of movement thus becomes a way to potentially
decrease overall information cost, by providing the ability to increase
variance only in the direction where this will be useful. General geodesics
must balance this effect against the orientation of eigenvectors at the
geodesic’s start and end, as well as the curved path they take through
$\boldsymbol{\mu}$-space. However, when the direction to the target already
aligns with a covariance eigenvector, there is no longer any motivation to
rotate the principal directions of variance. The geodesic solution optimises
for increasing the variance associated with that direction to ease movement,
and reaching the correct target, while all other eigenvalues have no movement
through $\boldsymbol{\mu}$-space to consider and simply transform according to
the equal means solution, such that $\log\lambda$ linearly approaches the
target covariance eigenvalue $\lambda_{t}$,

$\lambda(t)=(\lambda_{t})^{t}.$ (11)

Even when the direction to the target in $\boldsymbol{\mu}$-space does not
align with an eigenvector of the target covariance, equation (11) still serves
as a lower bound for the values of covariance eigenvalues even while their
associated eigenvectors shift direction (Fig. 3). This is geometrically
understood by considering that equation (11) represents the lowest-cost path
eigenvalues follow to reach a given target, ignoring influence from any
movement through $\boldsymbol{\mu}$-space. When we do consider that movement,
there may be a motivation to increase variance in that direction in order to
lower the information cost associated with it. However, there is never any
reason to decrease variance below the path defined by equation (11).
Interestingly, this result emerges naturally from our geometric understanding,
but seems far more difficult to demonstrate mathematically.

Figure 3: Paths traced out by the larger (blue) and smaller (red) eigenvalue
of $\boldsymbol{\Sigma}$, for geodesics in the space of bivariate MVNs
connecting the origin to random points with final covariance eigenvalues
$\log\lambda_{1}=1$, $\log\lambda_{2}=-1$. The straight dotted lines depict
$\log\lambda=(\log\lambda_{t})t$, the lower bound for eigenvalues along all
geodesic paths.

##  3 Non-Geodesic Paths

###  3.1 Commonly-Used Paths

Given the observation that paths with low information cost improve samplers
using distribution sequences [9, 10], we summarise here how commonly-used
paths in such settings function in terms of information cost, using our
distribution family of MVNs as an accessible way to understand how these
different choices of paths behave. Some of these paths along the MVN manifold
were discussed in a very recent publication [27], however in keeping with the
character of this review we concentrate on practical performance and an
interpretation external to information geometry. For notational consistency
with the geodesic paths already discussed, we parameterise these paths using
$t\in[0,1]$, with $t=0$ corresponding to the starting point
$\boldsymbol{\theta}_{0}=(\boldsymbol{\mu}_{0},\boldsymbol{\Sigma}_{0})$ and
$t=1$ to the endpoint,
$\boldsymbol{\theta}_{1}=(\boldsymbol{\mu}_{1},\boldsymbol{\Sigma}_{1})$.

A predominant choice in a variety of contexts is the annealing path (e.g. [37,
38]), which takes the form

$p(t;\mathbf{x})=\bigl{[}p_{0}(\mathbf{x})\bigr{]}^{1-t}\,\bigl{[}p_{1}(\mathbf{x})\bigr{]}^{t},\qquad\quad\mbox{equivalently}\quad\begin{aligned}
\boldsymbol{\delta}(t)&=(1-t)\boldsymbol{\delta}_{0}+t\boldsymbol{\delta}_{1}\\\
\boldsymbol{\Delta}(t)&=(1-t)\boldsymbol{\Delta}_{0}+t\boldsymbol{\Delta}_{1}.\end{aligned}$

This path represents a linear interpolation between the log densities
associated with the two distributions, and as MVNs are a member of the
exponential family, the annealing path also corresponds to a linear
interpolation in terms of their natural parameters [12], and equivalently,
$\boldsymbol{\delta}$ and $\boldsymbol{\Delta}$. Observing that the geodesic
paths (8) are non-linear in these canonical co-ordinates immediately
demonstrates that annealing paths are non-geodesic.

A second choice of path for distributions in a single exponential family is
the moment-averaged path, which can offer significantly improved performance
in annealed importance sampling [12]. For MVNs, this path takes the form

$\displaystyle\boldsymbol{\mu}(t)$
$\displaystyle=(1-t)\boldsymbol{\mu}_{0}+t\boldsymbol{\mu}_{1}$ (12)
$\displaystyle\boldsymbol{\Sigma}(t)$
$\displaystyle=(1-t)\boldsymbol{\Sigma}_{0}+t\boldsymbol{\Sigma}_{1}+\frac{1}{4}t(1-t)(\boldsymbol{\mu}_{1}-\boldsymbol{\mu}_{0})(\boldsymbol{\mu}_{1}-\boldsymbol{\mu}_{0})^{T}.$

These paths move in Euclidean fashion through $\boldsymbol{\mu}$-space, which
immediately renders them non-geodesic. However, these paths hold the
attractive property that this movement in $\boldsymbol{\mu}$-space is eased by
the temporary injection of variance in that direction, mimicking one of the
key behaviours of geodesic paths.

A somewhat separate type of path are those defined by the Wasserstein
distance. The Wasserstein distance is formally a distance with an associated
metric [6], albeit one of different form to that of the Fisher metric.
Informally, the Wasserstein distance measures the amount of physical work
required to move the probability mass of one density function, in the fashion
that minimises the work performed, so that it becomes the density function of
another. Paths that are geodesic with respect to the Wasserstein metric are
not in general geodesic with respect to the Fisher metric, but are interesting
to consider given the use of Wasserstein distance as a statistical distance
measure [4] and in averaging distributions via the Fréchet mean [39].

In the space of MVNs, the 2-Wasserstein distance is typically used as it
provides a closed form for the geodesic path (and distance) between two MVNs
[40]. This path takes the form

$\displaystyle\boldsymbol{\mu}(t)$
$\displaystyle=(1-t)\boldsymbol{\mu}_{0}+t\boldsymbol{\mu}_{1}$ (13)
$\displaystyle\boldsymbol{\Sigma}(t)$
$\displaystyle=\Bigl{(}(1-t)\mathbf{I}+t\mathbf{C}\Bigr{)}\boldsymbol{\Sigma}_{0}\Bigl{(}(1-t)\mathbf{I}+t\mathbf{C}\Bigr{)},\qquad\mathbf{C}=\boldsymbol{\Sigma}_{1}^{1/2}\left(\boldsymbol{\Sigma}_{1}^{1/2}\boldsymbol{\Sigma}_{0}\boldsymbol{\Sigma}_{1}^{1/2}\right)^{-1/2}\boldsymbol{\Sigma}_{1}^{1/2}.$

Paths minimising the 2-Wasserstein distance also move in Euclidean fashion
through $\boldsymbol{\mu}$-space, but without any explicit injection of
variance — equations (13) may be interpreted as applying a linearly-varying
affine transformation to the starting point of the path, starting with the
identity transformation and ending with an affine transformation that maps
$(\boldsymbol{\mu}_{0},\boldsymbol{\Sigma}_{0})$ to
$(\boldsymbol{\mu}_{1},\boldsymbol{\Sigma}_{1})$. As such, the Wasserstein
metric does not directly acknowledge the information about the random variable
contained in the distribution [41], and indeed for the univariate case
equations (13) reduce to simply a linear interpolation of mean and standard
deviation, $\sqrt{\Sigma(t)}=(1-t)\sqrt{\Sigma_{0}}+t\sqrt{\Sigma_{1}}$ [3].
Despite this information-unaware property, the Wasserstein distance has
numerous applications in statistics, many of which are reviewed by Panaretos
and Zemel [6]. It is also well suited to comparing distributions with disjoint
support, in contrast to information theoretic approaches [41].

In Fig. 4a, we visualise these different choices of closed-form path for a few
examples in the space of bivariate MVNs. We immediately see that
2-Wasserstein-optimal paths move in direct fashion through the space, and are
perhaps best interpreted as optimal absent the informational context.
Annealing paths strongly favour travelling the directions of high variance at
the start and end point, as they do not inject additional variance. These
paths curve too much in $\boldsymbol{\mu}$-space relative to the geodesic path
and hence become information-inefficient. Conversely, moment-averaged paths
inject additional variance to ease movement through $\boldsymbol{\mu}$-space,
but do not curve this movement to take advantage of the variance at the path’s
endpoints.

As observed by Grosse et al. [12], the behaviours of annealing and moment-
averaged paths can be understood through a variational lens, as the former can
be shown to be reverse-KL optimal (mode-seeking, avoiding additional variance)
and the latter forwards-KL optimal (mean-seeking, adding additional variance).
Nielsen [27] demonstrates the potential of paths formed by combining these two
paths, trading off their respective strengths and weaknesses.

a)

b)

Figure 4:  a) The paths taken through the space of bivariate MVNs, and
associated manifold distances, of commonly-used non-geodesic paths for the
examples shown in Fig. 1. b) Performance (in terms of manifold distance
relative to the geodesic distance $d_{F}$) for different choices of path
connecting the origin to a large number of random points in the space of MVNs.
Moment-averaged paths prove superior.

To test each type of path’s performance more generally, we generate a large
number of random points in MVN space by independently sampling components of
$\boldsymbol{\mu}$ uniformly, $\mu_{i}\sim U[-10,10]$, sampling covariance
eigenvalues log-uniformly, $\log\lambda\sim U[-10,10]$, and creating
covariance matrices from these eigenvalues by generating random orthonormal
matrices through the Gram–Schmidt process applied to randomly selected
direction vectors. We connect each of these points to the origin with the
different types of paths, and calculate their length by summing the lengths of
each of their segments (2) after applying a fine discretisation. Boxplots of
these lengths across the many test problems are given in Fig. 4b. Despite
taking very differently-shaped paths through the space of MVNs, annealing
paths and 2-Wasserstein-optimal paths perform extremely similarly in terms of
the information distance associated with their traversal, at least for test
problems generated in this manner. Moment-averaged paths are much better at
travelling close-to-information-optimal paths through MVN space, in accordance
with their demonstrated ability to improve annealed importance sampling [12].

###  3.2 MVN-Specific Approximations

Ideas that approximate the geodesic path, either by approximating the correct
choice of $\mathbf{v}$ to reach a given target or by generating paths that
reach the target and are close to geodesic, are useful both for approximating
or bounding the Fisher–Rao distance [36, 27], and for improving numerical
efforts to find geodesic paths (Section 4). Here we briefly discuss both types
of approximation, in all cases considering the case where equation (6) has
already been used to convert the problem of connecting two MVNs to the
equivalent problem of connecting the origin to a given target,
$\boldsymbol{\theta}_{t}=(\boldsymbol{\mu}_{t},\boldsymbol{\Sigma}_{t})$.

Nielsen [27] observes that approximately-geodesic paths can be found by
projecting from a higher-dimensional Riemannian manifold for symmetric
positive definite matrices, ${\cal P}$, for which Calvo and Oller [20]
demonstrated a submanifold is isometric to the MVN Riemannian manifold.
Geodesics connecting the origin to where the target point is embedded in
${\cal P}$ are given simply by

$\mathbf{p}(t)=\begin{bmatrix}\boldsymbol{\Sigma}_{t}+\boldsymbol{\mu}_{t}\boldsymbol{\mu}_{t}^{T}&\boldsymbol{\mu}_{t}\\\
\boldsymbol{\mu}_{t}^{T}&1\end{bmatrix}^{t}=\begin{bmatrix}p_{\boldsymbol{\Sigma}}(t)&p_{\boldsymbol{\mu}}(t)\\\
p_{\boldsymbol{\mu}}(t)^{T}&p_{\beta}(t)\end{bmatrix},$ (14)

where the latter matrix is defined purely for notational convenience. The sub-
manifold corresponding to MVNs requires $p_{\beta}=1$, but $\mathbf{p}(t)$
violates this for intermediate values of $0<t<1$ and so leaves the MVN
manifold. Nevertheless, approximately-geodesic paths through MVN space can be
obtained by projecting back onto the MVN manifold [27], via

$\boldsymbol{\mu}(t)=\frac{p_{\boldsymbol{\mu}}(t)}{p_{\beta}(t)},\qquad\boldsymbol{\Sigma}(t)=\mathbf{p}_{\boldsymbol{\Sigma}}(t)-\frac{p_{\boldsymbol{\mu}}(t)p_{\boldsymbol{\mu}}(t)^{T}}{\mathbf{p}_{\beta}(t)}.$
(15)

These paths are generally highly competitive, although Nielsen [27] observed
in numerical experiments that for higher dimensions ($d>11$) moment-averaged
paths (12) became superior. We have also observed this to be the case for
geodesics between points separated very far in (variance-weighted)
$\boldsymbol{\mu}$-space, even in lower dimensions.

We may also approximate geodesics in the sense of finding choices for
$\mathbf{v}=(\mathbf{x},\mathbf{B})$ that approximate those of the desired
geodesic, and these prove useful for numerically finding the true geodesic.
Some of these approximations can be obtained by extending the special cases in
Table 1 beyond their range of application. The “equal means” solution
corresponds to the case $\mathbf{B}=\mathbf{G}$ and as detailed in the
supplement, we gain an approximation for small $\mathbf{x}$,

$\displaystyle\mathbf{x}$
$\displaystyle\approx\Bigl{(}\mathbf{I}-(\mathbf{I}-\boldsymbol{\Sigma}_{t})^{-1}\Bigr{)}\log(\boldsymbol{\Sigma}_{t})\boldsymbol{\Sigma}^{-1}_{t}\boldsymbol{\mu}_{t}$
(16) $\displaystyle\mathbf{B}$
$\displaystyle\approx\log\boldsymbol{\Sigma}_{t},$

by substituting $\mathbf{B}^{2}\approx\mathbf{G}^{2}$ into equations (8).

We also obtain an approximate solution by relaxing the requirement that
$\boldsymbol{\mu}_{t}$ be aligned with a covariance matrix eigenvector.
Rather, we project the target mean onto each eigenvector separately, form the
optimal univariate solution for each component and combine them,

$\mathbf{B}\approx\sum_{j=1}^{d}\frac{\alpha_{j}}{\gamma_{j}}g_{j}\mathbf{v}_{j}\mathbf{v}_{j}^{T},\qquad\quad\mathbf{x}\approx\sum_{j=1}^{d}\operatorname{sign}(\mathbf{v}_{j}^{T}\boldsymbol{\delta}_{t})g_{j}\sqrt{\frac{1}{2}\left(1-\left(\frac{\alpha_{j}}{\gamma_{j}}\right)^{2}\right)}\mathbf{v}_{j},$
(17)

where $(\lambda_{j},\mathbf{v}_{j})$ are the eigenpairs of
$\boldsymbol{\Sigma}_{t}^{-1}$,
$\boldsymbol{\delta}_{t}=\boldsymbol{\Sigma}_{t}^{-1}\boldsymbol{\mu}_{t}$ as
before, and

$\alpha_{j}=(\mathbf{v}_{j}^{T}\delta)^{2}-2\lambda_{j}^{2}+2\lambda_{j},\quad\gamma_{j}=\sqrt{\alpha_{j}^{2}+8(\mathbf{v}_{j}^{T}\boldsymbol{\delta}_{t})\lambda_{j}^{2}},\quad
g_{j}=\operatorname{acosh}\left(1+\gamma_{j}^{2}/8\lambda_{j}^{3}\right).$

This approximation is exact for any of the cases in Table 1. A related idea
was used to bound the Fisher–Rao distance by Strapasson et al. [36].

Finally, we may also generate approximate initial velocities by considering
the initial velocity of the geodesics (14) on ${\cal P}$. As we demonstrate in
the supplementary material, some blocks of the logarithm of the matrix in
equation (14) provide an approximation for the required $\mathbf{x}$ and
$\mathbf{B}$,

$\mathbf{p}^{\prime}(0)=\log\begin{bmatrix}\boldsymbol{\Sigma}_{t}+\boldsymbol{\mu}_{t}\boldsymbol{\mu}_{t}^{T}&\boldsymbol{\mu}_{t}\\\
\boldsymbol{\mu}_{t}^{T}&1\end{bmatrix}\approx\begin{bmatrix}\mathbf{B}&\mathbf{x}\\\
\mathbf{x}^{T}&\cdot\end{bmatrix}.$ (18)

As with the small-$\mathbf{x}$ approximation above, this approximation is
exact in the case of equal means, as can be seen by direct substitution of
$\mathbf{x}=\mathbf{0}$ and $\boldsymbol{\mu}_{t}=\mathbf{0}$.

Our focus is on the use of these approximations inside the shooting algorithms
discussed subsequently. However, in the supplementary we display a few example
geodesics and the corresponding approximate geodesics (Fig. S1) and their
general performance across randomly-generated test problems (Fig. S2). The
projection-based approximation (18) emerges as the best overall performer,
especially in higher dimensions.

##  4 Shooting to Find Geodesic Paths

###  4.1 Fundamental Concepts

Outside the cases listed in Table 1, the choice of $\mathbf{x}$ and
$\mathbf{B}$ in equation (8) that produce the geodesic that reaches a given
target point remains unknown. Instead, geodesics between arbitrary points in
MVN space are produced numerically, via shooting [23]. Shooting is the process
of taking an initial guess for the velocity
$\mathbf{v}=(\mathbf{x},\mathbf{B})$, firing out the geodesic with this
velocity, and then adjusting this velocity according to where the geodesic ray
ends up relative to the target point. These adjustments are made leveraging
the known geometric properties of the space. Convergence of shooting is
nontrivial, and Han and Park [23] suggest that for manifold distances greater
than about seven the algorithm demonstrates issues with both stability and
convergence rate. For these more challenging problems, they propose an
alternative shooting via path refinement (Section 4.3).

Here, we first review the fundamental pieces of the algorithm of Han and Park
[23] to complete this work’s differential geometry-free treatment of these
statistical manifolds. Then, in Section 4.2, we demonstrate the use of MVN-
specific approximations (Section 3.2) inside the shooting process.

After shifting a problem to the origin using equation (5), we can obtain the
geodesic associated with a given initial velocity using equation (9). To
formulate some correction to this velocity, we require a “residual” vector
that compares the geodesic’s endpoint and the target point. This vector
becomes an update to the geodesic’s initial velocity via the differential
geometric process known as parallel transport, which moves a vector along a
Riemannian manifold whilst transforming it so as to keep it pointing in the
same direction relative to the manifold.

The parallel transport equations for the case of MVNs are given in Han and
Park [23], and for geodesics beginning at the origin they take the form

$\displaystyle\frac{d\mathbf{u}_{\boldsymbol{\mu}}}{dt}$
$\displaystyle=\frac{1}{2}\boldsymbol{\Sigma}\left(\mathbf{B}-\mathbf{x}\boldsymbol{\mu}^{T}\right)\boldsymbol{\Sigma}^{-1}\mathbf{u}_{\boldsymbol{\mu}}+\frac{1}{2}\mathbf{u}_{\boldsymbol{\Sigma}}\,\mathbf{x}$
(19) $\displaystyle\frac{d\mathbf{u}_{\boldsymbol{\Sigma}}}{dt}$
$\displaystyle=\frac{1}{2}\boldsymbol{\Sigma}\left(\mathbf{B}-\mathbf{x}\boldsymbol{\mu}^{T}\right)\boldsymbol{\Sigma}^{-1}\mathbf{u}_{\boldsymbol{\Sigma}}+\frac{1}{2}\mathbf{u}_{\boldsymbol{\Sigma}}\left(\mathbf{B}-\mathbf{x}\boldsymbol{\mu}^{T}\right)-\frac{1}{2}\boldsymbol{\Sigma}\,\mathbf{x}\,\mathbf{u}_{\boldsymbol{\mu}}^{T}-\frac{1}{2}\mathbf{u}_{\boldsymbol{\mu}}\,\mathbf{x}^{T}\,\boldsymbol{\Sigma}.$

Integrating equations (19) backwards along a geodesic (that is, from $t=1$ to
$t=0$), we can thus shift the residual vector back to the start of the
geodesic, where it can be used to update our initial velocity guess.
Integration of these equations is achieved numerically, with the values of
$\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$ at each timepoint given by (8).

The issue is, for a non-Euclidean manifold the traditional concept of a
residual vector is unavailable. Han and Park [23] use the Euclidean residual
as a “good enough” choice in their shooting approach, reducing the magnitude
of resulting velocity updates when the magnitude of this vector grows too
large and the resulting approximation unsatisfactory. Observing that the
residual vector in Euclidean space is also the velocity of the geodesic
connecting the two points (a straight line), here we also consider the
approximate velocities in Section 3.2 as alternative choices of residual
vector (Section 4.2). Note, however, that even if the exact geodesic velocity
is available, this does not produce a perfect update.

In order to decide an appropriate step length for a velocity update, shooting
makes use of the Riemannian geometry concept known as the Jacobi field. The
Jacobi field describes the difference between infinitesimally close geodesics,
and hence in a fashion analogous to the derivative or Jacobian in other
contexts, can be used to obtain the linearised sensitivity of a geodesic to
some perturbation. To explicitly define the Jacobi field and the numerical
shooting algorithms, we first introduce the compact notation ${\cal
G}(t\mathbf{v})$ to denote the action of tracing out a geodesic through MVN
space from the origin $(\mathbf{0},\mathbf{I})$ with initial velocity
$\mathbf{v}=(\mathbf{x},\mathbf{B})$ out to distance $t$. This function ${\cal
G}(t\mathbf{v})$ is given by either of equations (8) or (9). The problem of
finding the geodesic between two points in MVN space may thus be expressed as
seeking an initial velocity $\mathbf{v}$ such that

$\boldsymbol{\theta}_{t}={\cal G}(\mathbf{v}),\qquad\qquad\mbox{or
equivalently,}\quad\mathbf{v}={\cal G}^{-1}(\boldsymbol{\theta}_{t}),$

where again
$\boldsymbol{\theta}_{t}=(\boldsymbol{\mu}_{t},\boldsymbol{\Sigma}_{t})$ is
the target point after shifting to the origin. We note that in differential
geometry, ${\cal G}(\mathbf{v})$ and ${\cal G}^{-1}(\boldsymbol{\theta})$ are
often described as the exponential map,
$\exp_{(\mathbf{0},\mathbf{I})}(\mathbf{v})$, and the logarithmic map,
$\log_{(\mathbf{0},\mathbf{I})}(\boldsymbol{\theta})$, respectively.

Using this notation, the Jacobi field associated with a given velocity update
$\Delta\mathbf{v}$ is given by

${\cal J}(t)=\left.\frac{d}{d\tau}{\cal
G}\Bigl{(}t(\mathbf{v}+\tau\Delta\mathbf{v})\Bigr{)}\right|_{\tau=0}=\lim_{\epsilon\rightarrow
0}\frac{{\cal G}\Bigl{(}t(\mathbf{v}+\epsilon\Delta\mathbf{v})\Bigr{)}-{\cal
G}(t\mathbf{v})}{\epsilon}.$ (20)

Estimating ${\cal J}(1)$ using finite differencing, we obtain the linearised
effect of a velocity update, $\Delta\mathbf{v}$, on the end point of the
geodesic,

${\cal J}(1)\approx\frac{{\cal G}(\mathbf{v}+s\Delta\mathbf{v})-{\cal
G}(\mathbf{v})}{s}.$ (21)

Seeking a scalar stepsize $s$ such that ${\cal G}(\mathbf{v}+s\Delta
v)\approx\boldsymbol{\theta}_{t}$, we can substitute this relation into
equation (21) and then scalarise it by taking an inner product of both sides,

$s\approx\frac{\Bigl{\langle}\boldsymbol{\theta}_{t}-{\cal
G}(\mathbf{v}),{\cal J}(1)\Bigr{\rangle}_{{\cal
G}(\mathbf{v})}}{\Bigl{\langle}{\cal J}(1),{\cal J}(1)\Bigr{\rangle}_{{\cal
G}(\mathbf{v})}}.$ (22)

Equation (22) mirrors the choice of Han and Park [23] to use the inner product
with the Jacobi field (projecting both sides of the equation onto one side).
The inner product used is the one that induces the metric (2), evaluated at
the geodesic endpoint. Just like the metric, this inner product depends also
on the location on the manifold $\boldsymbol{\theta}$, and for two vectors
$\mathbf{u}=(\mathbf{u}_{\boldsymbol{\mu}},\mathbf{u}_{\boldsymbol{\Sigma}})$
and
$\mathbf{u}^{\prime}=(\mathbf{u}^{\prime}_{\boldsymbol{\mu}},\mathbf{u}^{\prime}_{\boldsymbol{\Sigma}})$
takes the form

$\Bigl{\langle}\mathbf{u},\mathbf{u}^{\prime}\Bigr{\rangle}_{\mathbf{p}}=\mathbf{u}_{\boldsymbol{\mu}}^{T}\boldsymbol{\theta}_{\Sigma}^{-1}\mathbf{u}^{\prime}_{\boldsymbol{\mu}}+\frac{1}{2}\operatorname{tr}\Bigl{(}\boldsymbol{\theta}_{\boldsymbol{\Sigma}}^{-1}\mathbf{u}_{\boldsymbol{\Sigma}}\boldsymbol{\theta}_{\boldsymbol{\Sigma}}^{-1}\mathbf{u}^{\prime}_{\boldsymbol{\Sigma}}\Bigr{)}.$
(23)

With all of these pieces in place, shooting can proceed by repeatedly
iterating the process of choosing a vector connecting the current geodesic
endpoint to the actual target, parallel transporting this vector to the origin
via backwards integration of equations (19), and then choosing a stepsize for
this velocity update via equation (22) to obtain a new geodesic. Termination
can be triggered using the norm of residual vectors, or the symmeterised
Kullback–Leibler divergence as an easily-calculated proxy for squared manifold
distance,

$D_{\mathrm{sKL}}=\frac{1}{4}\Bigl{(}\operatorname{tr}\bigl{(}\boldsymbol{\Sigma}_{1}^{-1}\boldsymbol{\Sigma}_{0}+\boldsymbol{\Sigma}_{0}^{-1}\boldsymbol{\Sigma}_{1}-2\mathbf{I}\bigr{)}+(\boldsymbol{\mu}_{1}-\boldsymbol{\mu}_{0})^{T}\bigl{(}\boldsymbol{\Sigma}_{0}^{-1}+\boldsymbol{\Sigma}_{1}^{-1}\bigr{)}(\boldsymbol{\mu}_{1}-\boldsymbol{\mu}_{0})\Bigr{)}$

Two iterations of the algorithm are demonstrated graphically in Fig. 5, and
algorithm pseudocode is given in Algorithm S1 in the supplement.

Figure 5: The process of finding geodesics via shooting and parallel
transport, visualised for bivariate MVNs. Left: A geodesic is fired using an
initial velocity guess (red curve). A residual vector (magenta arrow) joining
the geodesic’s endpoint to the target (black) is found, and parallel
transported back along the curve to become a correction to its initial
velocity. The update is scaled using equation (22) to get the new velocity
vector (blue arrow) defining the improved geodesic (blue curve). Right: A
second shooting iteration improves the geodesic to the black one. Only the
$\boldsymbol{\mu}$-component of each vector is visualised, due to the nature
of bivariate MVN visualisation.

###  4.2 Alternate Choices of Residual

We find that approximating the geodesic velocity between a pair of points
produces significantly better residual vectors for shooting. In Fig. 6, we
demonstrate iteration counts for shooting across 5000 test problems generated
in the same fashion as in Section 3.1. As might be expected,
small-$\mathbf{x}$ approximations begin to hurt convergence as the distance to
the target increases and their performance during the initial shooting
iterations deteriorates. Defining residuals using the component-based or
projection-based approximations consistently and significantly improves mean
performance, the latter performing best.

a)

b)

Figure 6: Boxplots demonstrating the convergence rates of shooting for
different choices of residual vector over 5000 randomly-generated scenarios.
Lengths are binned to aid in summarising the results (numbers in parentheses
denote bin sample counts). a) Good approximations improve convergence,
especially as problem difficulty (distance to target) increases. The
projection-based approximation (18) performs best. b) Initialising shooting
using a geodesic approximation is generally beneficial, especially using
approximation (18).

Han and Park [23] initialise geodesic shooting using a zero velocity, a choice
that we find can help with convergence by avoiding pinning to a location that
generates bad shooting updates. However using the geodesic approximations we
have reviewed here, we may initialise shooting using a velocity that is
already close to the one being sought, potentially improving convergence and
reducing the chance of pinning. In Fig. 6b we compare between zero-initialised
and approximation-initialised shooting, for the two most competitive
approximations (component-wise, (17), and projection-based, (18)). Pinning was
not observed in these test cases, and initialisation generally speeds up
performance. However, as distance to the target grows larger and approximation
quality decreases, initialisation can become hurtful (observed here for
component-based initialisation). Initialisation using the projection-based
approximation remains beneficial for the target distances considered here
($d_{F}<8$).

###  4.3 Shooting via Path Refinement

Addressing the issues with convergence of shooting for pairs of MVNs separated
by larger Fisher–Rao distances, Han and Park [23] also proposed a method in
which a set of intermediate points between the start point and target are
refined until all lie along the desired geodesic path. This is achieved by
using the shooting described in Section 4 to find geodesic connections between
all of the even-numbered points and all of the odd-numbered points, and using
the midpoints of these geodesics to update the locations of the opposite set
of points. Iterating this process causes the path defined by these points to
settle onto the geodesic (see Fig. 7, and Algorithm S2 in the supplementary
material), and each individual set (odd or even) may be run in parallel. Once
the locations of the set of points have settled onto the geodesic, the stored
initial velocity for the geodesic between the first and third points can be
scaled according to the number of points to obtain the desired velocity. The
Fisher–Rao distance is given simply by totalling the lengths of each
individual geodesic segment.

Figure 7: An example iteration of path refinement, beginning with the
Euclidean path (grey) and attempting to find the geodesic path (black,
dashed). The geodesic between the even-index points (red) is used to position
the central odd-index points. These geodesics (blue), together with their
midpoints, provide an updated path that is closer to geodesic. Using more path
points (right panel) produces a smaller improvement in path quality per
iteration.

Path refinement breaks down a longer-distance shooting problem into many
shorter-distance shooting problems, avoiding the convergence issues with the
former. The algorithm requires a full application of the base shooting
algorithm (Algorithm S1) per point, per iteration, although we find this cost
can be significantly ameliorated by initialising each shooting call using the
geodesic found in the previous path refinement iteration. Instead, the issue
with path refinement is that as the number of points defining the path
increases, it becomes too easy for individual segments of the path to be
approximately geodesic, causing little update in their midpoints (Fig. 7)
while the overall path remains non-geodesic. This also requires us to specify
a stricter tolerance for the shooting subproblems than we desire for the
overall geodesic being sought.

These effects are demonstrated by using path refinement with a large number of
points ($N=129$) to find the geodesic connecting two distant MVNs,
$\left(\begin{bmatrix}1\\\ 2\end{bmatrix},\begin{bmatrix}1&0.1\\\
0.1&10\end{bmatrix}\right)$ and $\left(\begin{bmatrix}70\\\
35\end{bmatrix},\begin{bmatrix}10&-0.8\\\ -0.8&1\end{bmatrix}\right)$ with
$d_{F}=10.417$. Given their observed superiority in Section 4.2, in the
results that follow all calls to Algorithm S1 as part of the path refinement
process use projection-based residuals to formulate shooting updates. However,
we find this choice is far less important inside path refinement.

Fig. 8a demonstrates the failure of standard path refinement to converge
(dashed lines). After more than 8000 shooting iterations the endpoint of the
geodesic using the $\mathbf{v}$ returned by standard path refinement (with a
Euclidean initialisation) is still $D_{sKL}>10^{15}$ away from the target MVN,
highlighting the difficulty of estimating geodesic velocities over large
distances (for which small changes have great influence on the endpoint).
Performance in terms of Fisher–Rao distance is slightly better (Fig. 8b), but
shows a false convergence where subsequent iterations cease to improve an
incorrect length estimate. These effects are most pronounced when the number
of path points is high, however the failure to reach the tolerance used in
subproblems can be observed even for as few as five path points (Fig. S3).

a)

b)

Figure 8: Performance of standard (dashed lines) and adaptive (solid lines)
path refinement for different initialisations (colours), for a large number of
path points (N=129) and bivariate MVNs separated by $d_{F}=10.42$. a) Geodesic
velocity accuracy, in terms of distance $D_{\mathrm{sKL}}$ between geodesic
endpoint and target. Adaptive path refinement is required for convergence, and
initialising with a path closer to geodesic speeds convergence. b) Similar
performance occurs for Fisher–Rao distance calculation. The standard algorithm
can produce false convergence to an incorrect value.

Path refinement may be repaired by adaptively culling points as the algorithm
proceeds (Algorithm S3). Initialising path refinement with $2^{m}+1$ points
allows us to repeatedly delete all even-indexed points whenever the path
ceases to sufficiently improve, using the geodesics found using the higher
number of points as initialisations to help the shooting problems now between
more distant points to converge. Once the path is reduced to only three
points, a final run of basic shooting produces the geodesic, using an
approximately-geodesic path as a very strong initial guess. This restores the
ability of path refinement to converge to the desired geodesic (solid lines in
Fig. 8), and to directly specify the desired tolerance (Fig. S3).

Once using adaptive path refinement, we observe that the paths discussed in
Sections 3.1 and 3.2 significantly improve the algorithm’s performance. For
two means separated by a large distance, as in the test case here, the moment-
averaged path is closest to geodesic and provides the fastest convergence.
Projection-based paths (15) also offer a significant improvement in
convergence over the Euclidean initialisations used by Han and Park [23].
Although the best choice of path will be problem-specific, it is comparatively
cheap to form paths and calculate their lengths before performing any
shooting, informing the likely best choice of path to use in joining a given
pair of distant MVNs.

##  5 Conclusions

Although they are among the most mathematically tractable of distributional
families, MVNs have continued to defy a complete information geometric
treatment. Here we have presented a more pragmatic understanding of how the
paths of shortest distance along this manifold behave, and how this may inform
efforts to determine such a path between two given MVNs. Finding these paths,
or their length (the Fisher–Rao distance), is the primary way in which such
information geometry is practically applied [9, 23, 24, 25]. Despite the non-
trivial Riemannian geometry, we have demonstrated how these geodesics permit a
more intuitive, traditional geometric understanding that explains when
geodesics do or not travel straight lines in $\boldsymbol{\mu}$-space (and
hence, when a closed-form solution is available). This geometric understanding
highlights the importance of repeated eigenvalues of the transformed target
covariance matrix
$\boldsymbol{\Sigma}_{0}^{-1/2}\boldsymbol{\Sigma}_{1}\boldsymbol{\Sigma}_{0}^{-1/2}$,
which correspond to more expansive subspaces of points in
$\boldsymbol{\mu}$-space that can be reached via such geodesics.

Treating the distance along the statistical manifold associated with MVNs as
an information distance also provides a powerful means of interpreting
commonly-used paths in sequential or tempered sampling methods. Annealing
paths select intelligent, variance-informed paths along which to move the
location of a distribution, but fail to incorporate additional variance to
reduce the cost of this movement. Moment-averaged paths incorporate additional
variance to compensate for the naive path of movement for their location,
mimicking a feature of geodesic paths and hence explaining their generally
lower information cost. Wasserstein-optimal paths completely ignore the
informational context.

Without a closed-form solution available, geodesics on the MVN statistical
manifold are found via a numerical shooting process. Here we have thoroughly
reviewed this process and examined the concerns associated with its practical
usage. In particular, we have highlighted the use of geodesic velocity
approximations to accelerate shooting, as well as issues with convergence
inherent to the shooting via path refinement proposed by Han and Park [23] and
the use of path point culling to resolve them. Considering path information
costs also provides good initialisations for the path refinement process.

This work has concentrated on the information geometry of MVNs, acknowledging
their ubiquity in statistics and position as a comparatively tractable family
that still demonstrates the additional geometric complexities that come with
vector-valued random variables. However, many of the ideas we have summarised
and discussed here are extensible to other distributional families and
contexts. The intuition of temporarily increasing distribution scale in order
to lower the information cost associated with changes in location should be
applicable to many other distributions, and indeed (univariate) Cauchy and
Student’s $t$ distributions possess a Fisher geometry very similar to
univariate normal distributions [42]. The moment-averaged paths that also
capture this desirable property are well-defined (though not always in closed
form) for any exponential family distribution [12]. The MVN-manifold geodesic
connecting the Laplace approximations [43] for a pair of distributions might
also be used to inform information-efficient paths between the distributions
themselves.

A closed-form solution for the geodesic connecting two members of the MVN
family has continued to elude researchers in information geometry despite
sophisticated treatments of the topic via differential geometry and Lie theory
[20, 31]. By presenting here the information geometry of MVNs with only
minimal use of differential geometry, our hope is that these fascinating
concepts might become more broadly considered, and that this might increase
the opportunity for these remaining questions in the space of MVNs to be
solved. Until a closed-form solution for the geodesic (and Fisher–Rao
distance) between an arbitrary pair of MVNs is found, the numerical approaches
we have reviewed and extended here remain necessary.

## Acknowledgements

BAJL acknowledges infrastructural and research support from the QUT Centre for
Data Science. KB is supported by the ARC Centre of Excellence for Plant
Success in Nature and Agriculture (CE200100015). RWS is supported by the
Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES), Fundação
de Amparo à Pesquisa do Estado de Minas Gerais (FAPEMIG) and Conselho Nacional
de Desenvolvimento Científico e Tecnológico (CNPq).

## Code Availability

MATLAB implementations of routines for generating, finding and visualising
geodesics or other common paths through the family of multivariate normals,
including code used to generate the figures in this manuscript, is available
via
https://www.github.com/betalawson/MVNgeodesics.

## References

  * DasGupta [2008] Anirban DasGupta.  _Asymptotic Theory of Statistics and Probability_.  Springer, New York, 2008.
  * Bishop [2006] Christopher M. Bishop.  _Pattern Recognition and Machine Learning_.  Springer, 2006.
  * Li and Zhao [2023] Wuchen Li and Jiaxi Zhao.  Wasserstein information matrix.  _Inf. Geom._ , 2023.
  * Arjovsky et al. [2017] Martin Arjovsky, Soumith Chintala, and Léon bottou.  Wasserstein generative adversarial networks.  In Doina Precup and Yee Whye Teh, editors, _Proceedings of the 34th International Conference on Machine Learning_ , volume 70 of _Proceedings of Machine Learning Research_ , pages 214–223. PMLR, 06–11 Aug 2017.
  * Srivastava et al. [2018] Sanvesh Srivastava, Cheng Li, and David B. Dunson.  Scalable Bayes via barycenter in Wasserstein space.  _Journal of Machine Learning Research_ , 19(8):1–35, 2018.
  * Panaretos and Zemel [2019] Victor M. Panaretos and Yoav Zemel.  Statistical aspects of Wasserstein distances.  _Ann. Rev. Stat. Appl._ , 6:405–431, 2019.
  * Jeffreys [1946] Harold Jeffreys.  An invariant form for the prior probability in estimation problems.  _Proc. R. Soc. Lond. A_ , 186:453–461, 1946.
  * van Erven and Harremoës [2014] Tim van Erven and Peter Harremoës.  Rényi divergence and Kullback–Leibler divergence.  _IEEE Trans. Inf. Theory_ , 60:3797–3820, 2014.
  * Sim et al. [2012] Aaron Sim, Sarah Fillipe, and Michael P. H. Stumpf.  Information geometry and sequential Monte Carlo.  _arXiv_ , page arXiv:1212.0764, 2012.
  * Syed et al. [2021] Saifuddin Syed, Vittorio Romaniello, Trevor Campbell, and Alexandre Bouchard-Côté.  Parallel tempering on optimized paths.  _arXiv_ , page 2102.07720, 2021.
  * Gelman and Meng [1998] Andrew Gelman and Xiao-Li Meng.  Simulating normalizing constants: From importance sampling to bridge sampling to path sampling.  _Stat. Sci._ , 13:163–185, 1998.
  * Grosse et al. [2013] Roger B Grosse, Chris J Maddison, and Russ R Salakhutdinov.  Annealing between distributions by averaging moments.  In C.J. Burges, L. Bottou, M. Welling, Z. Ghahramani, and K.Q. Weinberger, editors, _Advances in Neural Information Processing Systems_ , volume 26. Curran Associates, Inc., 2013.  URL https://proceedings.neurips.cc/paper/2013/file/fb60d411a5c5b72b2e7d3527cfc84fd0-Paper.pdf.
  * Rao [1945] C. R. Rao.  Information and accuracy attainable in the estimation of statistical parameters.  _Bull. Calcutta Math. Soc._ , 37:81–91, 1945.
  * Landau and Lifshitz [1971] Lev D. Landau and Evgeny M. Lifshitz.  _The Classical Theory of Fields_.  Pergamon Press, Oxford, 1971.
  * Shun’ichi [1985] Amari Shun’ichi.  _Differential-Geometrical Methods in Statistics_.  Springer-Verlag, Berlin, 1985.
  * Nielsen [2020] Frank Nielsen.  An elementary introduction to information geometry.  _Entropy_ , 22:1100, 2020.
  * Banerjee et al. [2005] Arindam Banerjee, Srujana Merugu, Inderjit S. Dhillon, and Joydeep Ghosh.  Clustering with Bregman divergences.  _J. Mach. Learn. Res._ , 6:1705–1749, 2005.
  * Ikeda [2006] Kazushi Ikeda.  On geometric structure of quasi-additive learning algorithms.  In _The 2006 IEEE International Joint Conference on Neural Network Proceedsings_ , pages 1123–1127, 2006.
  * Martens [2020] James Martens.  New insights and perspectives on the natural gradient method.  _J. Mach. Learn. Res._ , 21:1–76, 2020.
  * Calvo and Oller [1991] Miquel Calvo and Josep M. Oller.  An explicit solution of information geodesic equations for the multivariate normal model.  _Statistics and Decisions_ , 9:119–138, 1991.
  * Transtrum et al. [2011] Mark K. Transtrum, Benjamin B. Machta, and James P. Sethna.  Geometry of nonlinear least squares with applications to sloppy models and optimization.  _Phys. Rev. E_ , 83:036701, 2011.
  * Sharp et al. [2022] Jesse A. Sharp, Alexander P. Browning, Kevin Burrage, and Matthew J. Simpson.  Parameter estimation and uncertainty quantification using information geometry.  _J. R. Soc. Interface_ , 19:20210940, 2022.
  * Han and Park [2014] Minyeon Han and F. C. Park.  DTI segmentation and fiber tracking using metrics on multivariate normal distributions.  _J. Math. Imaging Vis._ , 49:317–334, 2014.
  * Pilté and Barbaresco [2016] Marion Pilté and Frédéric Barbaresco.  Tracking quality monitoring based on information geometry and geodesic shooting.  In _2016 17th International Radar Symposium (IRS)_ , pages 1–6. IEEE, 2016.
  * Du et al. [2020] Kaixuan Du, Pin Wan, Yonghua Wang, Xiongzhi Ai, and Huang Chen.  Spectrum sensing method based on information geometry and deep neural network.  _Entropy_ , 22:94, 2020.
  * Skovgaard [1984] Lene T. Skovgaard.  A Riemannian geometry of the multivariate normal model.  _Scand. J. Stat._ , 11:211–223, 1984.
  * Nielsen [2023] Frank Nielsen.  A simple approximation method for the Fisher-–Rao distance between multivariate normal distributions.  _Entropy_ , 25:654, 2023.
  * Herntier and Peter [2022] Trevoer Herntier and Adrian M. Peter.  Transversality conditions for geodesics on the statistical manifold of multivariate Gaussian distributions.  _Entropy_ , 24:1698, 2022.
  * Malagò and Pistone [2015] Luigi Malagò and Giovanni Pistone.  Information geometry of the Gaussian distribution in view of stochastic optimization.  In _Proceedings of the 2015 ACM Conference on Foundations of Genetic Algorithms XIII_ , FOGA ’15, pages 150–162. Association for Computing Machinery, 2015.
  * Mahalanobis [1936] Prasanta C. Mahalanobis.  On the generalised distance in statistics.  _Proc. Natl. Inst. Sci. India_ , 2:49–55, 1936.
  * Imai et al. [2011] Takuro Imai, Akira Takaesu, and Masato Wakayama.  Remarks on geodesics for multivariate normal models.  _J. Math-for-Ind._ , 3:125––130, 2011.
  * Calin and Udrişte [2014] Ovidiu Calin and Constantin Udrişte.  _Geometric Modeling in Probability and Statistics_.  Springer Cham, Switzerland, 2014.
  * Pinele et al. [2020] Julianna Pinele, João E. Strapasson, and Sueli I. R. Costa.  The Fisher–Rao distance between multivariate normal distributions: Special cases, bounds and applications.  _Entropy_ , 22:404, 2020.
  * Kessy et al. [2018] Agnan Kessy, Alex Lewin, and Korbinian Strimmer.  Optimal whitening and decorrelation.  _The American Statistician_ , 72(4):309–314, 2018\.
  * Eriksen [1987] P. S. Eriksen.  Geodesics connected with the Fisher metric on the multivariate normal manifold.  In _Proceedings of the GST Workshop, Lancaster_ , pages 225–229, 1987\.
  * Strapasson et al. [2016] João E. Strapasson, Julianna Pinele, and Sueli I. R. Costa.  A totally geodesic submanifold of the multivariate normal distributions and bounds for the Fisher–Rao distance.  In _2016 IEEE Information Theory Workshop (ITW)_ , pages 61–65, 2016\.
  * Hukushima and Nemoto [1996] Koji Hukushima and Koji Nemoto.  Exchange Monte Carlo method and application to spin glass simulations.  _J. Phys. Soc. Japan_ , 65:1604–1608, 1996.
  * Duan and Fulop [2015] Jin-Chuan Duan and Andras Fulop.  Density-tempered marginalized sequential Monte Carlo samplers.  _J. Bus. Econ. Stat._ , 33:192–202, 2015.
  * Cuturi and Doucet [2014] Marco Cuturi and Arnaud Doucet.  Fast computation of Wasserstein barycenters.  In Eric P. Xing and Tony Jebara, editors, _Proceedings of the 31st International Conference on Machine Learning_ , volume 32 of _Proceedings of Machine Learning Research_ , pages 685–693, Bejing, China, 2014. PMLR.
  * Delon and Desolneux [2020] Julie Delon and Agnès Desolneux.  A Wasserstein-type distance in the space of Gaussian mixture models.  _SIAM J. Imaging Sci._ , 13:936–970, 2020.
  * Khan and Zhang [2022] Gabriel Khan and Jun Zhang.  When optimal transport meets information geometry.  _Inf. Geom._ , 5:47–78, 2022.
  * Mitchell [1988] Ann F. S. Mitchell.  Statistical manifolds of univariate elliptic distributions.  _Int. Stat. Rev._ , 56:1–16, 1988.
  * Tierney and Kadane [1984] Luke Tierney and Joseph B. Kadane.  Accurate approximatinos for posterior moments and marginal densities.  _J. Am. Stat. Assoc._ , 81:82–86, 1984.

## Supplementary Material

###  S.1 Derivation of the “Axis-aligned” Solution

For a geodesic beginning at the origin, the geodesic equations for
multivariate normals take the form (8). Using the affine transformation to
shift the geodesic between two points to one beginning at the origin, the
target point is given by equation (5), or in terms of canonical co-ordinates
for MVNs,

$\boldsymbol{\delta}_{t}=\mathbf{R}^{T}\boldsymbol{\Sigma}_{0}^{-1/2}\boldsymbol{\Sigma}_{1}^{-1}\left(\boldsymbol{\mu}_{1}-\boldsymbol{\mu}_{0}\right),\qquad\qquad\boldsymbol{\Delta}_{t}=\mathbf{R}^{T}\boldsymbol{\Sigma}_{0}^{-1/2}\boldsymbol{\Sigma}_{1}\boldsymbol{\Sigma}_{0}^{-1/2}\mathbf{R}.$

The remaining problem for finding the geodesic is to identify the initial
velocity parameters $\mathbf{x}$ and $\mathbf{B}$, such that at $t=1$ the
geodesic passes through the target point,

$\displaystyle\boldsymbol{\delta}_{t}$
$\displaystyle=-\mathbf{B}\Bigl{(}\cosh(\mathbf{G})-\mathbf{I}\Bigr{)}(\mathbf{G}^{-})^{2}\mathbf{x}+\sinh(\mathbf{G})\mathbf{G}^{-}\mathbf{x},$
$\displaystyle\boldsymbol{\Delta}_{t}$
$\displaystyle=\mathbf{I}+\frac{1}{2}\Bigl{(}\cosh(\mathbf{G})-\mathbf{I}\Bigr{)}+\frac{1}{2}\mathbf{B}\Bigl{(}\cosh(\mathbf{G})-\mathbf{I}\Bigr{)}(\mathbf{G}^{-})^{2}\mathbf{B}-\frac{1}{2}\sinh(\mathbf{G})\mathbf{G}^{-}\mathbf{B}-\frac{1}{2}\mathbf{B}\sinh(\mathbf{G})\mathbf{G}^{-}.$

As discussed in the main text, these solutions occur when $\mathbf{B}$ is a
diagonal matrix, and $\mathbf{x}$ is aligned with a single co-ordinate axis
(which also ensures $\mathbf{G}$ is diagonal). Considering the behaviour of
geodesics in equation (8) under these conditions, we see that
$\boldsymbol{\delta}_{t}$ remains aligned with the same axis as $\mathbf{x}$,
and $\boldsymbol{\Delta}$ remains a diagonal matrix. If the target point falls
within this space of valid targets, then we have a solution to the geodesic
connecting the origin to that target. In this case where $\mathbf{B}$ and
$\mathbf{G}$ are diagonal, the equations decouple in the sense that they may
be solved component-wise,

$\displaystyle\delta_{k}$
$\displaystyle=\left(-\frac{B_{kk}}{G_{kk}^{2}}\left(\cosh
G_{kk}-1\right)+\frac{1}{G_{kk}}\sinh G_{kk}\right)x_{k},$
$\displaystyle\Delta_{kk}$ $\displaystyle=1+\frac{1}{2}\left(\cosh
G_{kk}-1\right)+\frac{1}{2}\frac{B_{kk}^{2}}{G_{kk}^{2}}(\cosh
G_{kk}-1)-\frac{B_{kk}}{G_{kk}}\sinh G_{kk}.$

Each pair of equations can be related to the geodesics of the Poincaré half-
plane model [26]. For completeness, we present here the direct approach given
in [9] in a little extra detail. For notational simplicity, we drop the
component notation in what follows, such that non-bolded quantities refer to
the value in one component (diagonal component for matrices).

Defining $R=B/G$, the second equation can be rewritten

$\Delta=\frac{1}{2}\left(1-R^{2}\right)+\frac{1}{2}\left(1+R^{2}\right)\cosh
G-R\sinh G.$

Substitution of the exponential definitions for $\cosh$ and $\sinh$ gives

$\displaystyle\Delta$
$\displaystyle=\frac{1}{2}\left(1-R^{2}\right)+\frac{1}{4}\left(1+R^{2}\right)\left(e^{G}+e^{-G}\right)-\frac{1}{2}R\left(e^{G}-e^{-G}\right)$
$\displaystyle=\frac{1}{2}\left(1-R^{2}\right)+\left(\frac{1}{4}+\frac{1}{4}R^{2}-\frac{1}{2}R\right)e^{G}+\left(\frac{1}{4}+\frac{1}{4}R^{2}+\frac{1}{2}R\right)e^{-G}$
$\displaystyle=\frac{1}{2}\left(1-R^{2}\right)+\frac{1}{4}\Bigl{(}1-R\Bigr{)}^{2}e^{G}+\frac{1}{4}\Bigl{(}1+R\Bigr{)}^{2}e^{-G}$
$\displaystyle=\frac{1}{4}\Bigl{(}(1-R)e^{\frac{1}{2}G}+(1+R)e^{-\frac{1}{2}G}\Bigr{)}^{2}.$

Given that $R^{2}=\dfrac{B^{2}}{G^{2}}=\dfrac{B^{2}}{B^{2}+2x^{2}}\leq 1$, the
right hand side equation is a positive term squared, and as such we can take
the positive branch when taking the square root of both sides. This then gives
a quadratic in $e^{\frac{1}{2}G}$, with solution

$\displaystyle e^{\frac{1}{2}G}$
$\displaystyle=\frac{\sqrt{\Delta}\pm\sqrt{\Delta-1+R^{2}}}{(1-R)}.$

Given that $R\leq 1$, both branches give positive values and are potentially
valid. Re-arranging for $G$ gives

$G=\log\left(\frac{\sqrt{\Delta}\pm\sqrt{\Delta-1+R^{2}}}{1-R}\right)^{2}.$

Now returning to the first equation of the component-wise solution and making
the same substitution for $R$ including $x=\pm G\sqrt{\frac{1}{2}(1-R^{2})}$,
we have

$\displaystyle\delta$ $\displaystyle=\Bigl{(}-\frac{B}{G^{2}}(\cosh
G-1)+\frac{1}{G}\sinh G\Bigr{)}x$ $\displaystyle=\pm\Bigl{(}-R\cosh G+R+\sinh
G\Bigr{)}\sqrt{\frac{1}{2}(1-R^{2})}$

Setting
$\nu=\left(\frac{\sqrt{\Delta}\pm\sqrt{\Delta-1+R^{2}}}{(1-R)}\right)^{2}$
such that $G=\log(\nu)$ for notational convenience, we have that

$\displaystyle\delta$
$\displaystyle=\pm\Biggl{(}-R\left(\frac{\nu}{2}+\frac{1}{2\nu}\right)+R+\frac{\nu}{2}-\frac{1}{2\nu}\Biggr{)}\sqrt{\frac{1}{2}(1-R^{2})}$
$\displaystyle\delta^{2}$
$\displaystyle=\frac{1}{2}\Biggl{(}-R\left(\frac{\nu}{2}+\frac{1}{2\nu}\right)+R+\frac{\nu}{2}-\frac{1}{2\nu}\Biggr{)}^{2}\left(1-R^{2}\right)$
$\displaystyle=\frac{1}{8\nu^{2}}\left(-R\nu^{2}-R+2R\nu+\nu^{2}-1\right)^{2}\left(1-R^{2}\right)$
$\displaystyle=\frac{1}{8\nu^{2}}(\nu-1)^{2}\Bigl{(}\nu(1-R)+R+1\Bigr{)}^{2}(1-R^{2}).$

Substituting the definition of $\nu$, and after some simplification,

$\displaystyle\delta^{2}$
$\displaystyle=\frac{\Bigl{(}\Delta\pm\sqrt{\Delta}\sqrt{\Delta+1-R^{2}}-(1-R)\Bigr{)}^{2}\Bigl{(}\Delta\pm\sqrt{\Delta}\sqrt{\Delta+1-R^{2}}\Bigr{)}^{2}}{2\Bigl{(}\Delta\pm\sqrt{\Delta}\sqrt{\Delta-1+R^{2}}+\frac{1}{2}\bigl{(}1-R^{2}\bigr{)}\Bigr{)}^{2}}\frac{1+R}{1-R}.$

From here, computer algebra gives the solution

$R=\frac{\pm(\delta^{2}+2\Delta-2\Delta^{2})}{\sqrt{4\Delta^{4}+4\Delta^{2}\delta^{2}+\delta^{4}-8\Delta^{3}+4\Delta\delta^{2}+4\Delta^{2}}}.$

Note that this is true regardless of whether the positive or negative branch
in the definition of $\nu$ is selected.

As $G$ is defined only in terms of $G^{2}$ (consider the Taylor series
representation of the equations for $\delta$ and $\Delta$), we select $G$ to
always be positive. Then, considering that in the scenario that the mean does
not need to move ($\delta=0$), a decrease in covariance ($\Delta>1$) must
correspond to a negative value for $B$, we know to choose the positive branch.
As such, using the notation-simplifying variables from Table 1 in the main
document, we write

$R=\frac{\alpha}{\gamma},$

and comparing the expansion of $\alpha^{2}$ to the terms in $\gamma^{2}$
reveals the relation

$\gamma=\sqrt{\alpha^{2}+8\Delta^{2}\delta^{2}}.$

Substituting the definition of $R$ into the definition for $\nu$, and using
the relation
$G=\log(\nu)=\operatorname{acosh}\left(1+\dfrac{(\nu-1)^{2}}{2\nu}\right)$, it
can be shown that for one of the two branches of $\nu$,

$G=\operatorname{acosh}\left(1+\frac{\gamma^{2}}{8\Delta^{3}}\right),$

which matches the form for $g$ given in Table 1. Then, with $G$ found, we have
that

$B=RG=\frac{\alpha}{\gamma}G,$

and for $x$,

$\displaystyle x$ $\displaystyle=\pm
G\sqrt{\frac{1}{2}\Biggl{(}1-\left(\frac{\alpha}{\gamma}\right)^{2}\Biggr{)}}$
$\displaystyle=\pm G\sqrt{\frac{1}{2}\cdot
8\frac{\Delta^{2}\delta^{2}}{\gamma^{2}}}$
$\displaystyle=\frac{2G\delta\Delta}{\gamma}.$

The positive branch is selected so that the sign of $x$ matches the sign of
$\delta$, such that $\mu$ moves in the correct direction.

In the above, component notation was dropped for simplicity, and the above
procedure must be solved for each component separately. However, this solution
is only valid when $\mathbf{B}$ and $\mathbf{G}$ are diagonal matrices,
corresponding to the case where $\boldsymbol{\delta}=\delta\mathbf{e}_{k}$ is
aligned with a single co-ordinate axis. As such, only one component takes the
full form elucidated above. In the remainder of cases, $\delta=0$, which gives
$x=0$ and thus $B=G$ and $R=1$. The second equation for the form of the
geodesics then reduces simply to

$\displaystyle\Delta$ $\displaystyle=\cosh B-\sinh B$ $\displaystyle=e^{-B}$
$\displaystyle\log(\Delta^{-1})$ $\displaystyle=B$
$\displaystyle\log(\Sigma_{t})$ $\displaystyle=B,$

where $\Sigma_{t}$ is the target covariance for the current component. Using
this velocity in the geodesic equation gives

$\Sigma(t)=e^{t\log\Sigma_{t}}=\left(\Sigma_{t}\right)^{t},$

showing that the covariance transforms exponentially to its target value,
starting from the origin $\Sigma=1$ as we are in the transformed space. This
corresponds to the observation in the main text that the covariance
eigenvalues transform exponentially in directions other than the direction of
movement — that result is obtained simply by removing the diagonalising part
of the affine transform.

Under these conditions where $\mathbf{B}$ and $\mathbf{G}$ are diagonal and
the geodesic equations can be considered component wise, we obtain our
solutions for $\mathbf{x}$ and $\mathbf{B}$ by now adding these components,

$\mathbf{x}=\frac{2G\delta\Delta}{\gamma}\mathbf{e}_{k},\qquad\qquad\mathbf{B}=\frac{\alpha}{\gamma}G\mathbf{e}_{k}\mathbf{e}_{k}^{T}-\sum_{j\neq
k}\log D_{jj}\mathbf{e}_{j}\mathbf{e}_{j}^{T},$

where the index $k$ is defined by the condition
$\boldsymbol{\delta}=\delta\mathbf{e}_{k}$, that is, the component in which
the mean must move to reach the target.

###  S.2 A Small-$\mathbf{x}$ Velocity Approximation via Taylor Series

When $\mathbf{x}$ is small,
$\mathbf{G}^{2}=\mathbf{B}^{2}+2\mathbf{x}\mathbf{x}^{T}\approx\mathbf{B}^{2}$
and equations (8) for the form of the geodesics simplify considerably. The
second equation becomes

$\boldsymbol{\Delta}(t)\approx\mathbf{I}+\frac{1}{2}\Bigl{(}\cosh(t\mathbf{B})-\mathbf{I}\Bigr{)}+\frac{1}{2}\Bigl{(}\cosh(t\mathbf{B})-\mathbf{I}\Bigr{)}-\frac{1}{2}\sinh(t\mathbf{B})-\frac{1}{2}\sinh(t\mathbf{B})=e^{-t\mathbf{B}},$

and treating this as an equality and choosing $\mathbf{B}$ to satisfy
$\boldsymbol{\Delta}(1)=\boldsymbol{\Delta}_{t}$ gives the approximation
$\mathbf{B}\approx-\log(\boldsymbol{\Delta}_{t})=\log(\boldsymbol{\Sigma}_{t})$.
Expanding the Taylor series for $\cosh$ and $\sinh$ in the first of equations
(8), one obtains

$\boldsymbol{\delta}(t)=-\mathbf{B}\sum_{j=1}^{\infty}\frac{1}{2j!}(t\mathbf{G})^{2j-2}\mathbf{x}+\sum_{j=1}^{\infty}\frac{1}{(2j-1)!}(t\mathbf{G})^{2j-2}\mathbf{x},$

which after substituting in $t=1$ and again applying
$\mathbf{G}^{2}\approx\mathbf{B}^{2}$ gives

$\displaystyle\delta_{t}$
$\displaystyle\approx\Bigl{(}-\sum_{j=1}^{\infty}\frac{1}{2j!}\mathbf{B}^{2j-1}+\sum_{j=1}^{\infty}\frac{1}{(2j-1)!}\mathbf{B}^{2j-2}\Bigr{)}\mathbf{x}$
$\displaystyle\approx\mathbf{B}^{-1}\Bigl{(}-\sum_{j=1}^{\infty}(-1)^{j}\frac{1}{j!}\mathbf{B}^{j}\Bigr{)}\mathbf{x}$
$\displaystyle\approx\mathbf{B}^{-1}\Bigl{(}\mathbf{I}-e^{-\mathbf{B}}\Bigr{)}\mathbf{x}.$

Substituting the approximation for $\mathbf{B}$,
$\mathbf{B}\approx-\log\boldsymbol{\Delta}_{t}$, this gives for $\mathbf{x}$

$\mathbf{x}\approx\Bigl{(}\mathbf{I}-\boldsymbol{\Delta}\Bigr{)}^{-1}-\log(\boldsymbol{\Delta}_{t})\boldsymbol{\delta}_{t}.$

Expressing these approximate choices for $\mathbf{x}$ and $\mathbf{B}$ in
terms of the mean and covariance parameters of the target multivariate normal
distribution, they become

$\displaystyle\mathbf{x}$
$\displaystyle\approx\Bigl{(}\mathbf{I}-(\mathbf{I}-\boldsymbol{\Sigma}_{t})^{-1}\Bigr{)}\log(\boldsymbol{\Sigma}_{t})\boldsymbol{\Sigma}^{-1}_{t}\boldsymbol{\mu}_{t}$
$\displaystyle\mathbf{B}$ $\displaystyle\approx\log\boldsymbol{\Sigma}_{t},$

where the Sherman–Morrison–Woodbury formula is used to write

$\Bigl{(}\mathbf{I}-\mathbf{A}^{-1}\Bigr{)}^{-1}=\mathbf{I}-(\mathbf{I}-\mathbf{A})^{-1}.$

###  S.3 A Velocity Approximation via Higher-dimensional Manifold Embedding

The approximately geodesic paths found by embedding the MVN statistical
manifold in the space of symmetric positive definite matrices also produce a
velocity approximation, simply by taking the initial velocity of these paths.
As discussed in the main document, these paths take the form

$\mathbf{p}(t)=\begin{bmatrix}\boldsymbol{\Sigma}_{t}+\boldsymbol{\mu}_{t}\boldsymbol{\mu}_{t}^{T}&\boldsymbol{\mu}_{t}\\\
\boldsymbol{\mu}_{t}^{T}&1\end{bmatrix}^{t}.$

We first seek $\mathbf{p}^{\prime}(0)$, which is found by taking the logarithm
of both sides and differentiating,

$\frac{d}{dt}\log\mathbf{p}(t)=\log\begin{bmatrix}\boldsymbol{\Sigma}_{t}+\boldsymbol{\mu}_{t}\boldsymbol{\mu}_{t}^{T}&\boldsymbol{\mu}_{t}\\\
\boldsymbol{\mu}_{t}^{T}&1\end{bmatrix}$

Then, as for the origin-transformed problem we have
$\mathbf{p}(0)=\mathbf{I}$, $\mathbf{p}(0)$ and $\mathbf{p}^{\prime}(0)$
commute. As such we may use
$\frac{d}{dt}\left(\log\mathbf{p}\right)=\mathbf{p}^{-1}\frac{d\mathbf{p}}{dt}$
to obtain simply

$\mathbf{p}^{\prime}(0)=\log\begin{bmatrix}\boldsymbol{\Sigma}_{t}+\boldsymbol{\mu}_{t}\boldsymbol{\mu}_{t}^{T}&\boldsymbol{\mu}_{t}\\\
\boldsymbol{\mu}_{t}^{T}&1\end{bmatrix}.$

Now we take the derivative of the embedding of an MVN path
$\boldsymbol{\theta}(t)=\bigl{(}\boldsymbol{\mu}(t),\boldsymbol{\Sigma}(t)\bigr{)}$
into the same space,

$\displaystyle\mathbf{p}_{{\cal N}}^{\prime}(t)$
$\displaystyle=\frac{d}{dt}\begin{bmatrix}\boldsymbol{\Sigma}(t)+\boldsymbol{\mu}(t)\bigl{[}\boldsymbol{\mu}(t)\bigr{]}^{T}&\boldsymbol{\mu}(t)\\\
\bigl{[}\boldsymbol{\mu}(t)\bigr{]}^{T}&1\end{bmatrix}$
$\displaystyle=\begin{bmatrix}\boldsymbol{\Sigma}^{\prime}(t)+\boldsymbol{\mu}^{\prime}(t)\bigl{[}\boldsymbol{\mu}(t)\bigr{]}^{T}+\boldsymbol{\mu}(t)\bigl{[}\boldsymbol{\mu}^{\prime}(t)\bigr{]}^{T}&\boldsymbol{\mu}^{\prime}(t)\\\
\bigl{[}\boldsymbol{\mu}^{\prime}(t)\bigr{]}^{T}&0\end{bmatrix}.$

At $t=0$, $\boldsymbol{\mu}(0)=\mathbf{0}$,
$\boldsymbol{\Sigma}(0)=\mathbf{I}$ and the initial velocity defines
$\boldsymbol{\mu}^{\prime}(0)=\mathbf{x}$ and
$\boldsymbol{\Sigma}^{\prime}(0)=\mathbf{B}$,

$\mathbf{p}^{\prime}_{{\cal N}}(0)=\begin{bmatrix}\mathbf{B}&\mathbf{x}\\\
\mathbf{x}^{T}&0\end{bmatrix}.$

Equating the two then provides the approximate choice of $\mathbf{x}$ and
$\mathbf{B}$,

$\begin{bmatrix}\mathbf{B}&\mathbf{x}\\\
\mathbf{x}^{T}&0\end{bmatrix}\approx\log\begin{bmatrix}\boldsymbol{\Sigma}_{t}+\boldsymbol{\mu}_{t}\boldsymbol{\mu}_{t}^{T}&\boldsymbol{\mu}_{t}\\\
\boldsymbol{\mu}_{t}^{T}&1\end{bmatrix},$

where $\mathbf{x}$ and $\mathbf{B}$ can be read out from the appropriate
matrix blocks on the right-hand side. When
$\boldsymbol{\mu}_{t}=\mathbf{x}=\mathbf{0}$, this approximation reduces to
$\mathbf{B}=\log\boldsymbol{\Sigma}_{t}$ and hence recovers the “equal means”
solution exactly (see Table 1 in main document).

###  S.4 Algorithms

procedure MVNshooting($\boldsymbol{\theta},\boldsymbol{\theta}^{\prime}$,
$\mathbf{v}_{0}$, tol, r_norm_max)

$\mathbf{O}\leftarrow\left(\mathbf{0},\quad\mathbf{I}\right)$
$\triangleright$Define origin

$\boldsymbol{\theta}_{t}\leftarrow\left(\boldsymbol{\theta}_{\boldsymbol{\Sigma}}^{-1/2}(\boldsymbol{\theta}^{\prime}_{\boldsymbol{\mu}}-\boldsymbol{\theta}_{\boldsymbol{\mu}}),\quad\boldsymbol{\theta}_{\boldsymbol{\Sigma}}^{-1/2}\boldsymbol{\theta}^{\prime}_{\boldsymbol{\Sigma}}\boldsymbol{\theta}_{\boldsymbol{\Sigma}}^{-1/2}\right)$
$\triangleright$Target point after $\boldsymbol{\theta}$ is mapped to origin

$\mathbf{v}\leftarrow\mathbf{v}_{0}$ $\triangleright$Initialise velocity to
provided guess (can use zero)

$\mathbf{g}\leftarrow{\cal G}(\mathbf{v})$ $\triangleright$Endpoint
$\mathbf{g}$ for the initial geodesic via equation (9)

while symKL($\mathbf{g},\boldsymbol{\theta}_{t}$) $>$ tol do

$\mathbf{r}\leftarrow$ residualVector($\mathbf{g},\mathbf{p}_{t}$)
$\triangleright$Residual between geodesic endpoint and target

r_norm
$\leftarrow\Bigl{[}\Bigr{.}\textsc{innerProduct}(\mathbf{r},\mathbf{r},\mathbf{g})\Bigl{.}\Bigr{]}^{1/2}$
$\triangleright$Used for stepsize considerations

$\Delta\mathbf{v}\leftarrow$ backTransport($\mathbf{r},\mathbf{v}$)
$\triangleright$Convert residual to velocity correction direction

$h\leftarrow$ tol /
$\Bigl{[}\textsc{innerProduct}(\Delta\mathbf{v},\Delta\mathbf{v},\mathbf{O})\Bigr{]}^{1/2}$
$\triangleright$Select finite differencing stepsize

$\mathbf{j}\leftarrow\bigl{(}{\cal
G}(\mathbf{v}+h\Delta\mathbf{v})-\mathbf{g}\bigr{)}/h$
$\triangleright$Evaluate geodesic endpoint sensitivity as per (20)

$s\leftarrow$
innerProduct($\boldsymbol{\theta}_{t}-\mathbf{g},\mathbf{j},\mathbf{g}$) /
innerProduct($\mathbf{j},\mathbf{j},\mathbf{g}$) $\triangleright$Step length

$s\leftarrow$ min($s$,r_norm_max / r_norm $s$) $\triangleright$Optionally
avoid over-ambitious updates

$\mathbf{v}\leftarrow\mathbf{v}+s\Delta\mathbf{v}$ $\triangleright$Apply
velocity update

$\mathbf{g}\leftarrow{\cal G}(\mathbf{v})$ $\triangleright$Endpopint of
geodesic with updated velocity via (9)

end while

$\mathbf{v}\leftarrow\left(\boldsymbol{\theta}_{\boldsymbol{\Sigma}}^{1/2}\mathbf{v}_{\boldsymbol{\mu}},\quad\boldsymbol{\theta}_{\boldsymbol{\Sigma}}^{1/2}\mathbf{v}_{\boldsymbol{\Sigma}}\boldsymbol{\theta}_{\boldsymbol{\Sigma}}^{1/2}\right)$
$\triangleright$Transform velocity back into original space

return $\mathbf{v}$ $\triangleright$And/or return Fisher–Rao distance,
$\Bigl{[}\textsc{innerProduct}(\mathbf{v},\mathbf{v},\boldsymbol{\theta})\Bigr{]}^{1/2}$

end procedure

Algorithm S1 Procedure for finding geodesics between MVNs
$\boldsymbol{\theta}$ and $\boldsymbol{\theta}^{\prime}$ via shooting.

procedure
MVNpathRefinement($\boldsymbol{\theta},\boldsymbol{\theta}^{\prime}$, tol,
r_norm_max, N)

$(\mathbf{a}_{1},\mathbf{a}_{2},\ldots,\mathbf{a}_{n})\leftarrow$
closedFormPath($\boldsymbol{\theta},\boldsymbol{\theta}^{\prime},N$)
$\triangleright$$N$-point path between $\boldsymbol{\theta}$ and
$\boldsymbol{\theta}^{\prime}$

looping $\leftarrow$ true

while looping do

for $i\leftarrow 3$ by $2$ to $N-1$ do $\triangleright$Connect even-index pts.
to update odd-index pts.

$\mathbf{v}_{i-1}\leftarrow$
MVNshooting($\mathbf{a}_{i-1},\mathbf{a}_{i+1},\mathbf{v}_{i-1}$, tol,
r_norm_max) $\triangleright$Find geodesic

$\mathbf{a}_{i}\leftarrow{\cal G}_{\mathbf{a}_{i-1}}(0.5\mathbf{v}_{i-1})$
$\triangleright$Place point halfway along geodesic

$L_{i-1}\leftarrow\Bigl{[}\textsc{innerProduct}(0.5\mathbf{v}_{i-1},0.5\mathbf{v}_{i-1},\mathbf{a}_{i-1})\Bigr{]}^{1/2}$
$\triangleright$Store $d_{F}(\mathbf{a}_{i-1},\mathbf{a}_{i})$

end for

for $i\leftarrow 2$ by $2$ to $N-1$ do $\triangleright$Connect odd-index pts.
to update even-index pts.

$\mathbf{v}_{i-1}\leftarrow$
MVNshooting($\mathbf{a}_{i-1},\mathbf{a}_{i+1},\mathbf{v}_{i-1}$, tol,
r_norm_max) $\triangleright$Find geodesic

$\mathbf{a}_{i}\leftarrow{\cal G}_{\mathbf{a}_{i-1}}(0.5\mathbf{v}_{i-1})$
$\triangleright$Place point halfway along geodesic

$L_{i-1}\leftarrow\Bigl{[}\textsc{innerProduct}(0.5\mathbf{v}_{i-1},0.5\mathbf{v}_{i-1},\mathbf{a}_{i-1})\Bigr{]}^{1/2}$
$\triangleright$Store $d_{F}(\mathbf{a}_{i-1},\mathbf{a}_{i})$

end for

Ltot $\leftarrow$ sum($L_{i}$) $\triangleright$Calculate current manifold
distance of the path

$\mathbf{g}\leftarrow{\cal G}_{\boldsymbol{\theta}}(0.5(N-1)\mathbf{v}_{1})$
$\triangleright$Fire geodesic from first point out to target

if $\textsc{symKL}(\mathbf{g},\mathbf{q})<$ tol then $\triangleright$Terminate
loop if target point reached

looping $\leftarrow$ false

end if

end while

return $0.5(N-1)\mathbf{v}_{1}$ $\triangleright$And/or return Fisher–Rao
distance, Ltot

end procedure

Note the use here of notation ${\cal G}_{\boldsymbol{\theta}}(\mathbf{v})$
refers to tracing out the geodesic with initial velocity $\mathbf{v}$ starting
from point $\boldsymbol{\theta}$ (by mapping to the origin, using equation
(9), and then mapping back). The base notation ${\cal G}(\mathbf{v})$ is
equivalent to ${\cal G}_{(\mathbf{0},\mathbf{I})}(\mathbf{v})$.

Algorithm S2 The path refinement process for finding the geodesic between
$\boldsymbol{\theta}$ and $\boldsymbol{\theta}^{\prime}$

procedure
MVNadaptivePathRefinement($\boldsymbol{\theta},\boldsymbol{\theta}^{\prime}$,
tol, r_norm_max, $m$)

req_improvement $\leftarrow 0.01$ $\triangleright$Cull if an iteration fails
to reduce relative length by 1%

$N\leftarrow 2^{m+1}$ $\triangleright$Number of points for level $l$

$(\mathbf{a}_{1},\mathbf{a}_{2},\ldots,\mathbf{a}_{n})\leftarrow$
closedFormPath($\boldsymbol{\theta},\boldsymbol{\theta}^{\prime},N$)
$\triangleright$$N$-point path between $\boldsymbol{\theta}$ and
$\boldsymbol{\theta}^{\prime}$

looping $\leftarrow$ true

Ltot_old $\leftarrow\infty$

while looping do

for $i\leftarrow 3$ by $2$ to $N-1$ do $\triangleright$Connect even-index pts.
to update odd-index pts.

$\mathbf{v}_{i-1}\leftarrow$
MVNshooting($\mathbf{a}_{i-1},\mathbf{a}_{i+1},\mathbf{v}_{i-1}$, tol,
r_norm_max) $\triangleright$Find geodesic

$\mathbf{a}_{i}\leftarrow{\cal G}_{\mathbf{a}_{i-1}}(0.5\mathbf{v}_{i-1})$
$\triangleright$Place point halfway along geodesic

$L_{i-1}\leftarrow\Bigl{[}\textsc{innerProduct}(0.5\mathbf{v}_{i-1},0.5\mathbf{v}_{i-1},\mathbf{a}_{i-1})\Bigr{]}^{1/2}$
$\triangleright$Store $d_{F}(\mathbf{a}_{i-1},\mathbf{a}_{i})$

end for

for $i\leftarrow 2$ by $2$ to $N-1$ do $\triangleright$Connect odd-index pts.
to update even-index pts.

$\mathbf{v}_{i-1}\leftarrow$
MVNshooting($\mathbf{a}_{i-1},\mathbf{a}_{i+1},\mathbf{v}_{i-1}$, tol,
r_norm_max) $\triangleright$Find geodesic

$\mathbf{a}_{i}\leftarrow{\cal G}_{\mathbf{a}_{i-1}}(0.5\mathbf{v}_{i-1})$
$\triangleright$Place point halfway along geodesic

$L_{i-1}\leftarrow\Bigl{[}\textsc{innerProduct}(0.5\mathbf{v}_{i-1},0.5\mathbf{v}_{i-1},\mathbf{a}_{i-1})\Bigr{]}^{1/2}$
$\triangleright$Store $d_{F}(\mathbf{a}_{i-1},\mathbf{a}_{i})$

end for

Ltot $\leftarrow$ sum($L_{i}$) $\triangleright$Calculate current manifold
distance of the path

$\mathbf{g}\leftarrow{\cal G}_{\boldsymbol{\theta}}(0.5(N-1)\mathbf{v}_{1})$
$\triangleright$Fire geodesic from first point out to target

if $\textsc{symKL}(\mathbf{g},\mathbf{q})<$ tol then $\triangleright$Terminate
loop if target point reached

looping $\leftarrow$ false

end if

if Ltot / Ltot_old $>1-\texttt{req\\_improvement}$ & $N>3$ then

$N\leftarrow(N-1)/2+1$ $\triangleright$Delete even-index points

for $i\leftarrow 1$ to $N$ do

$\mathbf{a}_{i}\leftarrow\mathbf{a}_{2i-1}$ $\triangleright$Locations of
retained points

$\mathbf{v}_{i}\leftarrow 2\mathbf{v}_{2i-1}$ $\triangleright$Retained
geodesics (doubled as each goes twice as far)

end for

end if

end while

$\mathbf{v}\leftarrow\textsc{MVNshooting}(\boldsymbol{\theta},\boldsymbol{\theta}^{\prime},\mathbf{v}_{1},\texttt{tol},\texttt{r\\_norm\\_max})$
$\triangleright$Perform final shooting

return $\mathbf{v}$ $\triangleright$And/or return Fisher–Rao distance,
$\Bigl{[}\textsc{innerProduct}(\mathbf{v},\mathbf{v},\boldsymbol{\theta})\Bigr{]}^{1/2}$

end procedure

Note the use here of notation ${\cal G}_{\boldsymbol{\theta}}(\mathbf{v})$
refers to tracing out the geodesic with initial velocity $\mathbf{v}$ starting
from point $\boldsymbol{\theta}$ (by mapping to the origin, using equation
(9), and then mapping back). The base notation ${\cal G}(\mathbf{v})$ is
equivalent to ${\cal G}_{(\mathbf{0},\mathbf{I})}(\mathbf{v})$.

Algorithm S3 Adaptive path refinement for finding the geodesic between
$\boldsymbol{\theta}$ and $\boldsymbol{\theta}^{\prime}$

procedure innerProduct($\mathbf{u},\mathbf{u}^{\prime},\boldsymbol{\theta}$)
$\triangleright$Inner product between two vectors at point
$\boldsymbol{\theta}$

return
$\mathbf{u}_{\boldsymbol{\mu}}^{T}\boldsymbol{\theta}_{\boldsymbol{\Sigma}}^{-1}\mathbf{u}^{\prime}_{\boldsymbol{\mu}}+\frac{1}{2}\mathrm{tr}\left(\boldsymbol{\theta}_{\boldsymbol{\Sigma}}^{-1}\mathbf{u}_{\boldsymbol{\Sigma}}\boldsymbol{\theta}_{\boldsymbol{\Sigma}}^{-1}\mathbf{u}^{\prime}_{\boldsymbol{\Sigma}}\right)$

end procedure

procedure symKL($\boldsymbol{\theta},\boldsymbol{\theta}^{\prime}$)
$\triangleright$Evaluates symmeterised KL divergence between two MVNs

return
$\frac{1}{4}(\boldsymbol{\theta}^{\prime}_{\boldsymbol{\mu}}-\boldsymbol{\theta}_{\boldsymbol{\mu}})^{T}\left(\boldsymbol{\theta}_{\boldsymbol{\Sigma}}^{-1}+{\boldsymbol{\theta}^{\prime}}_{\boldsymbol{\Sigma}}^{-1}\right)(\boldsymbol{\theta}^{\prime}_{\boldsymbol{\mu}}-\boldsymbol{\theta}_{\boldsymbol{\mu}})+\frac{1}{4}\operatorname{tr}\left(\boldsymbol{\theta}_{\boldsymbol{\Sigma}}^{-1}\boldsymbol{\theta}^{\prime}_{\boldsymbol{\Sigma}}+{\boldsymbol{\theta}^{\prime}}_{\boldsymbol{\Sigma}}^{-1}\boldsymbol{\theta}_{\boldsymbol{\Sigma}}-2\mathbf{I}\right)$

end procedure

procedure residualVector($\boldsymbol{\theta},\boldsymbol{\theta}^{\prime}$)
$\triangleright$Returns a vector joining $\boldsymbol{\theta}$ to
$\boldsymbol{\theta}^{\prime}$ in some sense

$\boldsymbol{\theta}_{t}\leftarrow\left(\boldsymbol{\theta}_{\boldsymbol{\Sigma}}^{-1/2}(\boldsymbol{\theta}^{\prime}_{\boldsymbol{\mu}}-\boldsymbol{\theta}_{\boldsymbol{\mu}}),\quad\boldsymbol{\theta}_{\boldsymbol{\Sigma}}^{-1/2}\boldsymbol{\theta}^{\prime}_{\boldsymbol{\Sigma}}\boldsymbol{\theta}_{\boldsymbol{\Sigma}}^{-1/2}\right)$
$\triangleright$Target point after transforming to origin

method $\leftarrow$ {euclid, taylor, eigen, projection}

switch method do

case euclid

$(\mathbf{x},\mathbf{B})\leftarrow\Bigl{(}(\boldsymbol{\theta}_{t})_{\boldsymbol{\mu}},\quad(\boldsymbol{\theta}_{t})_{\boldsymbol{\Sigma}}-\mathbf{I}\Bigr{)}$
$\triangleright$Euclidean residual between origin and
$\boldsymbol{\theta}_{t}$

case taylor

$(\mathbf{x},\mathbf{B})\leftarrow$ Evaluate equation (16)

case eigen

$(\mathbf{x},\mathbf{B})\leftarrow$ Evaluate equation (17)

case projection

$(\mathbf{x},\mathbf{B})\leftarrow$ Match components in equation (18)

$\mathbf{r}\leftarrow\left(\boldsymbol{\theta}_{\boldsymbol{\Sigma}}^{1/2}\mathbf{x},\quad\boldsymbol{\theta}_{\boldsymbol{\Sigma}}^{1/2}\,\mathbf{B}\,\boldsymbol{\theta}_{\boldsymbol{\Sigma}}^{1/2}\right)$
$\triangleright$Transform the vector back into the original space

return $\mathbf{r}$

end procedure

Algorithm S4 Subfunctions used by the geodesic shooting routines

procedure backTransport($\mathbf{r},\mathbf{v}$) $\triangleright$Transports
$\mathbf{r}$ along geodesic defined by $\mathbf{v}$

steps_per_length $\leftarrow 250$ $\triangleright$Specifies timestepping
resolution

$L\leftarrow\Bigl{[}\textsc{innerProduct}(\mathbf{v},\mathbf{v},\mathbf{I})\Bigr{]}^{1/2}$
$\triangleright$Find length of geodesic

$N\leftarrow\textsc{ceil}(\texttt{steps\\_per\\_length}*L)$
$\triangleright$Select number of numerical integration points

$h\leftarrow-1/N$ $\triangleright$Timestep for geodesic numeric integration
(backwards)

$t=1$ $\triangleright$Initialise numerical integration at geodesic endpoint

for $i=1,\ldots,N$ do $\triangleright$Loop from $t=1$ back to $t=0$

$\mathbf{r}\leftarrow\textsc{RK4Step}(\textsc{parallelRHS}(\mathbf{r},t,\mathbf{v}),h,t)$
$\triangleright$RK4 for given RHS function

$t\leftarrow t+h$

end for

return $\mathbf{r}$

end procedure

procedure parallelRHS($\mathbf{u},t,\mathbf{v}$) $\triangleright$RHS of
equations (19) for updating $\mathbf{u}$

$\mathbf{g}\leftarrow{\cal G}(t\mathbf{v})$ $\triangleright$Find geodesic
endpoint at this time

$\frac{d\boldsymbol{\mu}}{dt}\leftarrow\mathbf{g}_{\boldsymbol{\Sigma}}\mathbf{v}_{\boldsymbol{\mu}}$
$\triangleright$Find geodesic rate of change from equation (7)

$\frac{d\boldsymbol{\Sigma}}{dt}\leftarrow\mathbf{g}_{\boldsymbol{\Sigma}}\left(\mathbf{v}_{\boldsymbol{\Sigma}}-\mathbf{v}_{\boldsymbol{\mu}}\mathbf{p}_{\boldsymbol{\mu}}^{T}\right)$
$\triangleright$Find geodesic rate of change from equation (7)

return
$\Biggl{(}\frac{1}{2}\Bigl{(}\frac{d\boldsymbol{\Sigma}}{dt}\mathbf{g}_{\boldsymbol{\Sigma}}^{-1}\mathbf{u}_{\boldsymbol{\mu}}+\mathbf{u}_{\boldsymbol{\Sigma}}\mathbf{g}_{\boldsymbol{\Sigma}}^{-1}\frac{d\boldsymbol{\mu}}{dt}\Bigr{)},\quad\frac{1}{2}\Bigl{(}\frac{d\boldsymbol{\Sigma}}{dt}\mathbf{g}_{\boldsymbol{\Sigma}}^{-1}\mathbf{u}_{\boldsymbol{\Sigma}}+\mathbf{u}_{\boldsymbol{\Sigma}}\mathbf{g}_{\boldsymbol{\Sigma}}^{-1}\frac{d\boldsymbol{\Sigma}}{dt}-\frac{d\boldsymbol{\mu}}{dt}\mathbf{u}_{\boldsymbol{\mu}}^{T}-\mathbf{u}_{\boldsymbol{\mu}}\left(\frac{d\boldsymbol{\mu}}{dt}\right)^{T}\Bigr{)}\Biggr{)}$

end procedure

Algorithm S5 Additional subfunctions used by the geodesic shooting routines

###  S.5 Supplementary Figures

Figure S1: Visualisation of some example geodesics for bivariate MVNs, and the
geodesics found by approximating the initial velocity to the target. The
component-based approximation is exact for straight-line geodesics (top
right), but the projection-based approximation is the best overall performer
(see also Fig. S2).

d = 2 d = 10

Figure S2: Performance (in terms of symmetric Kullback–Leibler divergence, see
main text) of approximate geodesic velocities connecting the origin to 5000
randomly-generated target MVNs ($\boldsymbol{\mu}$ components selected
$\mu_{j}\sim U[-5,5]$ and covariance eigenvalues $\log\lambda_{j}\sim
U[-5,5]$), for bivariate and 10-variate MVNs. Lengths are binned to aid in
summarising the results (numbers in parantheses denote bin sample counts).
Manifold projection performs best overall, and is noticeably more consistent
for test problems generated in this fashion. However in the lower dimension,
the component-wise approximation is sometimes superior for distant targets.
Figure S3: Standard path refinement (dashed lines) fails to find a geodesic to
the tolerance used for shooting subproblems, and suffers greatly from choosing
too many path points. Adaptive path refinement (solid lines) rectifies these
issues, but still benefits from choosing as few points as possible. Lines are
drawn to aid figure interpretation only, and simply connect the points of data
for each method.