a simple implementation of

Kernel Mahalanobis Distance

from

A.Ruiz and P.E. Lopez-de-Teruel, "Nonlinear Kernel-Based Pattern Analysis", IEEE Transactions on Neural Networks, vol. 12, no.1, pp.16-32.

aruiz@um.es
http://dis.um.es/~alberto

Definitions

Needs["Statistics`NormalDistribution`"]

RG[m_:  0, s_:  1] := Random[NormalDistribution[m, s]]

Meanval[l_] := Plus @@ (l/Length)[l]

Normaliza[l_] := (l/Plus) @@ l

SumaVectorial[x__List] := Plus[x]

Norm[x_List] := (Sqrt[x . x])

TakeDiagonal[m_] := 
      Table[m[[i, i]], {i, 
          Length[m]}]

Kernel["lin"] = #1 . #2 &;

\(\(Kernel["rbf", s_:  1] = 
        Function @@ {E\^\(-\(\(#1 . #1 + #2 . #2 - 2  #1 . #2\)\/\(2  s\^2\)\)\)};\)\)

Kernel["pol", p_:  2] = (((#1 . #2 + 1))^p) &;

KernelExpand[x_List, k_Function] := 
      Function @@ {Table[
            k[x[[i]], #], {i, 
              Length[x]}]}

\(SymDec[m_] := 
      Module[{w, α, r}, 
		{w, α} = 
          Eigensystem[m] // Chop; 
		r = 
          Count[w, _?Positive]; 
		{DiagonalMatrix@Take[w, r], 
          Take[α, r]}]\)

\(SpPow[λ_, p_] := 
      Module[{ts}, 
		ts = 
          FoldList[Plus, 0, Normaliza[λ]] // 
            Rest; 
		\(Position[ts, x_ /; x ≥ p] // First\) // 
          First]\)

\(SDCut[1 | 1. ] := Identity
\)
\(\(SDCut[p_]\)[{λ_, v_}] := 
      Module[{n, l}, 
		l = TakeDiagonal[λ]; 
		n = 
          SpPow[l, p]; 
		Print[n]; 
		{DiagonalMatrix[Take[l, n]], 
          Take[v, n]}]\)

\(KMahStr[x_, k_:  κ\_l] := 
      Module[{n, K, α, W, δ, m, u, na}, 
		n = 
          Length[x]; 
		K = Outer[k, x, x, 1]; 
		u\  = 
          Table[1/n, {i, n}, {j, n}]; 
		{W, α} = 
          SymDec[K - \ u . K - K . \ u + \ u . K . u] // 
            SDCut[$EigenPrec]; 
		δ = \(\@n\) 
            Inverse[W] . α; 
		m = Meanval[K]; 
		{k, 
          x, δ, m, {Det[W/n], W\_L}}]\)

dm[{k_, x_, δ_, m_,_}] :=
	Composition[Norm[δ . SumaVectorial[# - m]]^2  &, 
        KernelExpand[x, k]]

KernelMahalanobis[x_, k_] := KMahStr[x, k] // dm

Examples

dat = 
        Table[{10  Cos[t] + RG[], 10  Sin[t] + RG[]}, {t, 
            0, π,  .05}];
ListPlot[dat, PlotStyle -> AbsolutePointSize[5], 
      AspectRatio -> Automatic]

$EigenPrec = 1;
q = KernelMahalanobis[dat, Kernel["pol", 2]][{x, y}] // 
        Expand

Plot3D[q // Sqrt, {x, -20, 20}, {y, -10, 20}];

ContourPlot[q // Log, {x, -20, 20}, {y, -10, 20}, 
        ContourShading -> True, Contours -> 20, 
        AspectRatio -> Automatic];

$EigenPrec =  .98;
q = KernelMahalanobis[dat, Kernel["rbf", 5]][{x, 
            y}];

ContourPlot[q // Evaluate, {x, -20, 20}, {y, -10, 20}, 
        ContourShading -> True, Contours -> 10, 
        AspectRatio -> Automatic];

Plot3D[q // Evaluate, {x, -20, 20}, {y, -10, 20}, 
        PlotPoints -> 40, Mesh -> False, 
        AspectRatio -> Automatic];