Proving triangle inequality for an affine invariant distance on $\mathcal{Sym}_n^+$ In the paper A Riemannian framework for tensor computing, by Pennec et al., on page 46 the authors state a "distance" function on the manifold of  positive definite matrices $\mathcal{Sym}_n^+$ given by
$$d(A,B) = \lVert\log A^{-1/2}BA^{-1/2}\rVert$$
where the norm on the right is the standard sum of squares euclidean norm (presumably on the eigenvalues of the argument to the log-term).
The authors then say that they could not prove (but had strong empirical evidence) that the above function $d$ satisfies the triangle inequality. I am not that familiar with differential geometry, so I don't know where to check whether this property holds (I need it in an application).
Where should I look for a proof / disproof of the triangle inequality for $d$?
 A: It seems that $d(A,B)$ is the distance function for the (unique up to scalar multiplication) Riemannian metric on $\mathcal{Sym}_n^+$ invariant by the $\operatorname{GL}_n(R)$ action $(g,A)\mapsto gAg^t$. Indeed $d(A,B)$ is then equal to $d(I,A^{-1/2}BA^{-1/2})$, and if $S$ is any symmetric matrix, $t\mapsto \exp(tS)$ is a geodesic through $I$ with speed $\lVert S\rVert=\operatorname{tr}(S^2)^{1/2}$.
To see that it is indeed geodesic, note that the inverse map $s_I:B\mapsto B^{-1}$ is a riemannian isometry ("symmetry about $I$"), as well as $s_A:B\mapsto A^{1/2}B^{-1}A^{1/2}$ ("symmetry about $A$"). But then when $A$ is near enough $I$, the (unique) geodesic segment from $I$ to $A^2$ has to be invariant (as a set) by  $s_A$ (which switches its end points), so to contain $A$. By iterating and passing to the limit, the geodesic has to coincide with a one parameter group $t\mapsto\exp(tS)$ for small $t$, hence for all $t$.
ADDED: in particular, $d$ satisfies the triangle inequality. The Riemannian manifold $\mathcal{Sym}_n^+$ is the prototypical symmetric space, and most books on riemannian geometry must treat it, for instance Gallot–Hulin–Lafontaine - Riemannian geometry (I've not checked).
A: Let $X$, $Y$ and $Z$ be positive definite Hermitian matrices (you ask about the real symmetric case, the Hermitian one includes it). Let the eigenvalues of $(X^*)^{-1/2} Y X^{-1/2}$ be $e^{\alpha_i}$, those of $(Y^*)^{-1/2} Z Y^{-1/2}$ be $e^{\beta_i}$ and those of $(X^*)^{-1/2} Z X^{-1/2}$ be $e^{\gamma_i}$. Set $A=Y^{1/2} X^{-1/2}$, $B=Z^{1/2} Y^{-1/2}$ and $C=Z^{1/2} X^{-1/2}$. The singular values of $A$ are then $e^{\alpha_i/2}$, and so forth. 
Note that $AB=C$. By a result of Klyachko, there exist Hermitian matrices $\mathfrak{a}$, $\mathfrak{b}$ and $\mathfrak{c}$ such that $\mathfrak{a}+\mathfrak{b}=\mathfrak{c}$, the eigenvalues of $\mathfrak{a}$ are $\alpha_i/2$, those of $\mathfrak{b}$ are $\beta_i/2$ and those of $\mathfrak{c}$ are $\gamma_i/2$. This result can be thought of as saying that, although $\log A + \log B \neq \log C$, and although $\log A$, $\log B$ and $\log C$ are not Hermitian, we can find matrices which do have those properties and have the same eigenvalues. 
The inequality you want to prove is that 
$$\left( \sum \alpha_i^2 \right)^{1/2} + \left( \sum \beta_i^2 \right)^{1/2} \geq \left( \sum \gamma_i^2 \right)^{1/2}$$
But $\sum \alpha_i^2 = 4 \mathrm{Tr} \ \mathfrak{a}^* \mathfrak{a}$. So this turns into the (standard) fact that $\mathrm{Tr} \ \mathfrak{a}^* \mathfrak{a}$ is a positive definite norm on the Hermitian matrices. 

I can't resist indulging in a little self promotion. This argument appears at the beginning of my paper Horn's Problem, Vinnikov Curves and the Hive Cone. There I consider the curve $\det(xX+yY+zZ)=0$ in $\mathbb{RP}^2$. This curve is hyperbolic and, by results of Helton and Vinnikov, all hyperbolic curves are of this form. It meets the three coordinate lines precisely at $-e^{\alpha_i}$, $-e^{\beta_i}$ and $-e^{\gamma_i}$. The point is to use results from the theory of hyperbolic curves to explain Horn's results on eigenvalues of matrix sums.
A: You can understand this geometrically by interpreting a positive definite symmetric matrix as a Euclidean metric on $\mathbb R^n$.  Imagine starting with a lump of some kind of moldable material, like clay, and reshaping it (we may as well restrict to linear deformations) until its shape is that defined by a given symmetric matrix. (Actually, for this particular metaphor, we need to restrict to the case that the volume element $\det(A)$ remains constant, but the determinant is easily separated out in both the formula and the process.)
What is the total energy needed to change the shape? Here, total energy of an infinitesimal change is the $L^2$ norm of the change in length of unit vectors (up to a constant, this can be measured by integrating over the unit sphere, or summing over any orthonormal basis).
It's intuitively obvious and easy to prove that the most efficient path is to find the principal directions and stretch them log-linearly to get to the desired endshape.
To measure the distance between two metrics, you can do the same thing: find the principal directions of one with respect to the other.  That's what the formula does: $A^{-1/2}$ is
a linear transformation that sends the metric defined by $A$ to the standard Euclidean metric, and $A^{-1/2} B A^{-1/2}$ is the metric $B$ transformed into the new coordinates.
A: This distance is related to the notion of Geometric Mean in $\mathcal{Sym}_n^+$ (see exercises 198/199 of my extra exercises for Matrices : Theory and Applications). The geometric mean of $A$, $B$ is given by (assume that $A$ is positive definite)
$$A\mathbin\sharp B=A^{1/2}\left(A^{-1/2}BA^{-1/2}\right)^{1/2}A^{1/2}.$$
Although this is not clear on the formula, $A\mathbin\sharp B=B\mathbin\sharp A$. This mean turns out to be the middle point of the geodesic segment $[A,B]$ for the metric $d$ of the question. The geodesic segment is unique, the underlying Riemannian manifold is hyperbolic. The segment is parametrized by
$$s\mapsto A^{1/2}\left(A^{-1/2}BA^{-1/2}\right)^{1-s}A^{1/2}.$$
The following inequality relates the geometric, arithmetic and harmonic means
$$\frac12\left(A^{-1}+B^{-1}\right)^{-1}\le A\mathbin\sharp B\le\frac12(A+B).$$
Finally, the geometric mean of the harmonic and arithmetic means is the geometric mean. Actually, the arithmetico-harmonic mean (defined as a limit by iterating both arithmetic and harmonic means) is the geometric mean.
