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

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 ith location is the cell on the matrix that is in distance d(i) in direction tet(i) from the ith 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'.

1. I have changed the modulo conversion so I won't have to index and find the zero values after it.
2. 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:

• Any specific reason you convert q to uint32, and back to doubles within your function? Other than that I don't see any obvious weird things. You could pre-compute cos and sin, I'm sure you don't need very high precision. The power is also expensive to compute. Interpolating into tables for those could save a lot of time. Dec 31, 2017 at 2:19
• I have a hard time believing the MEX version is slower. How did you convert? Is there an automated conversion? A hand-written C function is likely to be faster, but not necessarily by much. Dec 31, 2017 at 2:21
• Thanks for the comment @CrisLuengo, I actually tried to create a table of distances to perform the power only once per simulation, but it turns out to be a little bit slower somehow. I'll try that for both power and the trigonometric functions, and see if it helps. As for MEX, I use the automatic code generation in Matlab, so it's probably less efficient than a native C code, but I don't write C regularly and have no time to make sure it will work as expected. However, I'm not surprised by that Mex taking longer, since I use only built-in functions, that are already written at a lower level...
– EBH
Dec 31, 2017 at 8:18
• ... and Mex impose a lot of unnecessary variable validations and communication between Matlab and C, which slowing things up. This is a well-known problem with using MEX (I can't find a good reference for that right now).
– EBH
Dec 31, 2017 at 8:23
• And, as for converting to uint32 - this is related to the outer program, where I use large vectors and save them, so in order to facilitate operations on them, I use uint32. This function is the only case where I need the double precision. Also, I find that the slowest part is usually generating random numbers (which is needed even if I prepare tables of distances before), so maybe there is a more efficient way to do that?
– EBH
Dec 31, 2017 at 8:31

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! :)

• Many thank for the answer, I'll take a closer look at the relevant parts, especially at part 1. As for your suggestions in parts 2&3, here are some relevant comments: 1. Modulo operation turns out to be slower also in Matlab; 2. You are right about the bias problem, I have already noticed it and changed my code properly, though, in a slightly different way, I'll try your correction for the modulo operation; ...
– EBH
Jan 12, 2018 at 7:05
• @EBH, I've edited part 2, I was thinking about a faster alternative to mod, I couldn't let go. :) Jan 12, 2018 at 7:21
• 3. The calling function use only linear indexing and store all the data in a vectorized from to speed things up. If I take your suggestion about coords I'll have to convert q to subscripts every time before calling the function and convert destinations back to linear. I had no time to continue working on this function lately, but a soon as I'll get back to it I'll update the question with further information.
– EBH
Jan 12, 2018 at 7:21
• The idea in part 1 looks like an interesting new direction, and I'll also test your new suggestion for fast_modulo. Stay tuned ;)
– EBH
Jan 12, 2018 at 7:25
• And thanks for the typo correction (corrected in question), I'm not a native English speaker so that isn't hard for me...
– EBH
Jan 12, 2018 at 7:32