I am interested in finding numerical solutions to a Monge-Ampere type equation for applications in physics. Due to the close connection between Monge-Ampere and optimal transport and the well developed libraries for solving the latter, I was thinking it would be great if I could use optimal transport solvers to get Monge-Ampere solutions.

The fly in the ointment is that numerical optimal transport libraries usually present their solutions in the form of an optimal transport plan. My question is how to infer the solution to a corresponding Monge-Ampere equation from an optimal transport plan. What follows is some background and notation to clarify the question.

My Monge-Ampere equation is of the form $$\det\nabla^2\phi(x) = \frac{f(x)}{g(\nabla \phi(x))}$$ where $f$ and $g$ are given probability density functions on the unit square $[0,1]^2$ and we seek solutions $\phi: [0,1]^2\rightarrow \mathbb{R}$ which are convex and satisfy the usual "second boundary condition". The gradient of the solution, $\nabla\phi$, yields a transport map from $f$ to $g$ which is optimal with respect to the usual quadratic cost. The corresponding transport plan $\gamma$ is the probability measure on $[0,1]^2 \times [0,1]^2$ with marginals $f$ and $g$ and supported on points of the form $(x,\nabla\phi(x))$.

Numerical libraries (at least the ones I know of) give you the plan $\gamma$ to some level of precision. Inferring $\phi$ from $\gamma$ can be logically broken down into two steps (though practically this may not be the best way to do it):

1. First you get the transport map $\nabla\phi$ from the plan $\gamma$.
2. Then you integrate $\nabla\phi$ to get $\phi$.

Because of finite precision of the numerical algorithms, neither of these steps is likely to have an exact solution. The ideally perfectly sharp support of $\gamma$ will get smeared out, leading to non-existence of a perfect transport map, and an approximate transport map will not in general be the gradient of any function.

I can think of naïve ways to implement both of these steps. To estimate the transport map from the pdf of the transport plan (which by abuse of notation I'll also call $\gamma$), one could use something like $\overline{\nabla \phi}(x) := max_y \gamma(x,y)$. To then get $\phi$ one could do a Helmholtz decomposition and keep just the scalar part.

Is there a better, more principled method for recovering $\phi$ from the corresponding transport plan? Any references are welcome. I am an outsider to both PDE and optimal transport.