metric in the Wasserstein space of gaussian measures

Let $\mathcal N_0^d$ be the space of positive definite matrices of size $d$. Also let $\mathcal{G}=GL(d,\mathbb R)$ as in the paper. We have the natural surjection $\Pi: \mathcal{G}\to \mathcal N_0^d$ defined by $\Pi(A)=A^TA$. The general linear group $\mathcal{G}$ has a Riemannian metric $G$ as an open submanifold of $\mathbb R^{d\times d}$; that is, the tangent space at every point of $\mathcal{G}$ is identified with $\mathbb R^{d\times d}$ which carries its standard metric. The author wants to "push forward" the metric $G$ by $\Pi$, which means putting a metric $g$ on $\mathcal N_0^d$ such that $\Pi$ becomes a Riemannian submersion. Of course, metrics do not normally push forward. The authors pretends that this issue does not exist and writes

We define a Riemannian metric $g$ on $\mathcal N_0^d$ by $$g_{\Pi}(d\Pi(Z),d\Pi(W))=G(Z_{\mathcal H},W_{\mathcal H}).$$

We are really defining the metric at $\Pi(A)$ for some $A\in \mathcal{G}$, and since $\Pi$ is not a bijection, one should explain why the choice of $A$ does not matter. Instead, $A$ is simply omitted from the formula.

The reason why $\Pi$ is able to push $G$ forward can be stated as: (1) the orthohonal group $O(d,\mathbb R)$ acts on $\mathcal{G}$ by left multiplication; (2) this action is isometric; (3) the fibers of $\Pi$ are orbits under this action.

On the next page there is an explicit formula (3.4) for $g$. Maybe it would be more logical to introduce $g$ by (3.4) and then check that it makes $\Pi$ a Riemannian submersion. Whatever. Regarding (3.4) it is important to emphasize (and the author does it) that the tangent vectors $X$ and $Y$ in (3.4) are not written in coordinates of the ambient space $\mathbb R^{d\times d}$. This is why your computation in the 2nd bullet item did not give the right result.

Let's analyze the case $d=1$. Here $\mathcal N_0^1=(0,\infty)$, covered by the map $x\to x^2$ which pushes forward the Euclidean metric on $\mathcal G=\mathbb R\setminus \{0\}$. The result is $\mathcal N_0^1 = ((0,\infty), \frac12 x^{-1/2}dx)$. Applied carelessly, the formula (3.4) would give a wrong result $g(dx,dx)=x\,dx^2$ instead of the correct one $g(dx,dx)=\frac{1}{4x}dx^2$. One has to transform the tangent vectors (multiplying them by $A^{-1}$ or something) before sticking them into $\operatorname{Tr}(XVY)$.

You also mentioned the square root map $r: \mathcal N_0^d\to \mathcal G$. This is a right inverse of $\Pi$ in the sense that $\Pi\circ r=\operatorname{id}_{\mathcal N_0^d}$. However, $r$ is not an isometric embedding. The reason is that the range of $r$ (which is the set of all positive definite matrices) is not horizontal with respect to $\Pi$: its tangent space is not the subspace on which $d\Pi$ is isometric. Here is a simple illustration in two dimensions: $(x,y)\mapsto x$ is a Riemannian submersion from $\mathbb R^2$ onto $\mathbb R$, the map $x\mapsto (x,e^x)$ is its right inverse, but it is not an isometric embedding of $\mathbb R$ into $\mathbb R^2$.

Here is an explicit calculation that confirms the statement in the preceding paragraph. Let $A=\begin{pmatrix}2 & 1 \\ 1 & 2 \end{pmatrix}$. The tangent space of $\mathcal{G}$ at $A$ is decomposed into the kernel of $d\Pi$, which is the span of $B=\begin{pmatrix}1 & 2 \\ -2 & -1 \end{pmatrix}$, and its orthogonal complement, on which $d\Pi$ is isometric. Since $B$ is not skew-symmetric, its orthogonal complement does not coincide with the space of symmetric matrices, i.e., with the tangent space to the range of $r$.

[Added] It's tempting to think that the metric on $\mathcal{N}_0^d$ should be $g_V(X,Y)=\frac{1}{4}\operatorname{Tr}(XV^{-1}Y)$, written without any coordinate tricks: $X$ and $Y$ are tangent vectors to $\mathcal{N}_0^d$ at $V$, which means they are symmetric matrices. But it does not seem to work out this way...

I'll try follow the idea of computations in the paper without necessarily agreeing with its details. Fix $V\in \mathcal{N}_0^d$ and let $A=V^{1/2}\in \mathcal{G}$. Given two vectors $X$ and $Y$ in $T_V \mathcal{G}$ (my notation for tangent space), we want to pull them back to horizontal vectors in $T_A \mathcal{G}$, and take the inner product there.

1. What are the vertical vectors $Z$ in $T_A \mathcal{G}$? They satisfy $(A+Z)^T (A+Z)=A^TA$ up to second order, which means $AZ+Z^TA=0$. Rewriting this as $AZ+(AZ)^T=0$, we see that $AZ$ is skew-symmetric.
2. What are the horizontal vectors $W$ in $T_A \mathcal{G}$? They are orthogonal to vertical vectors: $\operatorname{Tr}(Z^TW)=0$. Writing the latter as $\operatorname{Tr}((AZ)^T (A^{-1}W))=0$, we conclude that $A^{-1}W$ must be symmetric. Let's record it by saying that $W=A\widehat W$ where $\widehat W$ is symmetric.
3. The image of horizontal vector $A=A\widehat W$ under $d\Pi_A$ is $AW+W^TA=A^2\widehat{W}+\widehat{W}A^2 = V\widehat{W}+\widehat{W}V$. This is a symmetric matrix, of course.
4. Given a symmetric matrix $X$, we want to solve $V\widehat{W}+\widehat{W}V=X$ for unknown symmetric matrix $\widehat{W}$. Unable to so do explicitly (?), we resign to denoting the solution as $\widehat{W}_X$. We also have $\widehat{W}_Y$ for our other tangent vector.
5. So, $g_V(X,Y) = \operatorname{Tr} ((A\widehat{W}_X)^T(A\widehat{W}_Y)) = \operatorname{Tr}(\widehat{W}_X V \widehat{W}_Y)$.

The result matches the paper, except that the "coordinate transformation" was spelled out. I wonder if $\widehat{W}_X$ can be reasonably written in terms of $X$...