Mahalanobis distance between two vectors

Question

I want to calculate the Mahalanobis distance between two vectors that represent two points.

For example:

u={1,2,4}; v={0,1,-2};

Mahalanobis[u_, v_] := Module[{cov, d}, (
   cov = Covariance[{u, v}];
   N@Sqrt[(u - v).PseudoInverse[cov].(u - v)]
   )]

I developed this function but I am not sure about the covariance matrix?

What seems to be the problem? Is the output not what you expect? — Mr.Wizard, Oct 14 '14 at 20:13
That's not what the Mahalanobis distance is. The covariance matrix should be a third parameter independent of u and v. — , Oct 14 '14 at 20:19
@RahulNarain, I don't understand what do you mean by your comment. — BetterEnglish, Oct 14 '14 at 22:47
For another explanation of Mahalanobis distance, see https://stats.stackexchange.com/questions/62092/bottom-to-top-explanation-of-the-mahalanobis-distance — Michael E2, May 22 '17 at 13:22
Example in the docs: http://reference.wolfram.com/language/ref/TTest.html#75259751 — Michael E2, May 22 '17 at 13:29

Gregory Klopper · Answer 1 · 2017-04-22T11:40:47.533

This is not what a Mahalanobis distance is. It isn't a distance between 2 vectors. It is defined as a distance between a vector and a cohort of vectors with a given mean and a covariance matrix (of the cohort).

Try this instead:

Dm = Compile[{{u, _Real, 1}, {\[Mu], _Real, 1}, {s, _Real, 2}}, 
  First@\[Sqrt]((u - \[Mu]).Inverse[s].Transpose[{u - \[Mu]}]), 
  CompilationOptions->{"ExpressionOptimization"->True}, 
  RuntimeOptions->"Quality", RuntimeAttributes->Listable, CompilationTarget->"C"]

cohort = RandomVariate[BinormalDistribution[{5, 5}, {.5, 1.5}, .3], 1000];
\[Mu] = Mean@cohort; s = Covariance@cohort;
Print["\[Mu] = ",{\[Mu]}\[Transpose]//MatrixForm, "   S = ",s//MatrixForm];

points = {\[Mu]}
  ~Join~Table[N@5+2{Cos[x],Sin[x]}, {x,1/16\[Pi],2\[Pi],1/8\[Pi]}]
  ~Join~Table[N@5+4{Cos[x],Sin[x]}, {x,1/16\[Pi],2\[Pi],1/8\[Pi]}];

ListPlot[{cohort, Labeled[#,Round[Dm[#,\[Mu],s],.01]]&/@points},
  PlotRange->{{0,10},{0,10}},AspectRatio->1,PlotStyle->{Darker@LightBlue,{Red,PointSize[.01]}}]

A bit of optimization:

(* Inverse[] cannot be parallelized and takes too long *)
manypoints=RandomVariate[NormalDistribution[5,3],{1000000,2}];
AbsoluteTiming[Dm[#,\[Mu],s]&[manypoints];]
(* {5.973094, Null} *)

(* using 2D matrix inverse formula as a special case *)
FastDm2D=Compile[{{u,_Real,1},{\[Mu],_Real,1},{s,_Real,2}},
  First@Sqrt[(u-\[Mu]).{{s[[2, 2]]/(-s[[1, 2]] s[[2, 1]] + s[[1, 1]] s[[2, 2]]), -(
  s[[1, 2]]/(-s[[1, 2]] s[[2, 1]] + s[[1, 1]] s[[2, 2]]))}, {-(
  s[[2, 1]]/(-s[[1, 2]] s[[2, 1]] + s[[1, 1]] s[[2, 2]])), 
  s[[1, 1]]/(-s[[1, 2]] s[[2, 1]] + s[[1, 1]] s[[2, 2]])}}.Transpose[{u-\[Mu]}]],
CompilationOptions->{"ExpressionOptimization" -> True},Parallelization->True,
RuntimeOptions->"Quality",RuntimeAttributes->Listable,CompilationTarget->"C"];
AbsoluteTiming[FastDm2D[#,\[Mu],s]&[manypoints];]
(* {0.222699, Null} *)

Of course that looks MUCH prettier when typed in Mathematica:

EDIT - OPTIMIZATION:

(* adding a wrapper to precompute inverse of S produces the fastest results *)
FastDmCompiled = 
  Compile[{{u, _Real, 1}, {\[Mu], _Real, 1}, {sInv, _Real, 2}},
   First@Sqrt[(u - \[Mu]).sInv.Transpose[{u - \[Mu]}]],
   CompilationOptions -> {"ExpressionOptimization" -> True}, 
   Parallelization -> True, RuntimeOptions -> "Quality",
   RuntimeAttributes -> Listable, CompilationTarget -> "C"];
FastDm[u_, \[Mu]_, s_] := FastDmCompiled[u, \[Mu], Inverse[s]];

AbsoluteTiming[FastDm[#,\[Mu],s] &[manypoints];]

(* {0.151167,Null} *)

Hope that helps. Good luck!

You might consider using LinearSolve[] instead of Inverse[] here. — J. M.'s missing motivation, Apr 22 '17 at 11:01
Have an example? I tried to plug in LinearSolve[s,IdentityMatrix@Length@s] and had no luck. Much slower than Inverse[]. — Gregory Klopper, Apr 22 '17 at 11:11
That's not the way to use LinearSolve[]; try Sqrt[(u - μ).LinearSolve[s, u - μ]]. You might want to see this. — J. M.'s missing motivation, Apr 22 '17 at 11:18
Still slower on my machine than Inverse, and still cannot be parallelized. I'm running 11.1 on Win7 with 8-thread i7 CPU. It's probably because it is already parallelized internally, and thus refuses to be run in parallel inside a compiled function. Says CompiledFunction::pext: Instruction 3 in CompiledFunction[...] calls ordinary code that can be evaluated on only one thread at a time. — Gregory Klopper, Apr 22 '17 at 11:30
Of course Inverse[s] can simply be passed into the function. ;-) That solution simply FLIES! for any vector dimension. Just edited the post to append this result. — Gregory Klopper, Apr 22 '17 at 11:41

Mahalanobis distance between two vectors

1 Answers1