I'm trying to speed up sampling from the projected Gaussian distribution (like Gaussian, but vectors normalized to have norm 1). Code below generates $d$ samples of size $d$ and it's about 300x slower than sampling from standard normal for $d=2000$.

Any suggestions how to reduce this overhead?

d = 2000;

normalizedGaussianSampler[diag_] := 
  With[{d = Length[diag]}, 
   Compile[{{n, _Integer}}, 
    vals = Sqrt[diag]*# & /@ 
      RandomVariate[NormalDistribution[], {n, d}];
    #/Sqrt[Total[#*#]] & /@ vals]];

sampler = normalizedGaussianSampler[Table[1/i, {i, 1, d}]];
sampler[d]; // AbsoluteTiming (*{21.170401`,Null}*)
   NormalDistribution[], {d, 
    d}]; // AbsoluteTiming (* {0.083068`,Null} *)


1 Answer 1


Take a look at the examples for Compile. Whenever you define some new symbols inside your function, you should make them locally scoped by using Module. Otherwise, the compiled function has to transfer the values (in your case vals) back to the main evaluator, which takes a lot of time.

normalizedGaussianSamplerFixed[diag_] := 
  With[{d = Length[diag]}, Compile[{{n, _Integer}},
     vals = 
      Sqrt[diag]*# & /@ RandomVariate[NormalDistribution[], {n, d}];
     #/Sqrt[Total[#*#]] & /@ vals]

samplerFixed = normalizedGaussianSamplerFixed[Table[1/i, {i, 1, d}]];
samplerFixed[d]; // AbsoluteTiming 
(* 0.165976 *)

RandomVariate[NormalDistribution[], {d, d}]; // AbsoluteTiming 
(* 0.135395 *)

Your Answer

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

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