I am interested in the residual covariance of a multivariate regression model. The regression is
$$ Y_t = X_t \beta + \varepsilon_t $$
and I have a log likelihood as follows $$ \mathcal{L}(\Sigma) = \frac{T}{2} \log | \Sigma^{-1} | - \frac{1}{2} \sum_{t=1}^T \varepsilon_t^{\top} \Sigma^{-1} \varepsilon_t $$ omitting constants and assuming $\Sigma$ is the residual covariance matrix. Additionally, $Y_t \in \mathbb{R}^{N}, X_t \in \mathbb{R}^{N}, \varepsilon_t \in \mathbb{R}^{N}, \beta \in \mathbb{R}^{N \times N}$ are the dimensions of the variables. I can obtain the gradient with respect to $\Sigma$ as follows $$ \frac{\partial}{\partial \Sigma} \mathcal{L}(\Sigma) = -\frac{T}{2} \Sigma^{-1} + \frac{1}{2} \Sigma^{-1} \varepsilon_t^\top \varepsilon_t \Sigma^{-1} $$ My question is, is there a closed form for the hessian of the log likelihood? More specifically, I want to obtain the maximum eigenvalue $\lambda_1$ of the hessian to obtain the Lipschitz constant for a gradient descent algorithm. This is now more difficult because we have to take derivatives with respect to inverses of matrices instead of scalar.
I believe the result is a four dimensional tensor, but this can be flatted to obtain a Kronecker product of dimension $N^2 \times N^2$. Then I want to obtain the maximum eigenvalue of this form.
