This question has a solution presented in this paper even if with the jargon and notation of theoretical physics. So, I will use a somewhat different notation and I will change
$${\bf A}(t)\rightarrow -i{\bf A}(t).$$
Then, I will compute the eigenvalues and eigenvectors of ${\bf A}(t)$ through
$${A}(t)|n;t\rangle=\lambda_n(t)|n;t\rangle.$$
Now, you get a series with a leading order term
$${\bf B}(r)=\sum_n e^{i\gamma_n}e^{-ir\int_{-1}^1 dt\lambda_n(t)}|n;1\rangle\langle n;-1| \qquad r\rightarrow\infty$$
being $\gamma_n=\int_{-1}^1dt\langle n;t|i\partial_t|n;t\rangle$ known as geometric phase. Then, an expansion in the inverse of $r$ can be obtained with the matrix
$$\tilde {\bf A}(t)=-\sum_{n,m,n\ne m}e^{i(\gamma_n(t)-\gamma_m(t))}e^{-ir\int_{t_0}^tdt[\lambda_m(t)-\lambda_n(t)]}\langle m;t|i\partial_t|n;t\rangle|m;t_0\rangle\langle m;t_0|$$
being in this case
$$\tilde {\bf B}(r)=\prod_{-1}^1e^{-i\tilde {\bf A}(t)dt}$$
so that
$$B(r)=\sum_n e^{i\gamma_n}e^{-ir\int_{-1}^1 dt\lambda_n(t)}|n;1\rangle\langle n;-1|\tilde {\bf B}(r).$$
This represents a solution of the Schroedinger equation
$$-ir{\bf A}(t)B(r;t,t_0)=\partial_tB(r;t,t_0)$$
in the interval $t\in [-1,1]$ and $r\rightarrow\infty$.
An example:
$$ A(t) = \frac{1}{1+t^2} \begin{pmatrix} 2 & t\\ -t & -2 \end{pmatrix} $$
and one has to solve the problem $$ \dot U(t)=rA(t)U(t) $$ with $r\gg 1$. We want to apply the technique outlined above. We note that $A(t)$ is not Hermitian and so, solving the eigenvalue problem, we get $\lambda_{\pm}=\pm r\frac{\sqrt{4-t^2}}{1+t^2}$ and $$ v_+=\frac{1}{2}\begin{pmatrix} \sqrt{2+\sqrt{4-t^2}}\\ -\frac{t}{\sqrt{2+\sqrt{4-t^2}}}\end{pmatrix} \qquad v_-=\frac{1}{2}\begin{pmatrix}-\frac{t}{\sqrt{2+\sqrt{4-t^2}}} \\ \sqrt{2+\sqrt{4-t^2}}\end{pmatrix}. $$ But $v_+^Tv_-\ne 0$ and so these vectors are not orthogonal. We need to solve also the eigenvalue problem $u^T(A-\lambda I)=0$ producing the following eigenvectors $$ u_+=\frac{1}{2}\begin{pmatrix} \sqrt{2+\sqrt{4-t^2}}\\ \frac{t}{\sqrt{2+\sqrt{4-t^2}}}\end{pmatrix} \qquad u_-=\frac{1}{2}\begin{pmatrix} \frac{t}{\sqrt{2+\sqrt{4-t^2}}} \\ \sqrt{2+\sqrt{4-t^2}}\end{pmatrix}. $$ It is easy to see that $u_+^Tv_-=u_-^Tv_+=0$. It is important to note that $\lambda(t)=\lambda(-t)$ and $u_+(-t)=v_-(t)$ and $u_-(-t)=v_+(t)$ and so, these eigenvectors are just representing a backward evolution in time. Now, we want to study the time evolution of a generic eigenvector $$ \phi(t)=\begin{pmatrix}\phi_+(t) \\ \phi_-(t)\end{pmatrix} $$ and this can be done by putting $$ \phi(t)=c_+(t)e^{r\int_0^tdt'\frac{\sqrt{4-t^{'2}}}{1+t^{'2}}}v_+(t)+ c_-(t)e^{-r\int_0^tdt'\frac{\sqrt{4-t^{'2}}}{1+t^{'2}}}v_-(t) $$ that will produce the set of equations $$ \dot c_+=\gamma_+c_++e^{-2r\int_0^tdt'\frac{\sqrt{4-t^{'2}}}{1+t^{'2}}}\frac{u_+^T\frac{dv_-}{dt}}{u_+^Tv_+}c_- $$
$$ \dot c_-=\gamma_-c_-+e^{2r\int_0^tdt'\frac{\sqrt{4-t^{'2}}}{1+t^{'2}}}\frac{u_-^T\frac{dv_+}{dt}}{u_-^Tv_-}c_+ $$ having set $\gamma_+=\frac{u_+^T\frac{dv_+}{dt}}{u_+^Tv_+}$ and $\gamma_-=\frac{u_-^T\frac{dv_-}{dt}}{u_-^Tv_-}$. These equations are interesting because they provide the way time evolution is formed in a non-hermitian case. But this is also saying to us that each component may evolve in time differently: One can be really smaller than the other for $r\gg 1$. But we can also understand the form of the higher order corrections:
$$ c_+(t)=c_+(0)+\int_0^tdt'e^{\int_0^{t'}dt''(\gamma_+(t'')-\gamma_-(t''))}e^{-2r\int_0^{t'}dt''\frac{\sqrt{4-t^{''2}}}{1+t^{''2}}}\frac{u_+^T\frac{dv_-}{dt''}}{u_+^Tv_+}c_-(0)+\ldots. $$
Using a saddle point technique, we can uncover here that the correction is exponentially small and cannot be stated that is something like $e^{r}/r^k$ in the general case.
