# Numerical accuracy of multivariate normal distribution

In MATLAB, I've written two snippets of code that compute the PDF of a multivariate normal distribution. However there's a difference in the values these two methods produce and I can't figure out why. I've narrowed the problem down to something having to do with computing the inverse of the covariance matrix.

Inaccurate code

function p = mvnpdf_inacc(X, mu, sigma)
xc = bsxfun(@minus, X, mu);
[n, k] = size(xc);
twopic = (2 * pi) ^ (-k / 2);
sqrtdetsig = sqrt(det(sigma)) ^ -1;
c = twopic * sqrtdetsig;
p = zeros(n, 1);
for i = 1:n
xci = xc(i, :);
p(i) = c * exp(-0.5 * (xci / sigma * xci'));
end
end

Accurate code

function p = mvnpdf_acc(X, mu, sigma)
[R, err] = cholcov(sigma, 0);

if err
error('%s', 'sigma is not both symmetric and positive definite');
end

X0 = bsxfun(@minus, X, mu) / R;
d = min(size(X));
slogdet = sum(log(diag(R)));
p = exp(-0.5 * sum(X0 .^ 2, 2) - slogdet - 0.5 * d * log(2 * pi));
end

Testing code

function iseq_func = test_mvnpdf(n)
x = linspace(-2, 2, n);
y = x;
[X, Y] = meshgrid(x, y);
XY = [X(:), Y(:)];
mu = [0, 0];
sigma = [1.0, 0.5; 0.5, 0.4];
p_inacc = mvnpdf_inacc(XY, mu, sigma);
p_acc = mvnpdf_acc(XY, mu, sigma);
p_diff = abs(p_inacc - p_acc);
iseq_func = nnz(p_diff) == 0;
end

I get a value of false from running iseq_func(25). What the heck is going on here? Thanks!

• @varty because it is 'built-in' to the \$tati\$tic\$toolbox, which is another thousand dollars or so. – shabbychef Nov 17 '11 at 19:06 • Where is the multiplication by$-1/2$in the exponent in the first calculation? Am I just not seeing it? – cardinal Nov 17 '11 at 19:10 • I would pay attention to @cardinal's comment first. Roundoff error might account for differences beyond the 12th significant figure or so, and the product-vs-logarithmic approach could account for overflow problems when they occur, but otherwise neither can explain noticeable differences in the results. – whuber Nov 17 '11 at 20:46 • With @cardinal's correction (edited into the question it appears), I do not get any difference between the two versions. However, I had to change the cholcov to chol, because the former is in the statistics toolbox. (See mathworks.com/help/toolbox/stats/cholcov.html ). This may be the rub, however: I was testing the two implementations with a positive definite sigma. If that condition does not hold, then cholcov does something odd, which might explain inconsistencies, if they still exist. – shabbychef Nov 18 '11 at 1:15 • @cardinal Yes, the whole point of my question is that there is any difference whatsoever. I'm wondering why that is (other than just simply minor differences in floating-point calculations). max(p_diff) is something like 5e-14. Why should there be any difference at all? – Phillip Cloud Nov 21 '11 at 0:48 ## 1 Answer I am guessing you tested the code using a sigma which is not positive semidefinite. The 'accurate' implementation (which I am guessing is from the statistics toolbox) is computing$y^{\top}y$where$y = C^{-\top} \left(x - \mu\right)$where$C$is the output of cholcov on$\Sigma$(sigma). Note that$y^{\top}y$must be non-negative by design. If you feed in a$\Sigma$which is not PD, cholcov silently returns some version of$C$and the rest of the code proceeds (I would have made this an error, I think). The code you wrote, however, is computing$\left(x - \mu\right)^{\top}\Sigma^{-1}\left(x - \mu\right)$. If$\Sigma$is not PD, this quantity can be negative. When$\Sigma\$ is positive definite, cholcov returns a proper Cholesky factorization, and the results are the same (up to round off).

• I used the matrix sigma = [1.0, 0.5; 0.5, 0.4];. – Phillip Cloud Nov 20 '11 at 23:47
• My matrix is plenty symmetric. Run sigma == sigma' in MATLAB. – Phillip Cloud Nov 21 '11 at 0:52
• @cpcloud : this is not a problem of symmetry, it is whether the matrix is positive definite. Try, for example sigma = [1.0, -1.4;-1.4, 0.2]. chol bonks, but cholcov returns some matrix (I am guessing. I do not have the statistics toolbox). – shabbychef Nov 21 '11 at 17:46