# Speeding up sampling of Gaussian entries

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}*)
RandomVariate[
NormalDistribution[], {d,
d}]; // AbsoluteTiming (* {0.083068,Null} *)



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}},
Module[{vals},
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 *)