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?

Answer 1

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!