12
$\begingroup$

Suppose I have $n$ data points $x_1,\dots,x_n$, each of which is $p$-dimensional. Let $\Sigma$ be the (non-singular) population covariance of these samples. With respect to $\Sigma$, what is the most efficient way known to compute the vector of squared Mahalanobis distances (from $\vec 0$) of the n data points.

That is we would like to compute the vector $(x_1^T\Sigma^{-1}x_1,\dots,x_n^T\Sigma^{-1}x_n)$.

Computing the inverse $\Sigma^{-1}$ seems to be quite slow for large matrices. Is there a faster way?

$\endgroup$
1
  • $\begingroup$ Can you give reference for this information please? $\endgroup$
    – baleo
    Feb 24, 2021 at 12:20

1 Answer 1

17
$\begingroup$
  1. Let $x$ be one of your data points.

  2. Compute the Cholesky decomposition $\Sigma=LL^\top$.

  3. Define $y=L^{-1}x$.

  4. Compute $y$ by forward-substitution in $Ly=x$.

  5. The Mahalanobis distance to the origin is the squared euclidean norm of $y$:

$$ \begin{align} x^\top\Sigma^{-1}x &= x^\top(LL^\top)^{-1}x \\ &= x^\top(L^\top)^{-1}L^{-1}x \\ &= x^\top(L^{-1})^\top L^{-1}x \\ &= (L^{-1}x)^\top(L^{-1}x) \\ &= \|y\|^2. \end{align} $$

$\endgroup$
4
  • $\begingroup$ I had an issue with the accuracy of this proposed algorithm and created a new question: stats.stackexchange.com/questions/147654/… $\endgroup$ Apr 22, 2015 at 1:01
  • $\begingroup$ Without wanting to take away anything from how insight full this answer is (I up-voted earlier); especially because the OP (@Lepidopteris) uses MATLAB I would argue that for implementation purposes *only one should use x'* Sigma \ x where effectively Sigma \ x will take care of the inv(A) * x (please never use this later command). $\endgroup$
    – usεr11852
    Apr 22, 2015 at 2:06
  • $\begingroup$ Why do you say for implementation purposes one should use x'* Sigma \ x? Why is chol not a good option in MATLAB? $\endgroup$ Apr 22, 2015 at 2:47
  • $\begingroup$ @Lepidopterist: I have answered this comment in the comments of the new question. $\endgroup$
    – usεr11852
    May 7, 2015 at 1:56

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.