Inverse of vector-valued function

Question

Having a differentiable function $f : \mathbb{R}^2 \longrightarrow \mathbb{R}^2$ such that $\det D(f)$ (Jacobi determinant) does not vanish. How can we get the inverse function $f^{-1} : \mathbb{R}^2 \longrightarrow \mathbb{R}^2$ such that it can be handled/plotted/differentiated like we do for $f$?

InverseFunction does not cover this case. NSolve, Root are unsuitable.

E.g.:

f[x_,y_]:={2x+Sin[x+y],2y+Cos[x+y]}
StreamDensityPlot[f[x,y],{x,-3,5},{y,-3,5}]

Answer 1

Not perfect, but it is a start.

f = X \[Function] {2 Indexed[X, 1] + 
     Sin[Indexed[X, 1] + Indexed[X, 2]], 
    2 Indexed[X, 2] + Cos[Indexed[X, 1] + Indexed[X, 2]]};
Df = X \[Function] Block[{x, y},
    D[f[{x, y}], {{x, y}, 1}] /. {x -> Indexed[X, 1], y -> Indexed[X, 2]}
    ];

g = Y \[Function] Block[{X}, X /. FindRoot[f[X] == Y, {X, 0. Y}]];
Dg = Y \[Function] Block[{x, y, X},
    X = g[Y];
    Inverse[Df[X]]
    ];

StreamDensityPlot[g[{x, y}], {x, -3, 5}, {y, -3, 5}]

Addendum

FindRoot can be very sensitive to the initial guess. By applying f to the points a sufficiently large and fine grid, one can employ Nearest to obtain a ``coarse inverse'' of f that can be refined with FindRoot:

ClearAll[f, Df];
Block[{X,x,y},
  f[X_] = {2 Indexed[X, 1] + Sin[Indexed[X, 1] + Indexed[X, 2]], 
    2 Indexed[X, 2] + Cos[Indexed[X, 1] + Indexed[X, 2]]};
  Df[X_] = Evaluate[
    Block[{x, y}, 
     D[f[{x, y}], {{x, y}, 1}] /. {x -> Indexed[X, 1], y -> Indexed[X, 2]}]
    ];
  ];

m = 100;
n = 100;
xgrid = Subdivide[-1., 3.5, m];
ygrid = Subdivide[-1., 3.5, n];
p = Tuples[{xgrid, ygrid}];
q = Map[f, p];
gcoarse = Nearest[q -> p];
ClearAll[g];
g[Y_?(VectorQ[#, NumericQ] &)] := 
  Block[{X}, X /. FindRoot[f[X] == Y, {X, gcoarse[Y][[1]]}]];

With g[Y_?(VectorQ[#, NumericQ] &)] := , I make also sure that g is evaluated only on numeric data. This should get you rid of most of the error messages.

I also added a few Blocks in order to make things more robust.