1
$\begingroup$

$$ \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
$\endgroup$
0
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.