Let have a function $f:R^{n}\longrightarrow R^{n}$ given by $x_{i}\stackrel{f}{\longrightarrow}y_{i}=y_{i}(x_{1},\ldots,x_{n})$, where $f$ is smooth and non-singular in the domain of interest. The components of the Jacobian of $f$ is given by $$\left[\mathrm{J}_{f}\right]_{ij}=\frac{\partial y_{i}}{\partial x_{j}},$$and the components Jacobian of the Inverse transformation are
$$\left[\mathrm{J}_{f^{-1}}\right]_{ij}=\frac{\partial x_{i}}{\partial y_{j}}.$$
According to the inverse function theorem, we have the nice property
$$\mathrm{J}_{f^{-1}}=\mathrm{J}_{f}^{-1}\leftrightarrow\Sigma\frac{\partial y_{i}}{\partial x_{k}}\frac{\partial x_{k}}{\partial y_{j}}=\delta_{ij}$$.
Now, the 3rd rank tensor formed with the Hessian of each entry of $f$ is
$$\mathrm{\left[H_{\mathit{f}}\right]}_{jk}^{i}=\frac{\partial^{2}y_{i}}{\partial x_{j}\partial x_{k}}.$$
The question is: what can be said about $\mathrm{\left[H_{\mathit{f^{-1}}}\right]}_{jk}^{i}=\frac{\partial^{2}x_{i}}{\partial y_{j}\partial y_{k}}$? does it have nice identities when the indices are contracted with $\mathrm{H}_{\mathit{f}}$? I've been unable to find any reference about it.
Suggestion. Try expanding each y coordinate as a Taylor expansion of order two in the variables x:
$ y^i = a^i_j x^j + q^i_{m,n} x^m x^n$ using the usual summation convention for repeated subscripts. (The quadratic form $q$ is the Hessian and the matrix $a$ is the Jacobian).
Then do the same after reversing the role of the variables. (Here capital letters indicate that we are discussing the inverse relations ):
$x^I = A^I_J Y^J + Q^I_{M,N} y^My^N$.
Then do some algebra: compose one set of equations with the other set, to represent the identity map as an algebraic combination of the two sets of Taylor coefficients (of course retaining only terms up to order two). You will get something like this multivariable equation:
$x= x$ + sum of two quadratics in $x$.
Explicitly, you get
$x^I=A^I_J ( a^J_j x^j + q^J_{m,n} x^mx^n) + Q^I_{M,N} ( a^M_n x^m +\ldots) (a^N_n x^n+\ldots)$
which after using the fact that the Jacobian matrices $A^I_J$ and $a^i_j$ are inverses of one another, simplifies to
$x^I= x^I + A^I_J q^J_{m,n} + Q^I_{M,N} a^M_n a^N_n x^m x^n +\ldots$.
Thus the algebraic compatibility condition you seek is that these two quadratics are self-cancelling:
$(A^I_J q^J_{m,n} + Q^I_{M,N} a^M_m a^N_n) x^m x^n =0$
Finally you can apply the inverse of the matrix $A$ to solve for $q$ in terms of $Q$.