Generating locations based on a power-law distribution of distance from their origin

Question

I have a small function that gets a vector of N linear indices locations in a matrix and returns N new locations in the same matrix. For each of the input locations, a random angle tet and a distance d from a power-law distribution are randomized. The new i^th location is the cell on the matrix that is in distance d(i) in direction tet(i) from the i^th input location. The matrix represents a torus map, so if the new location is out of the matrix boundaries, it flows back from the other side.

This function is a part of a much larger program that simulates ecological community. It is called hundreds of millions of times and it is the slowest part (by far) in the program, taking about 10 times more than the next slowest operation.

I have tried the following:

1. Convert it to a MEX file
2. Casting all variables to the minimum size (i.e. uintX instead of double)
3. Using the built-in GeneralizedParetoDistribution to make a PDF object for the randomization

None of them reduced execution time, all of them did the opposite. Also, I cannot use parallel computing as I already use it at a 'higher' level of the program.
I'll be glad for any advice/idea that will make this run faster.
Thanks.

EDIT

Thanks to several advise from @Cris Luengo I was able to cut some runtime:

1. I use another random generator called 'simdTwister' which is apparently significantly faster than the default 'twister'.

2. I have changed the modulo conversion so I won't have to index and find the zero values after it.
3. I use a predefined list of values from power-law distribution to avoid the power computation on every run.

All other ideas within the comments and the answer was slower than the original code.

This may look like an answer to the question, but it isn't. The reason I'm still puzzled about it is that when I implement the third part (using a list of distances) in my outer function and use it in within (part) of the full program, it suddenly takes longer with this pre-computed list then without it, although I compute the list only once for all runs. I have checked this in different ways and have no idea why it is so. Please let me know if there is something wrong (or inefficient) with the way I implement part 3 in the code below.
Also, I would like to know if there is a way to randomize numbers from a 2D Pareto distribution, so I can avoid randomizing angles and converting them.

I have updated the code below to reflect this changes.

CODE

Here is the function code:

function [destinations] = randisp_so(Origins,L,alpha,distances)
% select random direction with a distance from power-law distribution and
% convert back to vectorized form of the matrix
%
% ORIGINS is a list of origin cells
% L is the size of one side of the map
% ALPHA is the scale parameter of the distribution
% The destinations vector is the same size as ORIGINS

OriginSize = numel(Origins);
maxd = floor(sqrt(2)*L/2); % maximum distance for dispersal

% converting vetorized represatation to matrix:
c = ceil(double(Origins)./L); % ind2sub(q) - columns
r = double(Origins)-(c-1).*L; % ind2sub(q) - rows

if nargin==4
    % choose random length from a given distribution (as array):
    distance = maxd.*distances(randi(numel(distances),OriginSize,1));
else
    % List of distances to disperse from power-law distribution:
    e = (-1./(alpha-1));
    RND = (1-rand(OriginSize,1));
    distance = maxd.*RND.^e;
end

% choose a random angle from uniform distribution [0 2*pi]
angles = rand(OriginSize,1)*2*pi;

new_r = round(r + distance.*cos(angles)); % convert length to cell row
new_c = round(c + distance.*sin(angles)); % convert length to cell column

% transform rows and columns to a 'close-map'
new_r = mod(new_r-1,L)+1;
new_c = mod(new_c-1,L)+1;

% convert back to vectotized form:
destinations = new_r+(new_c-1).*L; %sub2ind(r,c);
end

And a little testbench for running this function:

% randisp_so testbench
n = 10;
L = 2^n+1; % the size of one side of the map
MapSize = L^2; % total number of cells in the map
N = MapSize*30;
alpha = 0.7;
rng('shuffle','simdTwister') % the fastest random generator
q = uint32(randi(MapSize,N,1));
% List of distances to disperse from power-law distribution:
distances = linspace(1,0,MapSize).^(-1./(alpha-1)).';

timeit(@() randisp_so(q,L,alpha)) % 2.5153 sec
timeit(@() randisp_so(q,L,alpha,distances)) % 1.9008 sec

And here is a comparison of the profile of a typical run:

Answer 1

Looking at your code again, there are a few things that you could try:

1. Avoid trigonometric functions

These seem to be the most expensive part of your code. You compute a single uniformly distributed random value, turn it into a random angle, then get coordinates on the circle from that:

angles = rand(qSize,1)*2*pi;
new_r = r + ceil(distance.*cos(angles)); % convert length to cell row
new_c = c + ceil(distance.*sin(angles)); % convert length to cell column

Instead you could compute two normally distributed values, and normalize the vector:

pos = randn(n,2);
new_coords = coords + ceil(distance(pos ./ sqrt(sum(pos.^2,2)));

(with coords would be equivalent to [c,r] for simplicity here).

The 2D normal distribution is isotropic, this leads to exactly the same probability distribution.

In my test this was only about 10% faster, but every bit counts, right? randn is more expensive than rand, and there's the sqrt as well.

Ideally there would be a way to convert the Gaussian distribution into the power law distribution. This would avoid normalizing pos, and then sampling another random value for the distance. I don't know if this is possible, or how costly it would be.

2. Modulo operation

I'm surprised by the cost of the modulo operation. But it makes sense, it's a floating-point operation. Your coordinates never extend past the end of the matrix by a whole lot. You could instead try comparisons:

% original code
new_r = mod(new_r,A)+1;

% alternative
I = new_r>A;
new_r(I) = new_r(I)-A;
I = new_r<1;
new_r(I) = new_r(I)+A;

This is actually slower in Octave, but it's possible that MATLAB does a better job at optimizing these indexing operations. Worth a try in any case. An alternative to mod is rem, but I don't think either is faster than the other.

Note that your original code is biased: you always add 1 to the index. This causes your movement to be biased by 1 to the right and bottom of the matrix. To be correct it would have to be:

new_r = mod(new_r-1,A)+1;

The other bias in your code is the use of ceil instead of round when you initially compute the coordinates. This causes your movement to be biased by an additional 0.5 to the right and bottom of the matrix. I recommend that you plot the output of randisp_so for a thousand points, all in the same position in the matrix. Make sure that the output positions are distributed isotropically around that point!

Edit:

This function is about 30% faster than mod (sorry for the poor variable names):

function v = fast_modulo(v,A)
persistent X;
if numel(X)~=3*A
   X = repmat(1:A,1,3);
end
v = X(v+A);

You'd call

new_r = fast_modulo(new_r,A);
new_c = fast_modulo(new_c,A);

The function expects integer inputs that are not too far outside of the 1:A range (should be OK in your application). As long as you call this function with the same A, it doesn't need to recompute X.

3. Prevent conversion linear index from coordinates

This is not the most expensive part of your code, but it seems illogical to me. You could pass coordinates to/from your function:

[destenations,coords] = randisp_so(coords,A,alpha)

The calling function would use destenations (I presume), and store coords for the next iteration. coords would be equivalent to [c,r].

BTW: did you notice that destenations is a typo? It's hard to type for me! :)