Consider a triangle whose vertices are \(x, y, z\). In a Euclidean space, the vector that brings us from \(x\) to \(z\) is the sum of that which brings us from \(x\) to \(y\) and that which brings us from \(y\) to \(z\): \[ (z - x) = (y - x) + (z - y). \] Curvature breaks the analogous statement on a Riemannian manifold. However, we may still expect that the equality would hold approximately if \(x, y, z\) are near each other. That is indeed the case, as shown below.
Let \(\calM\) be a Riemannian manifold. Pick a closed ball \(\bar{B}(\bar{x}, r) = \{ x \in \calM : \dist(x, \bar{x}) \leq r \}\) in such a way that \(\inj(x) > 2r\) for all \(x\) in the ball. (Given any \(\bar{x}\), it is possible to pick a valid \(r > 0\), owing to continuity of the injectivity radius \(\inj\) as a function on \(\calM\).) In particular, \(\inj(\bar{x}) > r\) and hence the ball is compact. (We could simplify and summon a g-convex set.)
Fix an arbitrary point \(z \in \bar{B}(\bar{x}, r)\). Consider the following function based on the Riemannian distance: \[ f(x) = \frac{1}{2} \dist(x, z)^2. \] Its gradient is \(\grad f(x) = -\Log_x(z)\). Its Hessian is \(L\)-Lipschitz continuous in the ball \(\bar{B}(\bar{x}, r)\) (for some \(L < \infty\)) because that ball is compact and \(f\) is \(C^\infty\) on that ball.
Pick any two points \(x, y\) in \(\bar{B}(\bar{x}, r)\). Together with \(z\), they form a geodesic triangle, with sides determined by the unique minimizing geodesics joining them. Consider the vectors \[\begin{align*} u_x & = \Log_x(y), & v_x & = \Log_x(z), & w_y & = \Log_y(z). \end{align*}\] We are interested in the vector \[ \xi = v_x - u_x - P_{x \leftarrow y}\,(w_y), \] where parallel transport \(P\) along minimizing geodesics is here well defined. (In a Euclidean space, we would have \(\Log_x(y) = y-x\) and \(P_{x \leftarrow y}\,\) is identity, so that \(\xi = 0\).) Notice that \[\begin{align} \xi & = P_{x \leftarrow y}\,\grad f(y) - \grad f(x) - \Log_x(y) \\ & = \Big(P_{x \leftarrow y}\,\grad f(y) - \grad f(x) - \Hess f(x)[\Log_x(y)]\Big) \\ & \qquad + \Big(\Hess f(x)[\Log_x(y)] - \Log_x(y)\Big). \end{align}\] Since \(\Hess f\) is \(L\)-Lipschitz continuous in the ball that contains the geodesic connecting \(x\) and \(y\), the first parenthesized term has norm bounded by \(\frac{L}{2} \dist(x, y)^2\) (Boumal 2023, Prop. 10.55). Moreover, \(\Hess f(z) = I_z\) (identity on \(\T_z \calM\)). Again using that \(\Hess f\) is \(L\)-Lipschitz continuous (this time along the geodesic connecting \(x\) and \(z\)), it follows that \[\begin{align} \Hess f(x) - I_x & = \Hess f(x) - P_{x \leftarrow z} \circ I_z \circ P_{z \leftarrow x} \\ & = \Hess f(x) - P_{x \leftarrow z} \circ \Hess f(z) \circ P_{z \leftarrow x} \end{align}\] is bounded in operator norm by \(L \dist(x, z)\). Overall, we have: \[ \|\xi\| \leq L\left( \frac{1}{2} \dist(x, y) + \dist(x, z) \right) \dist(x, y). \] This holds for all \(x, y, z\) in the ball \(\bar{B}(\bar{x}, r)\), with the same constant \(L\).