Hessian of Log of Matrix-t distribution

Question

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.

Answer 1

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 }$$

Answer 2

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.

Answer 3

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)

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