Hessian of matrix fractional function I am interested in the matrix fractional function (from $GL_n(\mathbb{R})$ to $\mathbb{R}$) : $S \mapsto x^\top S^{-1} x$, where $x \in \mathbb{R}^n$ is fixed.
By computing $f(S + E) - f(S)$ for a small $E$ I am able to show that its gradient at $S$ is $-S^{-1} x x^\top S^{-1}$, but I am stuck when it comes to computing its Hessian (which is a fourth order tensor). Any hint ?
EDIT: Pushing the Taylor expansion leads me to:
\begin{align}x^\top(S + E)^{-1} x - x^\top S ^{-1}  x &= x^\top(Id + S^{-1}  E)^{-1} 
 S^{-1}  x -x^\top S^{-1}  x 
\\&= x^\top(Id - S^{-1} E + S^{-1}  E S^{-1}E) S^{-1} x -x^\top S^{-1}  x +o(||E||^2)
\\&= -x^\top S^{-1} E S^{-1}  x + x^\top S^{-1}  ES^{-1} E S^{-1} x
\end{align}
and $-x^\top S^{-1} E S^{-1}  x = - tr _,x^\top S^{-1} E S^{-1} x = - tr \, S^{-1} x x^\top S^{-1} E$
which gives me the gradient, and the expression of the second order term, but not the Hessian.
 A: Note that $x\rightarrow 1/x$ is convex only when $x>0$.  Thus, the correct function to study is $f:A\in S^+\rightarrow x^TA^{-1}x$ where $S^+$ is the set of the $n\times n$ symmetric $>0$ matrices. Recall that $S^+$ is a convex open subset of $S$, the set of symmetric matrices..
Proposition. $f$ is convex.
Proof.
 There is a mistake in the calculation of Guy Fsone. Indeed $D^2f_A(H,H)=2x^TA^{-1}HA^{-1}HA^{-1}x$ for every $H\in S$ (recall that $S$ is the tangent space of $S^+$).
Note that $A^{-1}HA^{-1}HA^{-1}=(A^{-1}HA^{-1/2})(A^{-1}HA^{-1/2})^T$ is a symmetric $\geq 0$ matrix. Then, for every $(A,H)$, $D^2f_A(H,H)\geq 0$ and $f$ is convex. $\square$
A: Use the Chain rule $$f(S) = g\circ \mathcal{I}(S)$$
 with $\mathcal{I}(S) =S^{-1}$ and $g(S) =x^\top Sx$  and see that,
$$g: S\mapsto x^\top Sx $$ is linear so $Dg(S)(E) =x^\top Ex$.
We know that $$\mathcal{DI}(S)(V) = -S^{-1}VS^{-1}$$
Hence, for $H $ fixed,
$$\mathcal{D}f(S)(H) =\mathcal{D}g \circ \mathcal{DI}(S)(H) =\color{red}{x^\top S^{-1}HS^{-1} E x :=h\circ \mathcal{I}(S)}$$
Where $$h: A\mapsto -x^\top AHA x$$
Observe that $h$ is quadratic map hence 
$$\mathcal{D}h(A)(K) =  -x^\top KHA x-x^\top AHK x$$
Therefore, 
$$\mathcal{D}^2f(S)(H)(K) =\mathcal{D}h \circ\mathcal{DI}(S)(K) =-x^\top [\mathcal{DI}(S)(K)]HS x-x^\top SH [\mathcal{DI}(S)(K)] x \\
= \mathcal{D}^2f(S)(H)(K) = x^\top SH S^{-1}KS^{-1} x+ x^\top S^{-1}KS^{-1}HS x . $$
i.e 
$$\color{red}{\mathcal{D}^2f(S)(H)(K) = x^\top SH S^{-1}KS^{-1} x+ x^\top S^{-1}KS^{-1}HS x }$$
