Hessian of Log of Matrix-t distribution I am trying to calculate the hessian of the log of the matrix-t distribution. I know that the log of the matrix-t distribution can be written:
$$\log T_{N\times P}(X| \nu, M, \Sigma, \Omega) \propto -\frac{\nu + n+ p-1}{2}\log|I_N + \Sigma^{-1}(X-M)\Omega^{-1}(X-M)^T|$$
where $\nu$ is a scalar, $M$ is an $N\times P$ matrix, $\Sigma$ is a $N\times N$ covariance matrix and $\Omega$ is an $P\times P$ covariance matrix. I also believe that I can calculate the gradient and that it is given by:
$$\frac{d\log T_{N\times P}(X| \nu, M, \Sigma, \Omega)}{dX} \propto -c(\Sigma^{-1}+\Sigma^{-1}(X-M)\Omega^{-1}(X-M)^T\Sigma^{-1})^{-1}(X-M)\Omega^{-1}M^T$$
where I have let $c=\nu + n+ p-1$. 
My issue, I don't know how to massage this into a form to calculate the Hessian. In particular, I am interested in the hessian in the form 
$$\frac{d\log T_{N\times P}(X| \nu, M, \Sigma, \Omega)}{dvec(X)dvec(X)'}$$ although I am flexible. 
 A: For ease of typing, define the variables
$$\eqalign{
S &= \Sigma^{-1},\,\,\,\,W=\Omega^{-1},\,\,\,\,Y=(X-M) \cr
A &=  I+SYWY^T,\,\,\,\,c=\frac{1-\nu-n-p}{2} \cr
\lambda &= c\,\log\det A \cr
}$$
Note that $\,S^T\!=\!S,\,$ $W^T\!=\!W,\,$ $\,\,dY\!=\!dX,\,$ and the function we need to differentiate is $\lambda$.
Find the differential and gradient
$$\eqalign{
d\lambda
 &= cA^{-T}:dA \cr
 &= cA^{-T}:(S\,dY\,WY^T+SYW\,dY^T) \cr
 &= c\big(SA^{-T}YW+A^{-1}SYW\big):dY \cr
 &= c\big(SA^{-T}YW+A^{-1}SYW\big):dX \cr
G=\frac{\partial\lambda}{\partial X}
 &= c\big(SA^{-T}YW+A^{-1}SYW\big) \cr
}$$
Now find the differential of the gradient
$$\eqalign{
dG 
 &= c\big(S\,dA^{-T}YW+SA^{-T}\,dY\,W+dA^{-1}SYW+A^{-1}S\,dY\,W\big) \cr
}$$
Expand the differential terms embedded in the RHS, then apply the vectorization operation, i.e.
$$\eqalign{
{\rm vec}(AXB) &= (B^T\otimes A)\,{\rm vec}(X)\cr
{\rm vec}(X^T) &= K{\rm vec}(X)\cr
}$$ where $K$ is a permutation matrix which is called the Commutation Matrix.
Putting all the pieces together, the Hessian looks like this 
$$\eqalign{
&H = \frac{\partial^2\lambda}{\partial x\partial x^T} \cr
&= c(I\otimes SA^{-T})+c(W\otimes A^{-1}S) \cr
&-c(WY^TA^{-1}YW\otimes SA^{-T})
-c(WY^TSA^{-T}YW\otimes A^{-1}S) \cr
&-c(S\otimes SA^{-T}YW)K
-c(WY^TSA^{-T}\otimes A^{-1}SYW)K
\cr
}$$
A: You ought to be able to numerically evaluate the Hessian using a matrix-level automatic differentiator, such as ADiMat http://www.sc.informatik.tu-darmstadt.de/res/sw/adimat/ and http://adimat.sc.informatik.tu-darmstadt.de/ .
If you only need Hessian-vector products, as opposed to the full Hessian matrix, you can specify an appropriate seed matrix.
A: The key identity you need is
$$
\frac{\partial}{\partial x} A^{-1} = - A^{-1}\frac{\partial A}{\partial x} A^{-1}
$$
which can be obtained by product-ruling $A^{-1}A = 1$.
Now to turn you gradient in a Hessian follow these steps, and don't get discouraged if it gets messy (introducing extra variables should help)


*

*product rule

*apply the "key" identity on $A = (\Sigma^{-1}+\Sigma^{-1}(X-M)\Omega^{-1}(X-M)^T\Sigma^{-1})$, use chain rule

