# Slow application of Outer on a matrix

I have some performance issue with the calculation of a gaussian kernel applied to a not-that-large dataset. I'm not very good with performance tuning in Mma and I have the feeling this one could be made much faster.

Here are the methods I use :

gaussianKernel[sigma_] := Function[{x1, x2},
N[Exp[-(x1 - x2).(x1 - x2)/(2 N[sigma]^2)]]
];

K[sigma_] := With[{kernel = Outer[gaussianKernel[sigma], X, X, 1]},
ArrayFlatten[{{1, kernel}}] (* add bias column to K *)
];

X // Dimensions; (* {194,32} *)

Timing[K[1.]] (*{0.761152, ...}*)


So quite slow... Just to be clear, K is a matrix where $K_{i,j}$ = gaussianKernel[X[[i]],X[[j]]].

Plus I need to iterate over the argument of K to find the optimal value. Any way this could be sped up? There's one way I can see, but I'm afraid it would be quite hard to implement : since K is symmetric, I could theoretically compute only half of these values.

Thanks for your time !

## 2 Answers

X = RandomReal[1, {194, 32}];

With[{arr = ArrayFlatten[{{0.0, -Total[Outer[Plus, X, -X, 1]^2, {3}]}}]},
k[sigma_] := Exp[arr/(2 sigma^2)]]

AbsoluteTiming[Do[k[s], {s, 1, 1000}]]
(* {0.257026, Null} *)


Create the matrix once, substitute sigma...

(* fake some data *)
X = N@RandomInteger[{1, 10}, {194, 32}];

array = ArrayFlatten[{{1, E^(1/(2 sigma^2)
Map[Tr, -(Outer[Subtract, X, X, 1]^2), {2}])}}]; //Timing

(* do new for sigma 1 to 20 *)
sigs = Table[Replace[array, sigma -> N@sig, {5}], {sig, 1, 20}]; // Timing

(* op method *)
(* do old for sigma 1 to 20 *)
old = Table[K[N@sig], {sig, 1, 20}]; // Timing

old==sigs

(*

{1.716011,Null}

{7.909251,Null}

{75.114481,Null}

True

*)


Netbook timings, but ~ 9X faster. Probably a more sophisticated mathematical means of doing this, don't have time now, will ponder.

Update: Moving some work out allowing direct operations vs substitutions:

siger[sigma_, basearray_] := ArrayFlatten[{{1, E^((1/(2 sigma^2)) basearray)}}];

basearray = Map[Tr, -(Outer[Subtract, X, X, 1]^2), {2}]; // Timing

newsigs = Table[siger[sig, basearray], {sig, 1, 20}]; // Timing

newsigs == sigs == old

(*

{1.404009,Null}

{0.499203,Null}

True

*)


So ~150X faster (again, on netbook...)