spatial transformation $\Phi_A$ maps $A(t_i)$ to $I(t_0)$ and $\Phi_B$ maps $B(t_i)$ to $I_0$. Then $A(t_i)=I(t_0)\circ\Phi_A^{-1}$ and $B(t_i)=I(t_0)\circ\Phi_B^{-1}$. Furthermore, assume that $\Phi$ maps $A(t_i)$ to $B(t_i)$, i.e., $B(t_i) = A(t_i)\circ\Phi^{-1}$. Then $\Phi=\Phi_B\circ\Phi_A^{-1}$. Assuming that the geodesic between $I(t_0)$ and $A(t_i)$ is parameterized by the initial velocity $v^A$ and between $I(t_0)$ and $B(t_i)$ by the initial velocity $v^B$ and that we travel between $I(t_0)$ and $A(t_i)$ in time $t_i-t_0$ (and similarly for $B(t_i)$) we can rewrite the map between $A(t_i)$ and $B(t_i)$ based on the exponential map as $$\Phi = \mathrm{Exp_{Id}}((t_i-t_0)v^B)\circ  \mathrm{Exp_{Id}}(-(t_i-t_0)v^A),$$ which can be approximated to first order as $$\Phi \approx  \mathrm{Exp_{Id}}((t_i-t_0)(v^B-v^A)).$$ Hence, the squared geodesic distance between the two images can be approximated as $$d^2(A(t_i),B(t_i))\approx \frac{1}{2}(t_i-t_0)^2\langle K(m^B-m^A),m^B-m^A\rangle,$$ where $v^A = K m^A$ and $v^B = K m^B$. Hence, Eq.  becomes $$\begin{aligned}
& E(\overline{I}, \overline{m}) = \frac{1}{2}\langle \overline{m}, K\overline{m}\rangle \notag\\
& + \frac{1}{2\sigma^2}\sum_i (t_i-t_0)^2 \langle K(\overline{m}-m_i),\overline{m}-m_i \rangle,
\label{eq:approx}\end{aligned}$$ where $\overline{m}$ is the sought-for initial momentum of the regression geodesic and $m_i$ are the initial momenta corresponding to the geodesic connecting $\overline{I}$ (the starting image of the geodesic) and the measurements $Y_i$ in time $t_i-t_0$. Differentiating Eq.  w.r.t. $\overline{m}$ results in $$\begin{aligned}
\nabla_{\overline{m}} E = K [\overline{m} + \frac{1}{\sigma^2} \sum_i (t_i - 