40
$\begingroup$

I am interested in an efficient code to generate an $n$-D Gaussian random field (sometimes called processes in other fields of research), which has applications in cosmology.

Attempt

I wrote the following code:

fftIndgen[n_] := Flatten[{Range[0., n/2.], -Reverse[Range[1., n/2. - 1]]}]

Clear[GaussianRandomField];
GaussianRandomField::usage = "GaussianRandomField[size,dim,Pk] returns
a Gaussian random field of size size (default 256)  and dimensions dim 
(default 2) with a powerspectrum Pk";

GaussianRandomField[size_: 256, dim_: 2, Pk_: Function[k, k^-3]] := 
Module[{noise, amplitude, Pk1, Pk2, Pk3, Pk4},
Which[
dim == 1,Pk1[kx_] := 
If[kx == 0 , 0, Sqrt[Abs[Pk[kx]]]]; (*define sqrt powerspectra*)
noise = Fourier[RandomVariate[NormalDistribution[], {size}]]; (*generate white noise*)
amplitude = Map[Pk1, fftIndgen[size], 1]; (*amplitude for all frequels*)
InverseFourier[noise*amplitude], (*convolve and inverse fft*)
dim == 2,
Pk2[kx_, ky_] := If[kx == 0 && ky == 0, 0, Sqrt[Pk[Sqrt[kx^2 + ky^2]]]];
noise = Fourier[RandomVariate[NormalDistribution[], {size, size}]];
amplitude = Map[Pk2 @@ # &, Outer[List, fftIndgen[size], fftIndgen[size]], {2}];
InverseFourier[noise*amplitude],
dim > 2, "Not supported"]
]

Here are a couple of examples on how to use it in one and 2D

GaussianRandomField[1024, 1, #^(-1) &] // ListLinePlot
GaussianRandomField[] // GaussianFilter[#, 20] & // MatrixPlot

Question

The performance is not optimal — On other interpreted softwares, such Gaussian random fields can be generated ~20 times faster. Do you have ideas on how to speed things up/improve this code?

$\endgroup$
4
  • $\begingroup$ What is the most time-consuming part of your code? I see many functions that you use over and over again. What does putting in some well defined AbsoluteTiming and Print statements tell you? $\endgroup$
    – tkott
    Apr 27, 2012 at 18:06
  • $\begingroup$ size = 1024*1024; Pk1[k_] := If[k == 0 , 0, Sqrt[Abs[k^(-1)]]]; (*define sqrt powerspectra*); noise = Fourier[ RandomVariate[ NormalDistribution[], {size}]]; (*generate white noise*)// AbsoluteTiming amplitude = Map[Pk1, fftIndgen[size], 1]; (*amplitude for all frequels*)// AbsoluteTiming InverseFourier[noise*amplitude]; // AbsoluteTiming So I guess the Map Powerpsectra part ? {0.194598,Null} {4.388953,Null} {0.506372,Null} $\endgroup$
    – chris
    Apr 27, 2012 at 18:37
  • $\begingroup$ Great thread, as always ;) Too bad I wasn't around back then and have found it just now, but I guess it is better to work out something not being biased by such great answers :) $\endgroup$
    – Kuba
    Apr 8, 2014 at 18:27
  • 1
    $\begingroup$ @user26866 you could have a look at arxiv.org/pdf/0804.3536.pdf $\endgroup$
    – chris
    Jun 29, 2018 at 8:12

5 Answers 5

37
+100
$\begingroup$

Here's a reorganization of GaussianRandomField[] that works for any valid dimension, without the use of casework:

GaussianRandomField[size : (_Integer?Positive) : 256, dim : (_Integer?Positive) : 2,
                    Pk_: Function[k, k^-3]] := Module[{Pkn, fftIndgen, noise, amplitude, s2},
  Pkn = Compile[{{vec, _Real, 1}}, With[{nrm = Norm[vec]},
                                        If[nrm == 0, 0, Sqrt[Pk[nrm]]]], 
                CompilationOptions -> {"InlineExternalDefinitions" -> True}];
  s2 = Quotient[size, 2];
  fftIndgen = ArrayPad[Range[0, s2], {0, s2 - 1}, "ReflectedNegation"];
  noise = Fourier[RandomVariate[NormalDistribution[], ConstantArray[size, dim]]]; 
  amplitude = Outer[Pkn[{##}] &, Sequence @@ ConstantArray[N @ fftIndgen, dim]];
  InverseFourier[noise * amplitude]]

Test it out:

BlockRandom[SeedRandom[42, Method -> "Legacy"]; (* for reproducibility *)
            MatrixPlot[GaussianRandomField[]]
        ]

matrix plot of 2D Gaussian random field

BlockRandom[SeedRandom[42, Method -> "Legacy"];
            ListContourPlot3D[GaussianRandomField[16, 3] // Chop, Mesh -> False]
            ]

contour plot of 3D Gaussian random field

Here's an example the routines in the other answers can't do:

AbsoluteTiming[GaussianRandomField[16, 5];] (* five dimensions! *)
   {28.000959, Null}
$\endgroup$
3
  • $\begingroup$ Great answer! I only have one question: when you input size=256, what does it correspond in the physical space? Basically my question amounts to: what is the "physical" size of the system? Hope it is clear what I mean, thanks! $\endgroup$
    – lorenzop
    Dec 17, 2021 at 16:37
  • 1
    $\begingroup$ @lorenzo, the size argument just specifies the number of points taken in each dimension; for instance, GaussianRandomField[256, 2] should return a list with dimensions {256, 256}. $\endgroup$ Dec 17, 2021 at 16:41
  • $\begingroup$ OK thanks, so it is simply the field "resolution" and the "size" parameter won't affect the physical dimensions of your field. Then, probably my question changes into: what is the physical size of the system considered in your example? Specifically, how can I map the returned list (field) to actual space coordinates? $\endgroup$
    – lorenzop
    Dec 17, 2021 at 16:49
24
$\begingroup$

If the power spectrum is always of the form $1/k^p$ there are some optimisations you can make. The first is to compile a function to generate the reciprocal of the spectrum: $k^p$. This can be very fast since you don't need to check for the $k=0$ point (the reciprocal spectrum is well behaved at $k=0\ $).

We can then use ReplacePart to temporarily set the value at $k=0$ to something other than 0 (I've used 1), then take the reciprocal to get the spectrum we really wanted, and finally use ReplacePart again to set the value at $k=0$ back to 0.

The other optimisation is to create the random noise directly in frequency space, by having independent normally distributed real and imaginary components. Doing this allows you to get 2 fields at once, from the real and imaginary parts of the complex field which is generated.

The code below is only for 2D, but I think it should generalise to 1D or 3D quite straighforwardly.

The timing to generate 2 fields at 1024 x 1024 is about 280 ms on my PC.

fastpower=Compile[{{n,_Integer},{p,_Real}},
Block[{temp},
temp=Range[-n,n-2,2]^2;
Outer[Plus,temp,temp]^(p/2)]];

spectrum[n_,gamma_]:=ReplacePart[1/ReplacePart[
RotateRight[fastpower[n,-0.5gamma],{n/2,n/2}],
{1,1}->1.],{1,1}->0.];

complexrandoms[n_]:=RandomReal[NormalDistribution[0,1],{n,n}] 
+ I RandomReal[NormalDistribution[0,1],{n,n}]

fields[n_,gamma_]:={Re[#],Im[#]}&@Fourier[complexrandoms[n]spectrum[n,gamma]];

Example:

{f1,f2}=fields[1024,-3];

GraphicsRow[Image@Rescale@#&/@{f1,f2}]

enter image description here

$\endgroup$
5
  • $\begingroup$ Whaow! This is much faster! This could be generalized to any monotonic powerspectrum, couldn't it? $\endgroup$
    – chris
    May 1, 2012 at 8:07
  • $\begingroup$ I tried in 3D: fastpower = Compile[{{n, _Integer}, {p, _Real}}, Block[{temp}, temp = Range[-n, n - 2, 2]^2; Outer[Plus, temp, temp, temp]^(p/2)]]; spectrum[n_, gamma_] := ReplacePart[ 1/ReplacePart[ RotateRight[ fastpower[n, -0.5 gamma], {n/2, n/2, n/2}], {1, 1, 1} -> 1.], {1, 1, 1} -> 0.]; complexrandoms[n_] := RandomReal[NormalDistribution[0, 1], {n, n, n}] (Cos[#] + I Sin[#]) &[RandomReal[{0, 2 \[Pi]}, {n, n, n}]]; It takes 5 sec on 2 256^3 fields $\endgroup$
    – chris
    May 1, 2012 at 8:17
  • $\begingroup$ ... and 6 secs on 2 64^4 fields. This is truly efficient! $\endgroup$
    – chris
    May 1, 2012 at 8:30
  • $\begingroup$ Simon, I have a concern; it seems your field is not gaussian as demonstrated by e.g. {f1, f2} = fields[64*2, -3]; f1 // Flatten // Histogram $\endgroup$
    – chris
    May 3, 2012 at 21:21
  • $\begingroup$ @Chris, oops sorry. The complexrandoms need to have real and imaginary parts which are independently normally distributed (so the amplitude is Rayleigh distributed). I've edited my answer appropriately. Also, you probably need to look at several realisations of the field to see the statistics properly. Try Histogram@Flatten@Table[fields[64, -3], {20}] $\endgroup$ May 4, 2012 at 9:53
14
$\begingroup$

Here is a version of the code which incorporates the compilation suggestion of acl.

fftIndgen[n_] :=  Flatten[{Range[0., n/2.], -Reverse[Range[1., n/2. - 1]]}]

Clear[GaussianRandomField];
GaussianRandomField::usage = " GaussianRandomField[size,dim,Pk] 
returns a Gaussian random field of size size (default 256)
and dimensions dim (default 2) with a powerspectrum Pk default k^(-3))";
GaussianRandomField[size_: 256, dim_: 2, Pk_: Function[k, k^-3]] := 
Module[{noise, amplitude, Pk1, Pk2, Pk3,Pk4,kx,ky,kz,ku},
Which[
dim == 1,
Pk1 = Compile[{{kx, _Real}}, 
 If[kx == 0 , 0, Sqrt[Abs[Pk[kx]]]]]; (*define sqrt powerspectra*)
 noise = Fourier[RandomVariate[NormalDistribution[], {size}]]; (*generate white noise*)
 amplitude = Map[Pk1, fftIndgen[size], 1]; (*amplitude for all frequels*)
InverseFourier[noise*amplitude]//Chop, (*convolve and inverse fft*)
dim == 2,
Pk2 = Compile[{{kx, _Real}, {ky, _Real}},If[kx == 0 && ky == 0, 0,
  Sqrt[Pk[Sqrt[kx^2 + ky^2]]]]];
noise = Fourier[RandomVariate[NormalDistribution[], {size, size}]];
amplitude = Map[Pk2 @@ # &, Outer[List, fftIndgen[size], fftIndgen[size]], {2}];
InverseFourier[noise*amplitude]//Chop,
dim == 3,
Pk3 = Compile[{{kx, _Real}, {ky, _Real}, {kz, _Real}},
 If[kx == 0 && ky == 0 && kz == 0, 0, Sqrt[Pk[Sqrt[kx^2 + ky^2 + kz^2]]]]];
noise = 
Fourier[RandomVariate[NormalDistribution[], {size, size, size}]];
amplitude = Map[Pk3 @@ # &, 
 Outer[List,fftIndgen[size],fftIndgen[size],fftIndgen[size]], {3}];
InverseFourier[noise*amplitude]//Chop,
dim == 4,
Pk4 = Compile[{{kx, _Real}, {ky, _Real}, {kz, _Real}, {ku,_Real}},If[
kx == 0 && ky == 0 && kz == 0 && ku == 0, 0,  Sqrt[Pk[Sqrt[kx^2 + ky^2 + kz^2 + ku^2]]]]];
noise = 
Fourier[RandomVariate[NormalDistribution[], {size, size, size, size}]];
amplitude = 
Map[Pk4 @@ # &, 
 Outer[List, fftIndgen[size], fftIndgen[size], fftIndgen[size], 
  fftIndgen[size]], {4}];
InverseFourier[noise*amplitude]//Chop,
dim > 4, "Not supported"]
]

We can check that the GRF has the correct powerspectrum via the command

lm = {Range[128], 
  Map[Mean, 
    Table[Fourier[GaussianRandomField[256, 1, #^(-3/2) &]] // 
        Abs // #^2 &, {5}] // Transpose] // 
   Take[#, {2, 256/2 + 1}] &} // Transpose // Log // 
LinearModelFit[#, {1, x}, x] &;
lm["ParameterTable"]

An example

GaussianRandomField[] // MatrixPlot

a 2D GRF

It can also be used to generate cloud-like images

u = GaussianRandomField[] // GaussianFilter[#, 1] &;Image[u] // ImageAdjust

a 2D GRF

Finally a 3D Example

GaussianRandomField[16, 3] // Chop // ListContourPlot3D

a 2D GRF

In terms of performance we have

GaussianRandomField[128, 3]; // AbsoluteTiming
{7.169987,Null}
$\endgroup$
12
$\begingroup$

While I did not have time to do any major reorgnization, let me demonstrate one cheap way to speed this up, namely, compilation. I'll do it for the 1d case for simplicity.

fftIndgen[n_] := Flatten[{Range[0., n/2.], -Reverse[Range[1., n/2. - 1]]}]

grf1[size_,pk_] :=
    Module[ {noise,amplitude,Pk1},
    (*Pk1[kx_] :=If[ kx == 0,0,Sqrt[Abs[pk[kx]]]]; (*define sqrt powerspectra*)*)
        With[ {pkt = pk},
            Pk1 = Compile[{{kx,_Real}},
                If[ kx == 0,
                    0,
                    Sqrt[Abs[pkt[kx]]]
                ]
            ]
        ];
        noise = Fourier[RandomVariate[NormalDistribution[], {size}]]; (*generate white noise*)
        amplitude = Map[Pk1, fftIndgen[size], 1]; (*amplitude for all frequels*)
        InverseFourier[noise*amplitude]
    ]

Assuming I did not break anything, this is 5-6 times faster:

GaussianRandomField[2^19, 1, #^(-1) &]; // AbsoluteTiming
grf1[2^19, #^(-1) &]; // AbsoluteTiming
(*{1.755694, Null}
{0.281733, Null}*)

All I did here is compile the Pk1 function (the use of With is to inject the argument into the body to be compiled; maybe it's not necessary, but I did not check).

$\endgroup$
1
  • $\begingroup$ Uhm... I guess I should have thought about that one... Thanks ! and Thanks to J. M. for the Outer Trick. $\endgroup$
    – chris
    Apr 27, 2012 at 18:43
11
$\begingroup$

I know, it's long a ago. I just stumbled into this post and cannot resist. JM's code can be sped up even more, even without Compile: First, we can employ Outer in conjunction with Plus to compute the squared norms of frequencies at once. Second, only the zero frequency makes problems when applying Pk, so we can set its squared norm temporarily to 1. and use vectorization when we apply Pk. (Edit: I just realize that Simon Woods had the same idea...) That leads to a roughly 18-fold speed up on my machine for the whole routine (with size = 1024). (Another curiosity is that applying Sqrt twice is faster than applying Power[#,1/4]& once - probably due to hardware optimization.)

GaussianRandomField2[
  size : (_Integer?Positive) : 256, dim : (_Integer?Positive) : 2, 
  Pk_: Function[k, k^-3]] := Module[{Pkn, fftIndgen, noise, amplitude, s2},
  s2 = Quotient[size, 2];
  fftIndgen = N@ArrayPad[Range[0, s2], {0, s2 - 1}, "ReflectedNegation"];
  amplitude = Sqrt[Outer[Plus, Sequence @@ ConstantArray[fftIndgen^2, dim], dim]];
  amplitude[[Sequence @@ ConstantArray[1, dim]]] = 1.;
  amplitude = Pk[amplitude];
  amplitude[[Sequence @@ ConstantArray[1, dim]]] = 0.;
  noise = Fourier[RandomVariate[NormalDistribution[], ConstantArray[size, dim]]];
  Re[InverseFourier[noise amplitude]]
  ]

SeedRandom[666];
a = GaussianRandomField[1024]; // AbsoluteTiming
SeedRandom[666];
b = GaussianRandomField2[1024]; // AbsoluteTiming
Max[Abs[a - b]]

{1.29018, Null}

{0.069911, Null}

2.72278*10^-18

$\endgroup$
6
  • 1
    $\begingroup$ Well-done! $\phantom{}$ $\endgroup$ Feb 26, 2018 at 4:59
  • $\begingroup$ @J.M.! Good to see you again! $\endgroup$ Feb 26, 2018 at 6:40
  • $\begingroup$ I am completely confused. Your example works if the power spectrum is k^-3 with the Sqrt. But if I give it as an argument Function[k,Exp[-k^2]] it fails. $\endgroup$
    – chris
    Apr 25, 2020 at 14:34
  • $\begingroup$ as in SeedRandom[666]; R = 1/10; a = GaussianRandomField[1024*8, 1, Function[k, Exp[-R^2 k^2]]]; // AbsoluteTiming SeedRandom[666]; b = GaussianRandomField2[1024*8, 1, Function[k, Exp[-R^2 k^2]]]; // AbsoluteTiming Max[Abs[a - b]] $\endgroup$
    – chris
    Apr 25, 2020 at 14:39
  • $\begingroup$ Me too. It is a long time ago and I have other things to do at the moment. Also, after reading the post once more, I think that there was a point to apply Sqrt twice, but I cannot recall what it was... $\endgroup$ Apr 25, 2020 at 14:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.