# How to take derivative of multivariate normal density?

Say I have multivariate normal $N(\mu, \Sigma)$ density. I want to get the second (partial) derivative w.r.t. $\mu$. Not sure how to take derivative of a matrix.

Wiki says take the derivative element by element inside the matrix.

I am working with Laplace approximation $$\log{P}_{N}(\theta)=\log {P}_{N}-\frac{1}{2}{(\theta-\hat{\theta})}^{T}{\Sigma}^{-1}(\theta-\hat{\theta}) \>.$$
The mode is $\hat\theta=\mu$.

I was given $${\Sigma}^{-1}=-\frac{{{\partial }^{2}}}{\partial {{\theta }^{2}}}\log p(\hat{\theta }|y),$$ how did this come about?

What I have done:
$$\log P(\theta|y) = -\frac{k}{2} \log 2 \pi - \frac{1}{2} \log \left| \Sigma \right| - \frac{1}{2} {(\theta-\hat \theta)}^{T}{\Sigma}^{-1}(\theta-\hat\theta)$$

So, I take derivative w.r.t to $\theta$, first off, there is a transpose, secondly, it is a matrix. So, I am stuck.

Note: If my professor comes across this, I am referring to the lecture.

• part of your problem may be that your expression for the log-likelihood has an error - you have $|\Sigma|$ where you should have $\log(|\Sigma|)$. Also, by any chance did you mean ${\Sigma}^{-1}=-\frac{{{\partial }^{2}}}{\partial {{\theta }^{2}}}\log p(\theta|y)$? – Macro May 1 '12 at 12:19
• Yes, you are right, sorry. Why is there negative sign in front of the partial derivative? – user1061210 May 1 '12 at 14:38
• I was just clarifying about the negative sign because, the negative second derivative is the observed fisher information, which is usually of interest. Also, by my own calculation, I'm finding that $\frac{{{\partial }^{2}}}{\partial {{\theta }^{2}}}\log p(\theta|y) = -\Sigma^{-1}$ – Macro May 1 '12 at 14:48
• So, what is the general procedure for discrete/continuous function? Take log, write in Taylor expansion form, differentiate twice w.r.t. $\theta$. Fisher info is not generally true most other densities, right? – user1061210 May 1 '12 at 15:19
• @user As I pointed out, the second derivative of the logarithm must have non-positive eigenvalues. Yes, there are links between variances and negative second partial derivatives, as the theory of maximum likelihood estimation, Fisher information, etc., reveals--Macro has referred to that earlier in these comments. – whuber May 1 '12 at 19:19

In chapter 2 of the Matrix Cookbook there is a nice review of matrix calculus stuff that gives a lot of useful identities that help with problems one would encounter doing probability and statistics, including rules to help differentiate the multivariate Gaussian likelihood.

If you have a random vector ${\boldsymbol y}$ that is multivariate normal with mean vector ${\boldsymbol \mu}$ and covariance matrix ${\boldsymbol \Sigma}$, then use equation (86) in the matrix cookbook to find that the gradient of the log likelihood ${\bf L}$ with respect to ${\boldsymbol \mu}$ is

\begin{align} \frac{ \partial {\bf L} }{ \partial {\boldsymbol \mu}} &= -\frac{1}{2} \left( \frac{\partial \left( {\boldsymbol y} - {\boldsymbol \mu} \right)' {\boldsymbol \Sigma}^{-1} \left( {\boldsymbol y} - {\boldsymbol \mu}\right) }{\partial {\boldsymbol \mu}} \right) \nonumber \\ &= -\frac{1}{2} \left( -2 {\boldsymbol \Sigma}^{-1} \left( {\boldsymbol y} - {\boldsymbol \mu}\right) \right) \nonumber \\ &= {\boldsymbol \Sigma}^{-1} \left( {\boldsymbol y} - {\boldsymbol \mu} \right) \end{align}

I'll leave it to you to differentiate this again and find the answer to be $-{\boldsymbol \Sigma}^{-1}$.

As "extra credit", use equations (57) and (61) to find that the gradient with respect to ${\boldsymbol \Sigma}$ is

\begin{align} \frac{ \partial {\bf L} }{ \partial {\boldsymbol \Sigma}} &= -\frac{1}{2} \left( \frac{ \partial \log(|{\boldsymbol \Sigma}|)}{\partial{\boldsymbol \Sigma}} + \frac{\partial \left( {\boldsymbol y} - {\boldsymbol \mu}\right)' {\boldsymbol \Sigma}^{-1} \left( {\boldsymbol y}- {\boldsymbol \mu}\right) }{\partial {\boldsymbol \Sigma}} \right)\\ &= -\frac{1}{2} \left( {\boldsymbol \Sigma}^{-1} - {\boldsymbol \Sigma}^{-1} \left( {\boldsymbol y} - {\boldsymbol \mu} \right) \left( {\boldsymbol y} - {\boldsymbol \mu} \right)' {\boldsymbol \Sigma}^{-1} \right) \end{align}

I've left out a lot of the steps, but I made this derivation using only the identities found in the matrix cookbook, so I'll leave it to you to fill in the gaps.

I've used these score equations for maximum likelihood estimation, so I know they are correct :)

• Great reference - was going to recommend it myself. Not a good pedagogical reference for someone who doesn't know matrix algebra though. The real challenge comes from actually working out $\Sigma$. A real pain. – probabilityislogic May 1 '12 at 11:11
• Another good source on matrix calculus is Magnus & Neudecker, amazon.com/… – StasK May 1 '12 at 15:32
• The equation's reference number has been changed (maybe due to a new edition). The new reference equation is 86. – goelakash May 17 '16 at 9:22
• I could be off-base here but I don't think this formula is correct. I've been using this with real examples and looking at their finite differences. It seems that the formula for $\frac{ \partial {\bf L} }{ \partial {\boldsymbol \Sigma}}$ gives the correct values for the diagonal entries. However, the off-diagonal entries are half of what they should be. – jjet Apr 9 '18 at 19:27

You need to make sure you properly take care of the repeated elements in $\mathbf{\Sigma}$, otherwise you're derivatives will be incorrect. For example, (141) the Matrix Cookbook gives for a symmetric $\mathbf{\Sigma}$ the following derivatives

\begin{align} \frac{\partial \log|\mathbf{\Sigma}|}{\partial \mathbf{\Sigma}}&=2\mathbf{\Sigma}^{-1}-(\mathbf{\Sigma}^{-1}\circ I) \end{align}

And (14) of Differentiation of functions of covariance matrices gives \begin{align} \frac{\partial \textrm{trace}(\mathbf{\Sigma}^{-1}\mathbf{x}\mathbf{x}^\top)}{\partial \mathbf{\Sigma}}&=-2\mathbf{\Sigma}^{-1}\mathbf{x}\mathbf{x}^\top\mathbf{\Sigma}^{-1}+(\mathbf{\Sigma}^{-1}\mathbf{x}\mathbf{x}^\top\mathbf{\Sigma}^{-1}\circ I) \end{align}

where $\circ$ denotes the Hadmard product and for convenience we have defined $\mathbf{x}:=\mathbf{y}-\mathbf{\mu}$.

Note in particular this is not the same as when symmetricity of $\mathbf{\Sigma}$ is not imposed. As a result we have that

\begin{align} \frac{\partial \mathbf{L}}{\partial \mathbf{\Sigma}}&=-\frac{\partial }{\partial \mathbf{\Sigma}}\frac{1}{2}\left(D\log|2\pi|+ \log|\mathbf{\Sigma}| + \mathbf{x}^{\top}\mathbf{\Sigma}^{-1}\mathbf{x})\right)\\ &=-\frac{\partial }{\partial \mathbf{\Sigma}}\frac{1}{2}\left( \log|\mathbf{\Sigma}| + \textrm{trace}(\mathbf{\Sigma}^{-1}\mathbf{x}\mathbf{x}^\top)\right)\\ &=-\frac{1}{2}\left( 2\mathbf{\Sigma}^{-1}-(\mathbf{\Sigma}^{-1}\circ I) -2\mathbf{\Sigma}^{-1}\mathbf{x}\mathbf{x}^\top\mathbf{\Sigma}^{-1}+(\mathbf{\Sigma}^{-1}\mathbf{x}\mathbf{x}^\top\mathbf{\Sigma}^{-1}\circ I)\right) \end{align}

where $D$ denotes the dimension of $\mathbf{x}$, $\mathbf{y}$ and $\mathbf{\mu}$ and the derivative of $D\log|2\pi|$ is 0

This ensures the $i,j^{th}$ element of $\frac{\partial \mathbf{L}}{\partial \mathbf{\Sigma}}$ corresponds to $\frac{\partial \mathbf{L}}{\partial \mathbf{\Sigma}_{ij}}$.

I tried to computationally verify @Macro's answer but found what appears to be a minor error in the covariance solution. He obtained \begin{align} \frac{ \partial {\bf L} }{ \partial {\boldsymbol \Sigma}} &= -\frac{1}{2} \left( {\boldsymbol \Sigma}^{-1} - {\boldsymbol \Sigma}^{-1} \left( {\boldsymbol y} - {\boldsymbol \mu} \right) \left( {\boldsymbol y} - {\boldsymbol \mu} \right)' {\boldsymbol \Sigma}^{-1} \right) ={\bf A} \end{align} However, it appears that the correct solution is actually $${\bf B}=2{\bf A} - \text{diag}({\bf A})$$ The following R script provides a simple example in which the finite difference is calculated for each element of ${\boldsymbol \Sigma}$. It demonstrates that ${\bf A}$ provides the correct answer only for diagonal elements while ${\bf B}$ is correct for every entry.

library(mvtnorm)

set.seed(1)

# Generate some parameters
p <- 4
mu <- rnorm(p)
Sigma <- rWishart(1, p, diag(p))[, , 1]

# Generate an observation from the distribution as a reference point
x <- rmvnorm(1, mu, Sigma)[1, ]

# Calculate the density at x
f <- dmvnorm(x, mu, Sigma)

# Choose a sufficiently small step-size
h <- .00001

# Calculate the density at x at each shifted Sigma_ij
f.shift <- matrix(NA, p, p)
for(i in 1:p) {
for(j in 1:p) {
zero.one.mat <- matrix(0, p, p)
zero.one.mat[i, j] <- 1
zero.one.mat[j, i] <- 1

Sigma.shift <- Sigma + h * zero.one.mat
f.shift[i, j] <- dmvnorm(x, mu, Sigma.shift)
}
}

# Caluclate the finite difference at each shifted Sigma_ij
fin.diff <- (f.shift - f) / h

# Calculate the solution proposed by @Macro and the true solution
A <- -1/2 * (solve(Sigma) - solve(Sigma) %*% (x - mu) %*% t(x - mu) %*% solve(Sigma))
B <- 2 * A - diag(diag(A))

# Verify that the true solution is approximately equal to the finite difference
fin.diff
A * f
B * f

• Thank you for your comment. I believe you interpret the notation differently than everyone else has, because you simultaneously change pairs of matching off-diagonal elements of $\Sigma$, thereby doubling the effect of the change. In effect you are computing a multiple of a directional derivative. There does appear to be a small problem with Macro's solution insofar as a transpose ought to be taken--but that would change nothing in the application to symmetric matrices. – whuber Apr 10 '18 at 21:08