2
$\begingroup$

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?

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

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]}}]

enter image description here

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:

enter image description here

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!

$\endgroup$
  • $\begingroup$ You might consider using LinearSolve[] instead of Inverse[] here. $\endgroup$ – J. M. will be back soon Apr 22 '17 at 11:01
  • $\begingroup$ Have an example? I tried to plug in LinearSolve[s,IdentityMatrix@Length@s] and had no luck. Much slower than Inverse[]. $\endgroup$ – Gregory Klopper Apr 22 '17 at 11:11
  • $\begingroup$ That's not the way to use LinearSolve[]; try Sqrt[(u - μ).LinearSolve[s, u - μ]]. You might want to see this. $\endgroup$ – J. M. will be back soon Apr 22 '17 at 11:18
  • $\begingroup$ 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. $\endgroup$ – Gregory Klopper Apr 22 '17 at 11:30
  • $\begingroup$ 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. $\endgroup$ – Gregory Klopper Apr 22 '17 at 11:41

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.