A Fast and Numerically Stable Implementation of the Multivariate Normal PDF

Consider the multivariate normal probability density function (PDF) for $\mathbf{x} \in \mathbb{R}^D$ with parameters $\boldsymbol{\mu} \in \mathbb{R}^D$ (mean) and $\boldsymbol{\Sigma} \in \mathbb{R}^{D \times D}$ (covariance):

$f(\mathbf{x}; \boldsymbol{\mu}, \boldsymbol{\Sigma}) = \frac{1}{\sqrt{(2 \pi)^D \det(\boldsymbol{\Sigma})}} \exp\Big\{ -\frac{1}{2} (\mathbf{x} - \boldsymbol{\mu})^{\top} \boldsymbol{\Sigma}^{-1} (\mathbf{x} - \boldsymbol{\mu}) \Big\}.$

Computing this is tricky because we need to compute both the inverse and determinant of the $D \times D$ covariance matrix $\boldsymbol{\Sigma}$. Each of these operations has runtime complexity $\mathcal{O}(D^3)$. Furthermore, computing matrix inverses and determinants can be numerically unstable if done naively. SciPy has a fast and numerically stable implementation that is worth understanding. The big idea is to do one intensive operation, eigenvalue decomposition, and then use that decomposition to compute the matrix inverse and determinant cheaply.

Matrix inverse

Since $\boldsymbol{\Sigma}$ is Hermitian, it has an eigendecomposition

$\boldsymbol{\Sigma} = \mathbf{Q} \boldsymbol{\Lambda} \mathbf{Q}^{\top}$

where $\mathbf{Q}$ is an orthogonal matrix whose columns are the eigenvectors of $\boldsymbol{\Sigma}$ and where $\boldsymbol{\Lambda}$ is a diagonal matrix of the associated eigenvalues. Since $(\mathbf{A}\mathbf{B})^{-1} = \mathbf{B}^{-1} \mathbf{A}^{-1}$ in general and since $\mathbf{Q}^{\top} = \mathbf{Q}^{-1}$ due to orthogonality, we can easily compute the inverse of $\boldsymbol{\Sigma}$,

$\begin{aligned} \boldsymbol{\Sigma}^{-1} &= (\mathbf{Q} \boldsymbol{\Lambda} \mathbf{Q}^{\top})^{-1} \\ &= (\mathbf{Q}^{\top})^{-1} \boldsymbol{\Lambda}^{-1} \mathbf{Q}^{-1} \\ &= \mathbf{Q} \boldsymbol{\Lambda}^{-1} \mathbf{Q}^{\top}. \end{aligned}$

Since $\boldsymbol{\Lambda}$ is diagonal, $\boldsymbol{\Lambda}^{-1}$ is just one divided by each value along the diagonal.

Determinant

To compute the determinant of $\boldsymbol{\Sigma}$, SciPy uses the following mathematical fact: for any $D \times D$ matrix $\mathbf{M}$ with eigenvalues $\lambda_1, \lambda_2, \dots \lambda_D$, the determinant is the product of the eigenvalues or

$\det(\mathbf{M}) = \prod_{d=1}^{D} \lambda_d.$

However, computing the determinant this way is numerically unstable, as I’ve written about before. The upshot is that if $\lambda_n$ is small, the computed determinant might be zero due to machine precision. So SciPy computes the log of the PDF so that computing the determinant amounts to

$\log(\det(\mathbf{M})) = \sum_{d=1}^{D} \log(\lambda_d).$

To compute the PDF, SciPy first computes the log PDF and then computes the exponent of that quantity. For completeness, the log PDF for the multivariate normal is

$\begin{aligned} \log f(\mathbf{x}; \boldsymbol{\mu}, \boldsymbol{\Sigma}) &= -\frac{1}{2} (\mathbf{x} - \boldsymbol{\mu})^{\top} \boldsymbol{\Sigma}^{-1} (\mathbf{x} - \boldsymbol{\mu}) - \log\big(\sqrt{(2 \pi)^D \det(\boldsymbol{\Sigma})}\big) \\ &= -\frac{1}{2} \big[ (\mathbf{x} - \boldsymbol{\mu})^{\top} \boldsymbol{\Sigma}^{-1} (\mathbf{x} - \boldsymbol{\mu}) + \log\big((2 \pi)^D \det(\boldsymbol{\Sigma})\big) \big] \\ &= -\frac{1}{2} \big[ (\mathbf{x} - \boldsymbol{\mu})^{\top} \boldsymbol{\Sigma}^{-1} (\mathbf{x} - \boldsymbol{\mu}) + D\log(2 \pi) + \log(\det(\boldsymbol{\Sigma})) \big]. \end{aligned}$

Implementation

Below is an abbreviated version of SciPy’s implementation of multivariate_normal.pdf. For clarity, I’ve removed any code for modularity or standarization. SciPy’s implementation makes one additional optimization worth mentioning. Rather than computing

$(\mathbf{x} - \boldsymbol{\mu})^{\top} \boldsymbol{\Sigma}^{-1} (\mathbf{x} - \boldsymbol{\mu})$

naively, it computes $(\mathbf{x} - \boldsymbol{\mu})^{\top} \mathbf{Q} \sqrt{\boldsymbol{\Lambda}^{-1}}$ and then squares it.

def pdf(x, mean, cov):
    return np.exp(logpdf(x, mean, cov))


def logpdf(x, mean, cov):
    # `eigh` assumes the matrix is Hermitian.
    vals, vecs = np.linalg.eigh(cov)
    logdet     = np.sum(np.log(vals))
    valsinv    = np.array([1./v for v in vals])
    # `vecs` is R times D while `vals` is a R-vector where R is the matrix 
    # rank. The asterisk performs element-wise multiplication.
    U          = vecs * np.sqrt(valsinv)
    rank       = len(vals)
    dev        = x - mean
    # "maha" for "Mahalanobis distance".
    maha       = np.square(np.dot(dev, U)).sum()
    log2pi     = np.log(2 * np.pi)
    return -0.5 * (rank * log2pi + maha + logdet)

SciPy performs a bunch of other checks such as thresholding based on eigenvalues, raising an exception for singular matrices if specified, and so forth.