I think I might see what was confusing me. This is really a comment, but it's too long for the comment thread. As my example, let's take $$A(t) = \frac{1}{1+t^2} \begin{pmatrix} 2 & t \\ -t & -2 \end{pmatrix}$$ So we want to solve the differential equation $U'(t) = r A(t) U(t)$, where $U$ is a $2 \times 2$ matrix with initial condition $U(-1) = \mathrm{Id}$.

We can actually compute the eigenvalues of $A(t)$ explicitly: They are $\sqrt{4-t^2}/(1+t^2)$. We compute $\int_{-1}^1 \pm \sqrt{4-t^2}/(1+t^2) dt \approx \pm 3.03022$. So your formula, as I understand it, is $$U(1) = e^{3.03022 r} u_1 v_1^T + e^{-3.03022 r} u_2 v_2^T + \cdots$$ where $u_i$ and $v_i$ are the eigenvectors of $A(-1)$ and $A(1)$.

What I think was confusing me is that it is somewhat misleading to call this the leading terms. The later terms in the series look like $e^{3.03022 r} r^{-k} (\mbox{stuff})$, right? So they actually dominate the $e^{-3.03022 r}$ term.

I think I might see what was confusing me. This is really a comment, but it's too long for the comment thread. As my example, let's take $$A(t) = \frac{1}{1+t^2} \begin{pmatrix} 2 & t \\ -t & -2 \end{pmatrix}$$ So we want to solve the differential equation $U'(t) = r A(t) U(t)$, where $U$ is a $2 \times 2$ matrix with initial condition $U(-1) = \mathrm{Id}$.

We can actually compute the eigenvalues of $A(t)$ explicitly: They are $\sqrt{4-t^2}/(1+t^2)$. We compute $\int_{-1}^1 \pm \sqrt{4-t^2}/(1+t^2) dt \approx \pm 3.03022$. So your formula, as I understand it, is $$U(1) = e^{3.03022 r} u_1 v_1^T + e^{-3.03022 r} u_2 v_2^T + \cdots$$ where $u_i$ and $v_i$ are the eigenvectors of $A(-1)$ and $A(1)$.

What I think was confusing me is that it is somewhat misleading to call this the leading terms. The later terms in the series look like $e^{3.03022 r} r^{-k} (\mbox{stuff})$, right? So they actually dominate the $e^{-3.03022 r}$ term.

I wish I weren't having so much trouble getting good numerical data, it would probably clear up my confusion a lot. In the meantime, here is why I am worried.

Let $A(t)$, $B(t)$ and $C(t)$ be three $2 \times 2$ matrix-valued functions as above, with $A(1)=B(1)=C(1)$ (and hence the same at $-1$.) Let $X(r)$, $Y(r)$ and $Z(r)$ be the parallel transport from $-1$ to $1$ be the differential equations $\phi'(t) = r A(t) \phi(t)$, $\phi'(t) = r B(t) \phi(t)$ and $\phi'(t) = r C(t) \phi(t)$. As I understand it, your method gives asymptotic expansions $$X(r) \approx U \begin{pmatrix} e^{x_1 r} & 0 \\ 0 & e^{x_2 r} \end{pmatrix} V \quad Y(r) \approx U \begin{pmatrix} e^{y_1 r} & 0 \\ 0 & e^{y_2 r} \end{pmatrix} V \quad Z(r) \approx U \begin{pmatrix} e^{z_1 r} & 0 \\ 0 & e^{z_2 r} \end{pmatrix} V \quad (1)$$ where I have the SAME matrices $U$ and $V$ in each cases, because they depend only on the eigenvectors of $A(1)=B(1)=C(1)$ and of $A(-1)=B(-1)=C(-1)$.

Am I right about $(1)$?

If so, here is the issue. Look at the quadratic form $$\det(x X(r) + y Y(r) + z Z(r)) \approx \det(U) \left( e^{r x_1} x + e^{r y_1} y + e^{r z_1} z \right) \left( e^{r x_2} x + e^{r y_2} y + e^{r z_2} z \right) \det(V).$$

The matrix of this form has leading terms $$\begin{pmatrix} \exp(r(x_1+x_2)) & & \\ \exp(r\max(x_1+y_2, x_2+y_1)) & \exp(r(y_1+y_2)) & \\ \exp(r\max(x_1+z_2, x_2+z_1)) & \exp(r\max(y_1+z_2, y_2+z_1)) & \exp(r(z_1+z_2)) \\ \end{pmatrix}$$ as long as the approximations in $(1)$ are good enough that we don't get extra cross terms.

Unless I am very confused, I can construct $A(t)$, $B(t)$, $C(t)$ such that this quadratic form looks like $x^2+y^2+z^2 + (e^r+e^{-r}) (xy+xz+yz)$. And there are no real numbers $(x_1, x_2, y_1, y_2, z_1, z_2)$ with $x_1+x_2=y_1+y_2=z_1+z_2=0$ and $\max(x_1+y_2, x_2+y_1)=\max(x_1+z_2, x_2+z_1)=\max(y_1+z_2, y_2+z_1)=1$. So something is wrong...