Suppose $M$ is manifold (smooth, compact, without boundary if needed), and let $d(\cdot,\cdot)$ be the geodesic distance function. For a fixed $x_0\in M$, we can define $f(x):=d(x_0,x)$ to be the single-source distance function. Away from singularities, $f$ satisfies the eikonal equation $\|\nabla f\|_2\equiv 1.$
I'd like to recover $f$ using a convex optimization problem. Suppose I relax the eikonal condition to a convex condition $\|\nabla f\|_2\leq1$, and take $\mu$ to be an arbitrary measure supported on all of $M$.
Can I recover $f$ using the following optimization problem? $$ \left\{ \begin{array}{rl} \sup_{f\in C^\infty(M)} & \int f(x)\,d\mu(x)\\ \textrm{s.t.} & \|\nabla f(x)\|_2\leq1\ \forall x\in M\\ & f(x_0) = 0 \end{array}\right. $$ Empirically this appears to be the case via some numerical experiments, and it makes some sense since the viscosity solution of the eikonal equation roughly satisfies these constraints.
Any pointers to relevant theory or simple arguments are much appreciated!
Your claim is true, provided that $\mu$ is supported all over $M$.
If $\mu$ is only sparsely supported, then your claim about the optimal solution $f$ to your convex optimization problem might also solve the eikonal is not true in general. In that case $|\nabla f| = 1$ holds only along the Wasserstein-1 optimal transport path from the point source $x_0$ to the target measure $\mu$. In particular, if $\mu$ consists also of a few point measure, then away from the geodesics connecting $x_0$ to those points, one does not expect $|\nabla f| = 1$; instead one observes $|\nabla f|<1$. For example, if $\mu$ is a point measure, say $\mu=\delta_{x_1}$. Then modifying $f$ mildly (so that it doesn't violate $|\nabla f|\leq 1$) on the far side of $x_1$ from $x_0$ will not change the value of the objective function (due to the sparsity of $\mu$).
The relevant theory is the dual formulation to the Wasserstein-1/Beckmann problem: https://en.wikipedia.org/wiki/Wasserstein_metric#Dual_representation_of_W1
Your optimization problem can be reformulated as follows. Define $\tilde\mu:=\mu - \left(\int_Md\mu\right)\delta_{x_0}$, where $\delta_{x_0}$ is the Dirac-delta measure at $x_0$ with unit mass. So $\tilde\mu$ is a signed measure with zero mean. Consider the following problem (a dual Wasserstein-1 problem) $$ \begin{cases} \sup_{f\in C^\infty}\int_M f d\tilde\mu\\ \text{subject to }|\nabla f|\leq 1\text{ pointwise.} \end{cases}\qquad\mathrm{(*)} $$ In this new problem, there is one obvious degeneracy: Adding a constant to $f$ (i.e. applying the transformation $f\mapsto f+c$) does change the optimality or the value of the maximization (this uses the zero mean property of $\tilde\mu$). In particular, one may add an additional constraint $f(x_0) = 0$ to the problem ($\ast$), which will then recover your optimization problem. In short, your optimization (let's label it ($\dagger$)) is equivalent to ($\ast$): an optimal solution to ($\dagger$) is an optimal solution to ($\ast$), and conversely an optimal solution to ($\ast$) is an optimal solution to ($\dagger$) after a constant shift $f\mapsto f-f(x_0)$.
Now that your optimization is equivalent to ($\ast$) which is of the form of the dual Wasserstein-1 problem, we may apply what we know about the solutions to the dual Wasserstein-1 problem.
The solution $f$ of ($*$) satisfies the eikonal equation $|\nabla f| = 1$ only along the optimal transport paths. More precisely, consider the Beckmann problem (the dual problem of ($*$)) $$ \begin{cases} \inf_{X\in\Gamma(TM)}\int_M|X|\\ \text{subject to }\nabla\cdot X = \tilde \mu \end{cases}\qquad(\ddagger) $$ whose optimal solution will be a vector field $X$ concentrated in the paths connecting the negative part ($x_0$) of $\tilde \mu$ to the positive part (supports of $\mu$) of $\tilde\mu$. (For instance Fig. 2 of https://people.csail.mit.edu/jsolomon/assets/w1.pdf ) On $\operatorname{supp}(X)$ we have $|\nabla f| = 1$. Elsewhere $f$ can take arbitrary value as long as $|\nabla f|\leq 1$.
One can see this by recognizing that $f$ is the Lagrange multiplier for $(\ddagger)$. The Euler–Lagrange equation for $(\ddagger)$ is derived as follows. The subdifferential of the functional $E = \int_M|X|$ (note that it is not differentiable when $X=0$) is given by $$ \partial E = \left\{Y\in\Gamma(TM)\,\Big|\, Y(x)=X(x)/|X(x)|\text{ if $|X(x)|\neq 0$}, \text{ and }|Y(x)|\leq 1\text{ if $|X(x)|=0$}\right\}. $$ On the other hand, the functional gradient of the constraint paired with Lagrange multiplier is given by $$ {\delta\over\delta X} \int_M f\nabla\cdot X = -\nabla f $$ assuming the no-flux condition for $X$ on the boundary. Therefore the Euler–Lagrange equation for $(\ddagger)$ is given by $$ \begin{cases} \nabla f(x) = {X(x)\over |X(x)|},&\text{for } X(x)\neq 0\\ |\nabla f(x)|\leq 1,&\text{for all $x\in M$}\\ \nabla\cdot X = \tilde\mu. \end{cases} $$ Therefore, away from the (potentially very sparse) support of $X$, we do not necessarily get $|\nabla f|=1$.
So, as long as you make sure that the solution $X$ of the associated Beckmann problem ($\ddagger$) is almost everywhere non-vanishing, you have an eikonal equation solution $|\nabla f|=1$; moreover, $\nabla f = X/|X|$. For example, making the support of $\mu$ all over $M$ will fulfill such a condition.