3
$\begingroup$

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 !

$\endgroup$
2
$\begingroup$
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} *)
$\endgroup$
1
$\begingroup$

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...)

$\endgroup$

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.