Fast computation of entries of large matrix

Hey guys,

I'm working on a problem with creating very large matrices (in the real problem the size is nearly 2500 x 2400). The computation can be seen below in the small example. Because I have to do this computation several thousand times, I'm wondering, if there is any faster way to get the same result, maybe by using parallel computing, gpu computing, etc.

I tried to vectorize the calculation, but I found no proper way, which was faster in the end.

I'm glad for any help :)

x1=rand(2500,5);   % some matrix with high number of rows and small number of columns > 1
x2=rand(2400,5);   % another matrix with nearly same size as x1
sigma=[1,2,3,4,5]; % some parameter which is either a scalar or a vector, 
                   % whose length is identical to the number of columns as x1 & x2


% prepare the resulting matrix (2500 x 2400)
X=zeros(length(x1(:,1)),length(x2(:,1))); % preallocation of storage

% calculate the entries:
tic
for i=1:length(x1(:,1))
    for j=1:length(x2(:,1))

        X(i,j)=exp( sum( -1./(sigma.^2) * ((x1(i,:)-x2(j,:)).^2)' ) ); 
        % sum() is necessary because of the changing length of sigma

    end
end
toc

Elapsed time is 7.122652 seconds.

[Cameron]

I have two initial questions:

1) Are you doing your second loop through the number of rows or columns? They way you wrote it initially makes your loop go through the rows twice. For example:

for j=1:length(x2(:,1))
end
%or
for j=1:size(x2,2)    
end

2) Should you be looping over your your variable sigma like this:

X(i,j)=exp( sum( -1./(sigma(j)^2) * ((x1(i,:)-x2(j,:)).^2)' ) ); 

Additionally, I would not name your variable "sigma" as that is already the name of a built in function.

[Jan]

A general hint: Avoid length(x1(:,1)) but use the more efficient and nicer size(x1, 1) or in modern Matlab versions height(x1).

[JoRo]

Thank you both for you hints.

I will rename the variable to avoid name problems :)

I also made minor changes in the example. Now x2 has nearly the same size (this can happen in the real calculation).

So the second for-loop is also going through the number of rows, but this time rows of x2.

You have to consider the rows as multidimensional vectors, while the columns store different data points.

[JoRo]

Wow, thank you both for your hints.

I have to go through your advices to completely understand what you did, but this seems great!

Thank you :)

Answer 1

[Matt J]

x1=rand(2500,5);   % some matrix with high number of rows and small number of columns > 1
x2=rand(2400,5);   % another matrix with nearly same size as x1
sigma=[1,2,3,4,5]; % some parameter which is either a scalar or a vector, 
                   % whose length is identical to the number of columns as x1 & x2

tic;                 
                   
 X= exp(-pdist2(x1./sigma,x2./sigma).^2);
 
toc

Elapsed time is 0.155956 seconds.

[Jan]

You can save the time for the squareroots also:

x1=rand(2500,5);   % some matrix with high number of rows and small number of columns > 1
x2=rand(2400,5);   % another matrix with nearly same size as x1
sigma=[1,2,3,4,5]; % some parameter which is either a scalar or a vector, 
                   % whose length is identical to the number of columns as x1 & x2
tic;
X = exp(-pdist2(x1./sigma,x2./sigma).^2); 
toc

Elapsed time is 0.080255 seconds.

tic;
Y = exp(-pdist2(x1./sigma,x2./sigma, 'squaredeuclidean')); 
toc

Elapsed time is 0.047662 seconds.

max(abs(X-Y), [], 'all')

ans = 1.1102e-16

Answer 2

[Jan]

n = 2500;
x1=rand(n, 5);     % some matrix with high number of rows and small number of columns > 1
x2=rand(n, 5);     % another matrix with same size as x1
sigma=[1,2,3,4,5]; % some parameter which is either a scalar or a vector, 
                   % whose length is identical to the number of columns as x1 & x2
X = zeros(n, n);   % preallocation of storage

tic
v = -1./(sigma.^2);
for j = 1:n
    w = x2(j, :);
    for i = 1:n
        X(i,j) = exp(v * ((x1(i, :) - w).^2).');
    end
end
toc

Elapsed time is 3.445791 seconds.

This runs with about the double speed. Then comment "sum() is necessary because of the changing length of sigma" is strange: sigma seems to be a constant?! Then the argument of the sum is a scalar in all cases, such that the sum can be omitted.

Matt J's pdist2 apporach ist much faster.