5
$\begingroup$

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!

$\endgroup$
  • 1
    $\begingroup$ @varty because it is 'built-in' to the \$tati\$tic\$ toolbox, which is another thousand dollars or so. $\endgroup$ – shabbychef Nov 17 '11 at 19:06
  • 4
    $\begingroup$ Where is the multiplication by $-1/2$ in the exponent in the first calculation? Am I just not seeing it? $\endgroup$ – cardinal Nov 17 '11 at 19:10
  • 2
    $\begingroup$ 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. $\endgroup$ – whuber Nov 17 '11 at 20:46
  • 1
    $\begingroup$ 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. $\endgroup$ – shabbychef Nov 18 '11 at 1:15
  • 1
    $\begingroup$ @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? $\endgroup$ – Phillip Cloud Nov 21 '11 at 0:48
1
$\begingroup$

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).

$\endgroup$
  • $\begingroup$ I used the matrix sigma = [1.0, 0.5; 0.5, 0.4];. $\endgroup$ – Phillip Cloud Nov 20 '11 at 23:47
  • $\begingroup$ My matrix is plenty symmetric. Run sigma == sigma' in MATLAB. $\endgroup$ – Phillip Cloud Nov 21 '11 at 0:52
  • $\begingroup$ @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). $\endgroup$ – shabbychef Nov 21 '11 at 17:46

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.