Direct calculation of the Fisher distance via Riemannian geodesics

Question

I'm looking for a reference for a direct calculation of the Fisher distance (to avoid overloading the term "metric") $d_F(x,y) := 2 \cos^{-1} \sum_i \sqrt{x_i y_i}$ as the geodesic distance for the (Levi-Civita connection for the) Riemannian metric $g_{ij}(x) = \delta_{ij}/x_i$ restricted to the standard simplex.

I am aware of the usual trick wherein one considers $\sqrt{x_i}$ as coordinates and avoids even writing out the geodesic equation $\ddot{x}_i = -\sum_{j,k} \Gamma_{jk}^i \dot{x}_j \dot{x}_k$, much less solving it (for instance, I can read it in Ay et al.'s Information Geometry, or see it discussed here on MO). I could not care less about calculations/references that use tricks like this--I really and truly and only want to see a reference that goes through the gory details, if one exists.

I've been trying to do this calculation and it's eating me alive. I want to generalize it, but first I need to truly understand how to do it the ugly way.

EDIT: any future answers should elaborate/extend the partial answer I've been able to come up with below OR detail precisely how to numerically integrate the ODEs in some convenient language (say, MATLAB or [ugh] Python).

Answer 1

Consider the metric $g$ on the domain $$ Q = \{ (x^0, \dots, x^n)\ :\ x^0, x^1, \dots, x^n > 0\} $$ given by $$ g_{ij} = \frac{\delta_{ij}}{x^i},\ 0 \le i, j \le n. $$ Let $$ \Delta = \{ x^0 + \cdots x^n = 1 \} \subset Q. $$ On $\Delta$, $x^0 = 1-x^1-\cdots - x^n$. Therefore $$ dx^0 = -dx^1 - \cdots - dx^n, $$ and the metric $g$ restricted to $\Delta$ is \begin{align*} g &= \frac{(dx^1)^2}{x^1} + \cdots + \frac{(dx^n)^2}{x^n} + \frac{(dx^1 + \cdots + dx^n)^2}{x^0}\\ &= \left(\frac{\delta_{ij}}{x^i} + \frac{1}{x^0}\right)\,dx^i\,dx^j \end{align*} Therefore, the components of $g$ with respect to the coordinates $(x^1, \dots, x^n)$ are $$ g_{ij} = \frac{\delta_{ik}\delta_{jk}}{x^k} + \frac{1}{x^0},\ 1 \le i, j \le n. $$ Differentiating this, we get $$ \partial_kg_{ij} = -\frac{\delta_{ik}\delta_{jk}}{(x^k)^2} + \frac{1}{(x^0)^2}. $$ It follows that \begin{align*} 2g_{kp}\Gamma^p_{ij} &= \partial_ig_{kj} + \partial_jg_{ik} - \partial_kg_{ij}\\ &= -\frac{\delta_{ik}\delta_{ij}}{(x^i)^2} - \frac{\delta_{ji}\delta_{jk}}{(x^j)^2} + \frac{\delta_{ki}\delta_{kj}}{(x^k)^2} + \frac{1}{(x^0)^2}\\ &= -\frac{\delta_{ik}\delta_{jk}}{(x^k)^2} + \frac{1}{(x^0)^2}. \end{align*} Let $x: [0,1] \rightarrow \Delta$ be a constant speed geodesic. Then length $\ell$ of the geodesic is given by \begin{align*} \ell^2 &= g(\dot{x},\dot{x})\\ &= \frac{(\dot{x}^0)^2}{x^0} + \sum_k \frac{(\dot{x}^k)^2}{x^k}, \end{align*} and the geodesic equation is $$ 0 = \nabla_{\dot{x}}\dot{x} = (\ddot{x}^k + \dot{x}^i\dot{x}^j\Gamma^k_{ij})\partial_k. $$ Equivalently, \begin{align*} 0 &= g_{kp}(\ddot{x}^p + \dot{x}^i\dot{x}^j\Gamma^p_{ij})\\ &= \left(\frac{\delta_{kp}}{x^k} + \frac{1}{x^0}\right) \ddot{x}^p + \frac{1}{2}\dot{x}^i\dot{x}^j\left( -\frac{\delta_{ik}\delta_{jk}}{(x^k)^2} + \frac{1}{(x^0)^2}\right)\\ &= \frac{\ddot{x}^k}{x^k} + \frac{1}{x^0}\sum_p \ddot{x}^p + \frac{1}{2}\left(-\left(\frac{\dot{x}^k}{x^k}\right)^2 + \left(\frac{1}{x^0}\sum_i \dot{x}^i\right)^2\right)\\ &= \frac{\ddot{x}^k}{x^k} - \frac{\ddot{x}^0}{x^0} + \frac{1}{2}\left(-\left(\frac{\dot{x}^k}{x^k}\right)^2 + \left(\frac{\dot{x}^0}{x^0}\right)^2\right), \end{align*} which can be rewritten as \begin{align*} \ddot{x}^k &= \frac{1}{2}\frac{(\dot{x}^k)^2}{x^k} + x^k\left(\frac{\ddot{x}^0}{x^0} - \frac{1}{2}\left(\frac{\dot{x}^0}{x^0}\right)^2\right). \end{align*} Therefore, \begin{align*} 0 &= \ddot{x}^0 + \sum_k \ddot{x}^k\\ &= \ddot{x}^0 + \frac{1}{2}\sum_k \frac{(\dot{x}^k)^2}{x^k} + \left(\sum_k x^k\right)\left(\frac{\ddot{x}^0}{x^0} - \frac{1}{2}\left(\frac{\dot{x}^0}{x^0}\right)^2\right)\\ &= \ddot{x}^0 + \frac{1}{2}\sum_k \frac{(\dot{x}^k)^2}{x^k} + (1-x^0)\left(\frac{\ddot{x}^0}{x^0} - \frac{1}{2}\left(\frac{\dot{x}^0}{x^0}\right)^2\right)\\ &= \frac{\ddot{x}^0}{x^0} + \frac{1}{2}\left(\frac{(\dot{x}^0)^2}{x^0} + \sum_k \frac{(\dot{x}^k)^2}{x^k}\right) - \frac{1}{2}\left(\frac{\dot{x}^0}{x^0}\right)^2\\ &= \frac{\ddot{x}^0}{x^0} - \frac{1}{2}\left(\frac{\dot{x}^0}{x^0}\right)^2 + \frac{1}{2}g(\dot{x},\dot{x})\\ &= \frac{\ddot{x}^0}{x^0} - \frac{1}{2}\left(\frac{\dot{x}^0}{x^0}\right)^2 + \frac{\ell^2}{2}. \end{align*} This is a Riccati equation. If we let $u^0 = \sqrt{x^0}$, then \begin{align*} \dot{u}^0 &= \frac{1}{2}\frac{\dot{x}^0}{\sqrt{x^0}}\\ \ddot{u}^0 &= \frac{1}{2}\frac{\ddot{x}^0}{\sqrt{x^0}} - \frac{1}{4}\frac{(\dot{x}^0)^2}{x^0\sqrt{x^0}}\\ &= \frac{1}{2}\sqrt{x^0}\left(\frac{\ddot{x}^0}{x^0} - \frac{1}{2}\left(\frac{\dot{x}^0}{x^0}\right)^2\right)\\ &= -\frac{\ell^2}{4}u^0. \end{align*} It follows that $$ u^0(t) = a^0\cos \omega t + b^0\sin\omega t, $$ where $\ell = 2\omega$, and \begin{align*} x^0(t) &= (a^0\cos \omega t + b^0\sin\omega t)^2\\ &= (a^0)^2(\cos\omega t)^2 + 2a^0b^0 \cos\omega t \sin\omega t + (b^0)^2(\sin\omega t)^2. \end{align*}

THe argument above can be repeated with $x^0$ replaced by $x^k$, for each $1 \le k \le n$, yielding, for each $0 \le k \le n$, \begin{align*} u^k(t) &= a^k\cos\omega t + b^k\sin\omega t\\ x^k(t) &= (a^k\cos\omega t + b^k\sin\omega t)^2. \end{align*} If we denote $p = x(0)$ and $q = x(1)$, then, for each $0 \le k \le n$, \begin{align*} a^k &= \sqrt{p^k}\\ b^k &= \frac{\sqrt{q^k} - \sqrt{p^k}\cos\omega}{\sin\omega}. \end{align*} On the other hand, \begin{align*} 1 &= x^0 + \cdots + x^n\\ &= |a|^2(\cos\omega t)^2 + 2a\cdot b(\cos\omega t)(\sin\omega t) + |b|^2(\sin\omega t)^2, \end{align*} which implies that $|a|^2 = |b|^2 = 1$ and $a\cdot b = 0$. It follows that \begin{align*} 0 &= a\cdot b\\ &= \sum_{k=0}^n \frac{\sqrt{p^kq^k} - p^k\cos\omega}{\sin\omega}\\ &= \frac{\sum_{k=0}^n \sqrt{p^kq^k} - \cos\omega}{\sin\omega}. \end{align*} Therefore, given any two points $p, q \in \Delta$, the distance $d(p,q) = \ell = 2\omega$ satisfies $$ \cos\omega = \sum_{k=0}^n \sqrt{p^kq^k}, $$ and a straightforward calculation shows that the constant speed geodesic $x: [0,1] \rightarrow \Delta$ connecting them is given by $$ x(t) = \left(\frac{\sqrt{p^k}\sin((1-t)\omega) + \sqrt{q^k}\sin (t\omega)}{\sin\omega}\right)^2. $$

Answer 2

I was eventually able to get the answer (only) for $n = 1$, building on a calculation of the Christoffel symbols for arbitrary $n$. But this still seems far from a full answer.

Let $(\mu_1,\dots,\mu_n)$ parametrize the standard simplex $\Delta_{n} := \{(\mu_1,\dots,\mu_{n+1}) \in \mathbb{R}^{n+1}_{\ge 0}: \sum_i \mu_i = 1\}$, so that $\mu_{n+1} = 1 - \sum_{i = 1}^n \mu_i$. Henceforth, summations will be assumed to be over indices in $[n]$. This choice of notation for coordinates follows that in [Ay et al.]

In these coordinates, and recalling (2.21) of [Ay et al.], the Fisher-Rao metric is (in a form that makes symmetry manifest and thus also inhibits mistakes at the cost of lengthier expressions) \begin{equation} \label{eq:FiniteFisherRaoMetricInStandardProjection} g_{ij}(\mu) = \delta_{ij} \mu_i^{-1/2} \mu_j^{-1/2} + \mu_{n+1}^{-1}. \end{equation} Its inverse (recalling (2.22) of [Ay et al.]) is \begin{equation} \label{eq:InverseFiniteFisherRaoMetricInStandardProjection} g^{ij}(\mu) = \delta_{ij} \mu_i^{1/2} \mu_j^{1/2} - \mu_i \mu_j. \end{equation}

Recall that the Christoffel symbols for (the Levi-Civita connection for) $g_{ij}$ are \begin{equation} \Gamma_{ij}^k = \frac{1}{2} g^{k\ell} \left ( \partial_j g_{i\ell} + \partial_i g_{j\ell} - \partial_\ell g_{ij} \right ). \end{equation} Now since $\partial_\ell \mu_{n+1} = -1$, we have $\partial_\ell \mu_{n+1}^{-1} = \mu_{n+1}^{-2}$ and \begin{align} \label{eq:DerivativeOfFiniteFisherRaoMetricInStandardProjection} \partial_\ell g_{ij} & = \delta_{ij} \left ( \partial_\ell \left [\mu_i^{-1/2} \right ] \mu_j^{-1/2} + \mu_i^{-1/2} \partial_\ell \left [\mu_j^{-1/2} \right ] \right ) + \mu_{n+1}^{-2} \nonumber \\ & = \delta_{ij} \left ( -\frac{1}{2} \delta_{\ell i} \mu_i^{-3/2} \mu_j^{-1/2} - \frac{1}{2} \mu_i^{-1/2} \delta_{\ell j} \mu_j^{-3/2} \right ) + \mu_{n+1}^{-2} \nonumber \\ & = -\frac{1}{2} \delta_{ij} (\delta_{\ell i} + \delta_{\ell j}) \mu_i^{-1} \mu_j^{-1} + \mu_{n+1}^{-2} \nonumber \\ & = -\delta_{\ell i} \delta_{ij} \mu_i^{-1} \mu_j^{-1} + \mu_{n+1}^{-2}. \end{align} Permuting indices appropriately, we get the explicit form of the Christoffel symbols: \begin{align} \label{eq:ChristoffelOfFiniteFisherRaoMetricInStandardProjection} \Gamma_{ij}^k & = \frac{1}{2} \sum_\ell \left ( \delta_{k\ell} \mu_k^{1/2} \mu_\ell^{1/2} - \mu_k \mu_\ell \right ) \cdot \begin{pmatrix} -\delta_{ji} \delta_{i\ell} \mu_i^{-1} \mu_\ell^{-1} + \mu_{n+1}^{-2} \\ -\delta_{ij} \delta_{j\ell} \mu_j^{-1} \mu_i^{-1} + \mu_{n+1}^{-2} \\ +\delta_{\ell i} \delta_{ij} \mu_i^{-1} \mu_j^{-1} - \mu_{n+1}^{-2} \end{pmatrix} \nonumber \\ & = \frac{1}{2} \sum_\ell \left ( \delta_{k\ell} \mu_k^{1/2} \mu_\ell^{1/2} - \mu_k \mu_\ell \right ) \cdot \left ( -\delta_{ij} \delta_{j\ell} \mu_j^{-1} \mu_i^{-1} + \mu_{n+1}^{-2} \right ) \nonumber \\ & = \frac{1}{2} \sum_\ell \begin{pmatrix} -\delta_{k \ell} \delta_{j \ell} \delta_{ij} \mu_k^{1/2} \mu_\ell^{1/2} \mu_i^{-1} \mu_j^{-1} + \delta_{k \ell} \mu_k^{1/2} \mu_\ell^{1/2} \mu_{n+1}^{-2} \\ + \delta_{j \ell} \delta_{ij} \mu_i^{-1} \mu_j^{-1} \mu_k \mu_\ell - \mu_k \mu_\ell \mu_{n+1}^{-2} \end{pmatrix} \nonumber \\ & = \frac{1}{2} \left ( -\delta_{ij} \delta_{jk} \mu_k \mu_i^{-1} \mu_j^{-1} + \mu_k \mu_{n+1}^{-2} + \delta_{ij} \mu_i^{-1} \mu_k - \mu_k \left [ 1 - \mu_n+1 \right ] \mu_{n+1}^{-2} \right ) \nonumber \\ & = \frac{1}{2} \left ( -\delta_{ij} \delta_{jk} \mu_k \mu_i^{-1} \mu_j^{-1} + \delta_{ij} \mu_i^{-1} \mu_k + \mu_k \mu_{n+1}^{-1} \right ) \nonumber \\ & = \frac{\mu_k}{2} \left ( \delta_{ij} \mu_i^{-1} \left [ 1 - \delta_{jk} \mu_j^{-1} \right ] + \mu_{n+1}^{-1} \right ). \end{align}

Inserting this dog's breakfast into the geodesic equation yields a more symmetric mess: \begin{align} \label{eq:GeodesicOfFiniteFisherRaoMetricInStandardProjection} \ddot{\mu}_k & = -\sum_{ij} \Gamma_{ij}^k \dot{\mu}_i \dot{\mu}_j \nonumber \\ & = -\frac{\mu_k}{2} \sum_{ij} \left ( \delta_{ij} \mu_i^{-1} \left [ 1 - \delta_{jk} \mu_j^{-1} \right ] + \mu_{n+1}^{-1} \right ) \dot{\mu}_i \dot{\mu}_j \nonumber \\ & = -\frac{\mu_k}{2} \left ( \sum_{ij} \delta_{ij} \mu_i^{-1} \dot{\mu}_i \dot{\mu}_j - \sum_{ij} \delta_{ij} \delta_{jk} \mu_i^{-1} \mu_j^{-1} \dot{\mu}_i \dot{\mu}_j + \mu_{n+1}^{-1} \sum_{ij} \dot{\mu}_i \dot{\mu}_j \right ) \nonumber \\ & = -\frac{\mu_k}{2} \left ( \sum_i \mu_i^{-1} \dot{\mu}_i^2 - \sum_i \delta_{ik} \mu_i^{-2} \dot{\mu}_i^2 + \mu_{n+1}^{-1} \left [ \sum_i \dot{\mu}_i \right ]^2 \right ) \nonumber \\ & = -\frac{\mu_k}{2} \left ( \sum_i \mu_i^{-1} \dot{\mu}_i^2 - \mu_k^{-2} \dot{\mu}_k^2 + \mu_{n+1}^{-1} \left [ \sum_i \dot{\mu}_i \right ]^2 \right ). \end{align}

Before resorting to the elegant and "usual" change of variable $\lambda^2 := \mu$ that lets us work on the sphere and exploit the fact that geodesics are great circles without getting bogged down in the sorts of calculations above in the first place, we first execute a frontal assault on the simple case $n = 1$. Some of this ground is covered in section 2 of [Ciaglia et al.]. For now write $\mu \equiv \mu_1$, $g \equiv g_{11}$, and $\Gamma \equiv \Gamma_{11}^1$. We have $$g = \frac{1}{\mu} + \frac{1}{1-\mu} = \frac{1}{\mu(1-\mu)}$$ and $$\Gamma = \frac{\mu}{2} \left( \frac{1}{\mu} \left [1-\frac{1}{\mu} \right ] + \frac{1}{1-\mu} \right ) = \frac{2\mu-1}{2\mu(1-\mu)}.$$ Thus $$\ddot{\mu} = -\Gamma \dot{\mu}^2 = -\frac{2\mu-1}{2\mu(1-\mu)} \dot{\mu}^2.$$

The second order autonomous ODE $\ddot{\mu} = -\Gamma \dot{\mu}^2$ can be solved as follows. Writing $v:= \dot{\mu}$, we have that $\ddot{\mu} = \dot{v} = \frac{d\mu}{dt} \frac{dv}{d\mu} = v \frac{dv}{d\mu}$, so the ODE becomes $v \frac{dv}{d\mu} = -\Gamma v^2$. Thus $\int \frac{dv}{v} = -\int \Gamma \ d\mu$, so $v = \exp(-\int \Gamma \ d\mu)$. Meanwhile, $-\int \Gamma \ d\mu \equiv \log \sqrt{\mu(1-\mu)} + C_1$, so $v = \dot{\mu} = e^{C_1} \sqrt{\mu(1-\mu)}$. Integrating this in turn yields an equation for $t$: $$t = e^{-C_1} \int \frac{d\mu}{\sqrt{\mu(1-\mu)}} = -2e^{-C_1} \cos^{-1} \sqrt{\mu} + C_2.$$ Setting the integration constants $C_1, C_2 = 0$ yields $$t = -2 \cos^{-1} \sqrt{\mu}.$$

Let $0 \le \mu, \mu' \le 1$. Now $$\left | \cos^{-1} \sqrt{\mu} - \cos^{-1} \sqrt{\mu'} \right | = \cos^{-1} \left ( \sqrt{\mu}\sqrt{\mu'} + \sqrt{1-\mu} \sqrt{1-\mu'} \right )$$ because the expression on the left is the angle between $(\sqrt{\mu},\sqrt{1-\mu})^T$ and $(\sqrt{\mu'},\sqrt{1-\mu'})^T$, while the expression on the right is the angle between the same two vectors as computed using the Euclidean inner product.

This finally yields that for $n = 1$, the geodesic distance associated with the Fisher-Rao metric is indeed the Fisher-Rao distance. It is possible to play tricks to extend this observation directly to the case $n > 1$, but that would defeat my purpose and not teach me anything new, since I already know the general form of the Fisher-Rao distance. So I hope that any other answers will show how to proceed from the geodesic system of ODEs formulated above. On that point, I understand how the change of coordinates $\lambda := \sqrt{\mu}$ yields a system of ODEs of the form $\ddot{\lambda} = -\|\dot{\lambda}\|^2 \lambda$, which describes uniform circular motion and all that. However in the generalization I have in mind I don't expect to have this sort of step available to me.