My question is about to what extent can we extend the implicit function theorem to Riemannian manifolds. In the Euclidean space, consider a bivariate function $F \colon \Theta \times \mathcal{X} \rightarrow \mathbb{R}$ where $\Theta$ and $\mathcal{X}$ are Euclidean spaces. For any $\theta \in \Theta$, assume that $F(\theta, \cdot)$ is differentiable and strongly convex on $\mathcal{X}$ such that there exists a unique vector $x^\star(\theta)$ satisfying
$$ x^{\star }(\theta) = \arg\min_{x\in\mathcal{X} } F(\theta, x). $$ We are interested in finding $\nabla_{\theta} x^\star (\theta)$, i.e., how the global minimizer is sensitive with respect to the change in $\theta$. What we do is to differentiate the first-order condition. In particular, $x^{\star }(\cdot)$ satisfies $$ \nabla_{x} F(\theta, x) \big\vert_{x = x^\star (\theta)} = 0 \qquad \forall \theta\in \Theta. $$ Taking derivatives with respect to $\theta$ yields $$ \nabla^2 _{x \theta} F(\theta, x) \big\vert_{x = x^\star (\theta)} + \Big[ \nabla^2 _{x x} F(\theta, x) \big\vert_{x = x^\star (\theta)} \Big] \nabla_{\theta} x^\star (\theta) = 0. $$ Solving this equation yileds $\nabla_{\theta} x^\star (\theta)$.
Now suppose $\Theta \subseteq \mathbb{R}^d$, $\mathcal{X}$ is a Riemannian manifold, and $F(\theta, \cdot)$ is geodesically strongly convex on $\mathcal{X}$, can we generalize such reasoning using the Riemannian gradient and Hessian?
Furthermore, suppose now $\mathcal{X}$ is a Wasserstein space, i.e., $\mathcal{X}$ contains all probability densities defined on $\mathbb{R}^p$ and is equipped with the second-order Wasserstein distance $W_2$. And $F(\theta, \cdot)$ is a strongly convex functional with respect to $W_2$, how to compute $\nabla_{\theta} x^\star (\theta)$?
Finally, when $\Theta$ itself is a Riemannian manifold, can we generalize the above idea to compute $\mathrm{grad} (x^\star(\theta))$ where $\mathrm{grad}$ stands for the Riemannian gradient.
