Observed information matrix with multivariate normal distribution

Question

$$ \DeclareMathOperator\tr{tr} \DeclareMathOperator\vecOP{vec} \newcommand\di{\mathrm{d}} \newcommand\D{\mathrm{D}} \newcommand\Hess{\mathrm{H}} $$

I do not have much experience with matrix differentials and thus my question is whether my derivation of the observed information matrix is correct.

We consider a single log-likelihood term for a pair $(\vec y_i, \vec x_i)$ in a multivariate normal distribution

\begin{align*} \mathcal L &= - \frac 12 \log \lvert Q \rvert - \frac 12 \tr Q^{-1} (\vec y_i - F^\top \vec x_i)(\vec y_i - F^\top\vec x_i)^\top + \dots \\ &= - \frac 12 \log \lvert Q \rvert - \frac 12 \tr Q^{-1} Z + \dots \\ Z &= (\vec y_i - F^\top \vec x_i)(\vec y_i - F^\top\vec x_i)^\top \end{align*}

The first differential w.r.t. $F$ and $Q$ are

$$\di\mathcal L = \frac 12\tr (\di Q) Q^{-1} (Z - Q)Q^{-1} + \tr Q^{-1}(\vec y_i - F^\top \vec x_i)\vec x_i^\top (\di F) $$

Thus, the derivative is

$$ \D\mathcal L = \left( \vecOP\left(\vec x_i(\vec y_i - F^\top \vec x_i)^\top Q^{-1}\right)^\top, \frac 12\vecOP\left(Q^{-1} (Z - Q)Q^{-1}\right)^\top \right) $$

as in this answer. The second differential is

\begin{align*} \di^2\mathcal L &= \frac 12 \tr (\di Q)(\di Q^{-1})(Z - Q)Q^{-1} + \frac 12 \tr (\di Q)Q^{-1}(\di Z - \di Q)Q^{-1} \\ &\hspace{20pt} + \frac 12 \tr (\di Q)Q^{-1}(Z - Q)(\di Q^{-1}) + \tr (\di Q^{-1})(\vec y_i - F^\top \vec x_i)\vec x_i^\top (\di F) \\ &\hspace{20pt} - \tr Q^{-1}(\di F)^\top \vec x_i \vec x_i^\top (\di F) \\ % % % &= -\frac 12\tr (\di Q)Q^{-1}(\di Q)Q^{-1}(Z - Q)Q^{-1} - \tr (\di Q)Q^{-1}(\vec y_i - F^\top \vec x_i)\vec x_i^\top (\di F)Q^{-1} \\ &\hspace{20pt} - \frac 12 \tr (\di Q)Q^{-1}(\di Q)Q^{-1} - \frac 12 \tr (\di Q)Q^{-1}(Z - Q)Q^{-1}(\di Q)Q^{-1} \\ &\hspace{20pt} - \tr Q^{-1}(\di Q)Q^{-1}(\vec y_i - F^\top \vec x_i)\vec x_i^\top (\di F) - \tr Q^{-1}(\di F)^\top \vec x_i \vec x_i^\top (\di F) \\ % % % &= -\tr (\di Q)Q^{-1}(\di Q)Q^{-1}(Z - \frac 12 Q)Q^{-1} - 2\tr (\di Q)Q^{-1}(\vec y_i - F^\top \vec x_i)\vec x_i^\top (\di F)Q^{-1} \\ &\hspace{20pt} - \tr Q^{-1}(\di F)^\top \vec x_i \vec x_i^\top (\di F) \end{align*}

Thus, the Hessian is

$$ \Hess\mathcal L = \begin{pmatrix} -Q^{-1}\otimes \vec x_i \vec x_i^\top & -Q^{-1} \otimes \left(\vec x_i (\vec y_i - F^\top \vec x_i)^\top Q^{-1}\right) \\ -Q^{-1} \otimes \left(Q^{-1}(\vec y_i - F^\top \vec x_i)\vec x_i^\top \right) & - Q^{-1} \otimes \left(Q^{-1}\left(Z - \frac 12 Q\right)Q^{-1} \right) \end{pmatrix} $$

This code seems to confirm the above

# define parameters
options(digits = 3)
y <- c(-.1, .2)
x <- c(-.5, .25)
Qmat <- matrix(c(1 , .4 ,  .4, .5), 2)
Fmat <- matrix(c(.8, .2, .1, .5), 2)
library(mvtnorm)

func <- function(vals){
  Fmat[] <- vals[1:4]
  Qmat[upper.tri(Qmat, diag = TRUE)] <- vals[5:7]
  Qmat[lower.tri(Qmat)] <- t(Qmat)[lower.tri(Qmat)]

  if(any(abs(eigen(Fmat)$values) >= 1) || any(eigen(Qmat)$values <= 0))
    return(NA_real_)
  dmvnorm(y, crossprod(Fmat, x), Qmat, log = TRUE)
}

# compute Hessian with numerical differentiation
library(numDeriv)
o <- numDeriv::hessian(func, c(Fmat, Qmat[lower.tri(Qmat, diag = TRUE)]))

# compare block diagonal elements
library( matrixcalc)
D <- duplication.matrix(2L)
Z <- tcrossprod(y - crossprod(Fmat, x))
K <- solve(Qmat)
o[1:4, 1:4] - (- kronecker(K, tcrossprod(x)))
#R           [,1]      [,2]      [,3]      [,4]
#R [1,]  3.53e-12 -1.67e-13 -3.71e-12 -5.01e-12
#R [2,] -1.67e-13  8.34e-14  4.46e-11 -6.68e-14
#R [3,] -3.71e-12  4.46e-11 -1.62e-11 -3.34e-13
#R [4,] -5.01e-12 -6.69e-14 -3.34e-13  4.27e-12

lower_diag_block <- -kronecker(K, solve(Qmat, Z %*% K - diag(1/2, 2)))
o[5:7, 5:7] - crossprod(D, lower_diag_block %*% D)
#R           [,1]      [,2]      [,3]
#R [1,] -7.13e-12  8.80e-12 -9.98e-10
#R [2,]  8.80e-12 -1.99e-11  3.09e-11
#R [3,] -9.98e-10  3.09e-11 -4.40e-11

# then the off diagonal blocks
off_block <-
  - kronecker(K, tcrossprod(x, y - crossprod(Fmat, x)) %*% K)
o[1:4, 5:7] - off_block %*% D
#R           [,1]      [,2]      [,3]
#R [1,] -2.77e-12  3.81e-12  3.12e-12
#R [2,]  5.49e-12 -9.39e-12 -2.55e-12
#R [3,]  1.12e-11  1.36e-11  4.85e-12
#R [4,] -4.04e-12  3.24e-12  2.33e-12
off_block <-
  - kronecker(K, solve(Qmat, tcrossprod(y - crossprod(Fmat, x), x)))
o[5:7, 1:4] - crossprod(D, off_block)
#R           [,1]      [,2]     [,3]      [,4]
#R [1,] -2.77e-12  5.49e-12 1.12e-11 -4.04e-12
#R [2,]  3.81e-12 -9.39e-12 1.36e-11  3.24e-12
#R [3,]  3.12e-12 -2.55e-12 4.85e-12  2.33e-12

Answer 1

The original solution I gave is wrong. To illustrate this, we define the log-likelihood as

$$\newcommand{\fhalf}{\frac 12}\DeclareMathOperator{\tr}{tr} \begin{align*} L&\propto -\fhalf \log \lvert Q \rvert - \fhalf (y - Fx)^\top Q^{-1}(y - Fx) \\ &= -\fhalf \log \lvert Q \rvert - \fhalf \tr Q^{-1}\underbrace{(y - Fx)(y - Fx)^\top}_Z \\ \end{align*}$$

so we now use $F$ instead of $F^\top$. The differential is

\begin{align*} dL &= - \fhalf \tr Q^{-1} dQ + \fhalf \tr Q^{-1} (dQ) Q^{-1}Z - \fhalf \tr Q^{-1}dZ \\ &= -\fhalf \tr Q^{-1} dQ + \fhalf \tr Q^{-1} (dQ) Q^{-1}Z + \tr x (y - Fx)^\top Q^{-1}dF \\ &= \fhalf \tr Q^{-1}(Z - Q)Q^{-1}dQ + \tr x (y - Fx)^\top Q^{-1}dF \end{align*}

Next, let

\begin{align*} dL &= \underbrace{\fhalf \tr Q^{-1}ZQ^{-1}dQ}_{Z_1} - \underbrace{\fhalf \tr Q^{-1}dQ}_{z_2} + \underbrace{\tr x (y - Fx)^\top Q^{-1}dF}_{z_3} \end{align*}

then we derive the following to find the second differential

\begin{align*} dz_1 &= \fhalf \tr (dQ^{-1}) ZQ^{-1}dQ + \fhalf \tr Q^{-1}Z(dQ^{-1}) dQ + \fhalf \tr Q^{-1}(dZ) Q^{-1}dQ \\ % &= -\fhalf \tr Q^{-1}(dQ) Q^{-1}ZQ^{-1}dQ - \fhalf \tr Q^{-1}ZQ^{-1}(dQ) Q^{-1}dQ \\ &\hspace{20pt} - \tr Q^{-1}(dF) x(y - Fx)^\top Q^{-1}dQ \\ % dz_2 &= \fhalf \tr (dQ^{-1}) dQ = -\fhalf \tr Q^{-1}(dQ) Q^{-1}dQ \\ % dz_3 &= -\tr xx^\top(dF)^\top Q^{-1}dF + \tr x (y - Fx)^\top (dQ^{-1})dF \\ % &= -\tr xx^\top(dF)^\top Q^{-1}dF - \tr x (y - Fx)^\top Q^{-1}(dQ)Q^{-1}dF \end{align*}

Thus, the second differential is

\begin{align*} d^2L &= -\fhalf \tr Q^{-1}(dQ) Q^{-1}ZQ^{-1}dQ - \fhalf \tr Q^{-1}ZQ^{-1}(dQ) Q^{-1}dQ \\ &\hspace{20pt} - \tr Q^{-1}(dF) x(y - Fx)^\top Q^{-1}dQ \\ &\hspace{20pt} + \fhalf \tr Q^{-1}(dQ) Q^{-1}dQ \\ &\hspace{20pt} - \tr xx^\top(dF)^\top Q^{-1}dF - \tr x (y - Fx)^\top Q^{-1}(dQ)Q^{-1}dF \\ % &= -\fhalf \tr Q^{-1}(dQ) Q^{-1}(Z - Q)Q^{-1}dQ - \fhalf \tr Q^{-1}ZQ^{-1}(dQ) Q^{-1}dQ \\ &\hspace{20pt} - \tr Q^{-1}(dF) x(y - Fx)^\top Q^{-1}dQ - \tr x (y - Fx)^\top Q^{-1}(dQ)Q^{-1}dF \\ &\hspace{20pt} - \tr xx^\top(dF)^\top Q^{-1}dF \end{align*}

Hence,

\begin{align*} \frac{\partial^2 L}{\partial F\partial F^\top} &= -xx^\top \otimes Q^{-1} \\ \frac{\partial^2 L}{\partial Q\partial F^\top} &= - \fhalf \left((Q^{-1}(y-Fx)x^\top \otimes Q^{-1}) + (Q^{-1}\otimes Q^{-1}(y-Fx)x^\top)K_{pp}\right) \\ \frac{\partial^2 L}{\partial Q\partial Q^\top} &= -\fhalf K_{pp} \left(Q^{-1} \otimes (Q^{-1}(Z - Q)Q^{-1}) - (Q^{-1}ZQ^{-1}) \otimes Q^{-1}\right) \end{align*}

where $K_{pp}$ is the $p$ by $p$ commutation matrix and we assume that both $x$ and $y$ are $p$-dimensional. I am though slightly confused about $\frac{\partial^2 L}{\partial Q\partial Q^\top}$. However, the code below shows that the above is correct and the solution in my question is wrong

# 3D example
x <- c(-3, 2, 1)
y <- c(-1, 0, 2)
F. <- matrix(c(.8, .2, .3, .1, .6, .2, .1, .1, .5), 3L)
Q <- matrix(c(3, 1, 1, 1, 2, 3, 1, 3, 7), 3L)
library(mvtnorm)
library(numDeriv)

# log-likelihood (which will not fail even if Q is not symmetric -- slightly
# odd)
func <- function(fq){
  F.[] <- fq[1:9]
  Q[] <- fq[-(1:9)]
  rss <- y - F. %*% x
  rss <- crossprod(rss, solve(Q, rss))
  c(-log(2 * pi) * 3 / 2 - determinant(Q)$modulus / 2 - drop(rss) / 2)
}

# sanity check...
func(c(F., Q))
#R [1] -6.866
dmvnorm(y, F. %*% x, Q, log = TRUE)
#R [1] -6.866

# compute Hessian
hm <- matrix(hessian(func, c(F., Q)), 18L, 18L)

# check that upper left block is correct
isTRUE(all.equal(hm[1:9, 1:9], -kronecker(tcrossprod(x), solve(Q))))
#R [1] TRUE

# check that lower left block is correct
library(matrixcalc)
cum <- commutation.matrix(3L, 3L)
tmp <- solve(Q, y - F. %*% x)

isTRUE(all.equal(
  hm[-(1:9), 1:9],
  -.5 * (
    kronecker(tcrossprod(tmp, x), solve(Q)) +
      kronecker(solve(Q), tcrossprod(tmp, x)) %*% cum), tolerance = 1e-4))
#R [1] TRUE

# check that lower right block is correct
tmp <- tcrossprod(solve(Q, y - F. %*% x))
res <- cum %*% (
  kronecker(solve(Q), solve(Q) - tmp) / 2 - kronecker(tmp, solve(Q)) / 2)
isTRUE(all.equal(hm[10:18, 10:18], res, tolerance = 1e-4))
#R [1] TRUE

# the first answer is wrong. To show this, we redfine the log-likelihood 
# function to take F^\top as input
func <- function(fq){
  F.[] <- fq[1:9]
  Q[] <- fq[-(1:9)]
  rss <- y - crossprod(F., x)
  rss <- crossprod(rss, solve(Q, rss))
  -log(2 * pi) * 3 / 2 - determinant(Q)$modulus / 2 - drop(rss) / 2
}

# compute Hessian
hm <- matrix(hessian(func, c(t(F.), Q)), 18L, 18L)

# upper left block is correct
all.equal(hm[1:9, 1:9], -kronecker(solve(Q), tcrossprod(x)))
#R [1] TRUE

# lower left block is wrong
all.equal(
  hm[-(1:9), 1:9],
  R1 <- -kronecker(solve(Q), tcrossprod(solve(Q, y - F. %*% x), x)))
#R [1] "Mean relative difference: 0.1455"

# but correct if you only consider the lower diagonal part
D <- duplication.matrix(3L)
all.equal(
  crossprod(D, hm[-(1:9), 1:9]), crossprod(D, R1), tolerance = 1e-4)
#R [1] TRUE

# ditto w/ lower right block
all.equal(
  hm[-(1:9), -(1:9)],
  R2 <- -kronecker(solve(Q), solve(Q, tcrossprod(y - F. %*% x) %*% solve(Q) -
                                     diag(.5, 3))))
#R [1] "Mean relative difference: 0.2664"

all.equal(
  crossprod(D, hm[-(1:9), -(1:9)]) %*% D,
  crossprod(D, R2) %*% D, tolerance = 1e-4)
#R [1] TRUE

Sorry for the long answer.