Vectorize двойной для петель в Matlab

Question

Я представляю мой простой рабочий код Matlab и будет задавать вопросы:

tic 

nrand1 = 10000; 
nrand2 = 20000; 

% Location matrix 1: [longitude, latitude, w1] 
lmat1=[rand(nrand1,1)-75 rand(nrand1,1)+39 round(rand(nrand1,1)*1000)+1]; 

% Location matrix 2: [longitude, latitude, w2]  
lmat2=[rand(nrand2,1)-75 rand(nrand2,1)+39 round(rand(nrand2,1)*100)+1]; 

% The number of rows for each matrix = In fact it's nrand1 X nrand2, obviously 
nobs1 = size(lmat1,1); 
nobs2 = size(lmat2,1); 

% The number of pair-wise distances 
% between L1 locations X L2 locations 
ndist = nobs1*nobs2; 

% Initialization: Distance vector and weight vector 
hdist = zeros(ndist,1); 
weight = zeros(ndist,1); 

% Double for loop -- for calculating the pair-wise distances and weights 
k=1; 
for i=1:nobs1 
for j=1:nobs2 
    % distances in kilometers. 
    lonH = sin(0.5*(lmat1(i,1)-lmat2(j,1))*pi/180.0)^2; 
    latH = sin(0.5*(lmat1(i,2)-lmat2(j,2))*pi/180.0)^2; 
    hdist(k) = 0.001*6372797.560856*2 ... 
       *asin(sqrt(latH+(cos(lmat1(i,2)*pi/180.0) ... 
       *cos(lmat2(j,2)*pi/180.0))*lonH)); 
    weight(k) = lmat1(i,3)*lmat2(j,3); 
k=k+1; 
end 
end 

toc 

Код вычисляет 10000 X 20000 расстояния и веса.

Elapsed time is 67.124844 seconds. 

Есть ли способ для векторизации обработки с двойным циклом или для выполнения параллельных вычислений? Если в Matlab нет места для улучшения производительности, мне, возможно, придется написать двойные петли в C и вызвать его из Matlab. Я не знаю, как позвонить в C из Matlab, поэтому я задам отдельный вопрос. Благодаря!

Answer 1

bsxfun Использованием, вы можете устранить for петли и необходимость вычисления матриц для каждой комбинации (это должно уменьшить использование памяти). Ниже примерно в шесть раз быстрее, чем ваш исходный код на моем компьютере с помощью R2014b:

nrand1 = 10000; 
nrand2 = 20000; 

% Location matrix 1: [longitude, latitude, w1] 
lmat1=[rand(nrand1,1)-75 rand(nrand1,1)+39 round(rand(nrand1,1)*1000)+1]; 

% Location matrix 2: [longitude, latitude, w2]  
lmat2=[rand(nrand2,1)-75 rand(nrand2,1)+39 round(rand(nrand2,1)*100)+1]; 

p180 = pi/180; 
lonH = sin(0.5*bsxfun(@minus,lmat1(:,1).',lmat2(:,1))*p180).^2; 
latH = sin(0.5*bsxfun(@minus,lmat1(:,2).',lmat2(:,2))*p180).^2; 
hdist = 0.001*6372797.560856*2*asin(sqrt(latH+bsxfun(@times,cos(lmat1(:,2).'*p180),cos(lmat2(:,2)*p180)).*lonH)); 
hdist1 = hdist(:); 
weight1 = bsxfun(@times,lmat1(:,3).',lmat2(:,3)); 
weight1 = weight1(:); 

Обратите внимание, что при использовании переменной p180, математика немного изменилась, так что вы не получите точно такое же значение, но они будет очень близко.

Answer 2

Решение состоит в том, что ваши входы (lmat1 и lmat2) не должны быть такими матрицами, как у вас есть. Каждый из них действительно три вектора. После того, как вы разбили векторы, вы можете создавать массивы, каждая из которых состоит из каждой перестановки lmat1 и lmat2 (что и делает ваш двойной цикл). В этот момент вы можете назвать свою математику как одиночные, полностью векторизованными операции ...

%make your vectors 
lmat1A = rand(nrand1,1)-75; 
lmat1B = rand(nrand1,1)+39; 
lmat1C = round(rand(nrand1,1)*1000)+1 
lmat2A = rand(nrand2,1)-75; 
lmat2B = rand(nrand2,1)+39; 
lmat2C = round(rand(nrand2,1)*1000)+1 

%make every combination 
lmat1A = lmat1A(:)*ones(1,nrand2); 
lmat1B = lmat1B(:)*ones(1,nrand2); 
lmat1C = lmat1C(:)*ones(1,nrand2); 
lmat2A = ones(nrand1,1)*(lmat2A(:)'); 
lmat2B = ones(nrand1,1)*(lmat2B(:)'); 
lmat2C = ones(nrand1,1)*(lmat2C(:)'); 

%do your math 
lonH = sin(0.5*(lmat1A-lmat2A)*pi/180.0).^2; 
latH = sin(0.5*(lmat1B-lmat2B)*pi/180.0).^2; 
hdist = 0.001*6372797.560856*2 ... 
    .*asin(sqrt(latH+(cos(lmat1B*pi/180.0) ... 
    .*cos(lmat2B*pi/180.0)).*lonH)); %use element-wise multiplication 
weight = lmat1C.*lmat2C; 

%reshape your output into vectors (not arrays), which is what your original code does 
lonH = lonH(:) 
latH = latH(:) 
hdist = hdist(:); 
weight = weight(:);