7
$\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$
12
$\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. Mahalanobis distance to the origin:

$$ \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$
  • $\begingroup$ I had an issue with the accuracy of this proposed algorithm and created a new question: stats.stackexchange.com/questions/147654/… $\endgroup$ – Lepidopterist Apr 22 '15 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 says Reinstate Monic Apr 22 '15 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$ – Lepidopterist Apr 22 '15 at 2:47
  • $\begingroup$ @Lepidopterist: I have answered this comment in the comments of the new question. $\endgroup$ – usεr11852 says Reinstate Monic May 7 '15 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.