Re: efficiently adding many 2D Gaussians

• To: mathgroup at smc.vnet.net
• Subject: [mg91590] Re: efficiently adding many 2D Gaussians
• From: Jens-Peer Kuska <kuska at informatik.uni-leipzig.de>
• Date: Thu, 28 Aug 2008 03:18:49 -0400 (EDT)
• References: <g93b4e\$koc\$1@smc.vnet.net>

```Hi,

look at

http://www.umiacs.umd.edu/~yangcj/papers/siam_fgt_v11.pdf

or search "Fast Gauss transform" with Google

Regards
Jens

M.Roellig wrote:
> Dear group,
>
> I need some help on reducing run time and memory usage for the
> following problem:
>
> I want to place a large number N (10000+)  of 2D Gaussians in a plane
> (on a large array, m x m ) and add them up.
> I came up with a few approaches but they all fail, or become very slow
> once N >10000 or m>100.
> Here are my ideas so far (not bothering about dimensions, units, etc.
> and using 1000 positions only):
>
> SeedRandom[1111111];
> posTab =
>  RandomReal[NormalDistribution[0, 1], {1000, 2}];
>
> 1. idea (fastest I could get):
>
> gsum = Compile[{{x, _Real, 1}},
>    Evaluate[
>     Total[Map[(1.*
>          Exp[-(((x[[1]] - #[[1]])^2 + (x[[2]] - #[[2]])^2)/(
>             2. 0.25^2.))]) &, posTab]]]];
>
> (* Warning Output:
> -----------------------
> Part::partd: Part specification x[[1]] is longer than depth of \
> object. >>
>
> Part::partd: Part specification x[[2]] is longer than depth of \
> object. >>
> ------------------------*)
>
> xyArray = Array[{#1, #2} &, {100, 100}] // MatrixForm;
> map1e4 = Map[gsum, xyArray, {3}]; // Timing
>
> Out[22]= {11.843, Null}
>
>
>
> 2. idea (trying to use mapping, and came up with a very inefficient
> solution :-)
>
> Total[
>    Table[Exp[-(((x - #[[1]])^2 + (y - #[[2]])^2)/0.125)], {x, -5,
>        5, .1}, {y, -5, 5, .1}] & /@ posTab]; // Timing
>
> Out[9]= {85.375, Null}
>
> Using more than 1000 Gaussians doesn't fit into my memory and makes
> WinXp swapping memory to disk
> basically making the computation time unbearable. Of course this is no
> big surprise since I create one array for each individual Gaussian
>
>  3. completely different approach hoping to take advantage of internal
> mathematica optimization:
>
> SeedRandom[1111111]; posTab =
>  RandomReal[NormalDistribution[0, 100], {1000, 2}];
> clumpSum =
>   Total[Map[(1.*
>        Exp[-(((x - #[[1]])^2 + (y - #[[2]])^2)/(2. 10.^2.))]) &,
>     posTab]];
> clumpEm100 =
>    Plot3D[clumpSum, {x, -350, 350}, {y, -350, 350}, PlotRange -> All,
>     PlotPoints -> 100, MaxRecursion -> 0]; // Timing
>
> mapEm100 = clumpEm100[[1]] /. GraphicsComplex[pts_, __] -> pts;
>
>
> Out[25]= {53.485, Null}
>
> I am grateful for any advice on how to do this faster and more
> efficiently.
> I guess a clever usage of Fold/Nest could do it, but I failed to get
> it done.
>
> The code I use now is a little more complex, since I need to add up
> Gaussians with
> different sizes and amplitudes, which I take from a list:
>
> --------------------
>
> fullEns ={{1/1000, 63095, 374775., 0.00237989, 112.101, 0.818147},
> {1/100,
>   10000, 185962., 0.00647649, 20.143, 2.22645}, {1/10, 1584, 92273.2,
>   0.0176247, 5.91138, 6.05894}, {1, 251, 45785.5, 0.0479629, 0.837552,
>    16.4884}, {10, 39, 22718.5, 0.130523, 0.156455, 44.8706}, {100, 6,
>   11272.8, 0.355198, 0.220978, 122.108}, {1000, 1, 5593.51, 0.966614,
>   0.240513, 332.297}}
>
> (*In each sublist element 2 is the number, element 5 the amplitude,
> element 6 the stand.dev.
> I am discarding the first to rows to reduce total numbers. *)
>
> SeedRandom[1111111]; posEnsTab =
>  Append[Table[
>    RandomReal[NormalDistribution[0, 100], {fullEns[[j, 2]], 2}], {j,
>     1, Length[fullEns] - 1}], {{0., 0.}}];
> ensSum = Plus @@
>    Table[Total[
>      Map[( fullEns[[j, 5]]*
>          Exp[-(((x - #[[1]])^2 + (y - #[[2]])^2)/(
>             2. fullEns[[j, 6]]^2.))]) &, posEnsTab[[j]]]], {j, 3,
>      Length[fullEns]}];
> ensSum100 =
>    Plot3D[ensSum, {x, -350, 350}, {y, -350, 350}, PlotRange -> All,
>     PlotPoints -> 100, MaxRecursion -> 0]; // Timing
>
> (* Out: {116.453, Null} *)
>
> ensMap100 = ensSum100[[1]] /. GraphicsComplex[pts_, __] -> pts;
> ListPlot3D[ensMap100, PlotRange -> All]
>
> Hope somebody can help me.
>
> Cheers,
>
> Markus
>

```

• Prev by Date: Re: efficiently adding many 2D Gaussians
• Next by Date: Re: NotebookWrite Output Cells That Will Delete from Menu?