The radial part of a normal distribution

Question

I am reading a paper that asks me to sample $s_i$ from a distribution like this:

$s_i \sim (2\pi)^{-\frac{d}{2}}A^{-1}_{d-1}r^{d-1}e^{-\frac{r^2}{2}}$

"Here the normalization constant $A_{d−1}$ denotes the surface volume of the d-dimensional unit ball. Thus $s_i$ matches the radial part of a normal distribution"

What is the radial part of a normal distribution?

The distribution above is related to the inverse of the hyper-surface of the unit ball, but the radius part is not inverted. Can i just sample from a normal distribution and scale? I am trying to get a sample from this distribution in C++.

Thanks in advance for any help in understanding the distribution above.

Answer 1

I gauss you are working on the fast-food method. As it is pointed out in the paper the distribution matches the radial part of a normal distribution. For $d = 2$ it is equal to Rayleigh distribution so you can use the standard implementations to sample $s_i$. For higher value of d you can sample from a $d$ dimensional zero mean gaussian distribution and then compute the length of the random vector. This method works well in terms of correctness. However, generating n samples using this method needs $nd$ time which is much higher than the overall time complexity of the fast-food algorithm. Thus, it is fine just for testing the correctness of the algorithm.

Update: I recently realized that $s_i$ has chi distribution. Equivalently, you can sample from chi-squared distribution and compute the square root of the samples. I succsesfully tested it in the Matlab and following you can find the source code:

function [ ws ] = FF_samples( ndim, nsamples, sigm)
    ws = [];
    d = 2 ^ ceil(log2(ndim));
    x = ceil(nsamples/d);
    for bi = 1:x
        I = eye(d);

        B = diag(2*randi(2, d, 1) - 3);
        H = hadamard(d);
        P = I(randperm(d), :);
        G = diag(normrnd(0,1, d, 1));
        S = diag( chi2rnd(d,d,1) .^ .5) * sum(diag(G) .^ 2).^-.5;
        %%
        nws = (1/sigm / d)^.5 * S * H * G * P * H * B;
        ws = [ws; nws];
    end
    ws = ws(1:nsamples, 1:ndim);
end

P.S.: It still needs some improvement since I don't use the fast Fourier transform like algorithm for multiplication of the Hamard matrix to a vector.