An expression for computing second order partial derivatives of an implicitely defined function Let $\Phi(x,y)=0$ be an implicit function s.t. $\Phi:\mathbb{R}^n\times \mathbb{R}^k\rightarrow \mathbb{R}^n$ and $\det\left(\frac{\partial \Phi}{\partial 
 x}(x_0,y_0)\right)\neq 0$. This means that locally at $(x_0,y_0)$ we can express $x_i$ as functions of $y$.
Next, we can compute partial derivatives of $x$ as 
\begin{equation}\tag{*}\frac{\partial x_i}{\partial y_j}=-\frac{\det\left(\left[\frac{\partial \Phi}{\partial 
 x_1},\dots,\frac{\partial \Phi}{\partial 
 x_{i-1}}, \frac{\partial \Phi}{\partial 
 y_j}, \frac{\partial \Phi}{\partial 
 x_{i+1}},\dots, \frac{\partial \Phi}{\partial 
 x_n}\right]\right)}{\det\left(\frac{\partial \Phi}{\partial 
 x}\right)}.\end{equation}
This is known. What I wonder is: 
Q: is it possible to compute second order partial derivatives in a systematic way? 
I tried to differentiate determinants using the Jacobi formula, but this leads to very complicated expressions that I cannot handle. I also expanded the determinants in ($*$) along the $i$ column (by which the respective matrices differ) and tried some other approaches, but they do not seem to bring me any further. 
On the other hand, if a go a straightforward way and differentiate $\Phi(x,y)$ twice, I get expressions involving tensors or rather multiindex notations, because neither second order partial derivatives, nor the derivatives of type $\frac{\partial x^i}{\partial y^j}$ are actually tensors. 
My hope is that maybe it is still possible to extract some nice tractable expression similar to how we got ($*$) from $\frac{\partial x}{\partial y}=-\left[\frac{\partial \Phi}{\partial x}\right]^{-1}\frac{\partial \Phi}{\partial y}$?
Here is a related question.
UPDATE: It seems that the problem turned out to be more difficult than I expected (although many people told me that it must have been solved by somebody). Since the hope for getting a resolutive answer fades and the bounty will expire in a couple of days, I'd gladly grant it to anybody who could point out a way to approach (if not solve) this problem. 
UPDATE 2: Let me expand a bit on the above. To illustrate my problem let's differentiate $\left[\frac{\partial \Phi}{\partial x}\right]^{-1}$ w.r.t. $y_i$:
\begin{multline*}\frac{\partial}{\partial y_i}\left[\frac{\partial \Phi}{\partial x}\right]^{-1}=-\left[\frac{\partial \Phi}{\partial x}\right]^{-1}\frac{\partial}{\partial y_i}\left[\frac{\partial \Phi}{\partial x}\right]\left[\frac{\partial \Phi}{\partial x}\right]^{-1}\\
=-\left[\frac{\partial \Phi}{\partial x}\right]^{-1}\left[\frac{\partial^2 \Phi}{\partial x\partial x}\right]\frac{\partial x}{\partial y_i}\left[\frac{\partial \Phi}{\partial x}\right]^{-1}-\left[\frac{\partial \Phi}{\partial x}\right]^{-1}\left[\frac{\partial^2 \Phi}{\partial y_i\partial x}\right]\left[\frac{\partial \Phi}{\partial x}\right]^{-1}.\end{multline*}
So, what is $\left[\frac{\partial^2 \Phi}{\partial x\partial x}\right]\frac{\partial x}{\partial y_i}$? A 3D matrix multiplied with a vector? How to treat these expressions? To make the things even more complicated we should now substitute $\frac{\partial x}{\partial y_i}$ with the respective expression for the first order partial derivatives. It becomes completely obscure and I cannot recognize any structure in it.
 A: I'll post a partial answer. 
Pretend your functions are given by Taylor series to the needed (second) order. 
So, we write 
$$x_p=h_p(y)=\sum \frac{\partial h_p}{\partial y_k} y_k + \sum_{i,j} \frac{1}{2}\frac{\partial^2 h_p}{\partial y_i \partial y_j} y_i y_j$$
$$\Phi_l(x,y)=  \sum_k \frac{\partial \Phi_l}{\partial x_k} x_k + \sum_p \frac{\partial \Phi_l}{\partial y_p} y_p +  \sum_{i,j} \frac{1}{2}\frac{\partial^2 \Phi_l}{\partial x_i \partial x_j} x_i x_j+ \sum_{q,r} \frac{1}{2}\frac{\partial^2 \Phi_l}{\partial y_q  \partial y_r} y_q y_r+\sum_{p,k} \frac{1}{2} \frac{\partial^2 \Phi_l}{\partial y_p  \partial x_k} y_p x_k+
\sum_{p,k} \frac{1}{2} \frac{\partial^2 \Phi_l}{\partial x_k  \partial y_p} x_k y_p $$
Now plug in and keep equate the coefficients of $y_i y_j$ to zero. There are 5 terms in $\Phi_l$. They contribute (in the case of $i\neq j$, so summing the "$y_iy_j$" and the "$y_jy_i$" contributions; there are 1/2 factors throughout if $i=j$):
1) Nothing.
2) $\sum_p \frac{\partial \Phi_l}{\partial y_p} \frac{\partial^2 h_p}{\partial y_i \partial y_j} $ 
3) $\frac{\partial^2 \Phi_l}{\partial y_i \partial y_j}$
4)$ \sum_{q,r} \frac{\partial^2 \Phi_l}{\partial x_q  \partial x_r} \frac{\partial h_q}{\partial y_i} \frac{\partial h_r}{\partial y_j}$
5)$\sum_{p}  \frac{\partial^2 \Phi_l}{\partial x_p  \partial y_j} \frac{\partial h_p}{\partial y_i} $ 
6)$\sum_{p}  \frac{\partial^2 \Phi_l}{  \partial y_i \partial y_p} \frac{\partial h_p}{\partial y_j} $
Varying $l$, one gets $n$  linear equations (labeled by $l$) in $n$ unknowns $\frac{\partial^2 h_p}{\partial y_i \partial y_j} $ (labeled by $p$), which can therefore be written as $A z=b$, with the matrix $A$ of the linear system given by $\Phi_{x}$. Hence these equations can be solved. The only trouble is in writing the $b$ vector, which is the sum of terms 3-6, in a "vector" format.
Maybe a better way to do bookkeeping for this is to use tree-speak like here... 
