向量化代碼

Question

以下是我的代碼。循環部分非常慢。我不知道是否有一種方法來優化循環部分。

N = 1000000; 
A = rand(N,3); 
B = rand(N,3); 
Dist = sqrt(sum((A - B).^2,2)); 
R = 2; 
id = rangesearch(A,A,0.01); 
result = zeros(N,1); 
for i = 1:N 
    idx = id{i}'; 
    v1 = A(i,:) - A(idx,:); 
    v2 = A(i,:) - B(idx,:); 
    C = cross(v1,v2,2); 
    D = sqrt(sum(C.^2,2))./Dist(idx); 
    result(i) = sum(2 * sqrt(R^2 - D.^2)); 
end 

這裏，A和B是記錄N個點的3D座標的矩陣。首先，我想找到矩陣A中的一個點的鄰居，比如點Ai，並且它的一個鄰居是Aj。我想計算從Ai到Aj-Bj線的距離。這就是我計算交叉產品的原因。最後，我加起來所有的距離。現在，這段代碼在我的電腦上運行了500秒。那麼是否有辦法讓我的代碼更快地運行或以其他方式更快地實現此目標？謝謝。

Answer 1

你的for循環其實很快。

正如@EBH在評論中提到的，您的代碼應該在最新的2016版本中運行，但由於我使用的是早期的2015版本，因此不支持隱式擴展。

原始聲明：rangesearch(A,A,0.01)並不保證您可以爲每個點獲取單個鄰居。事實上，當我與N=10運行，那麼ID總是{1 2 3 4 5 6 7 8 9 10}

固定的循環方法：

tic 
result = zeros(N,1); 
for i = 1:N 
    idx = id{i}'; 
    v1 = bsxfun(@minus, A(i,:), A(idx,:)); 
    v2 = bsxfun(@minus, A(i,:), B(idx,:)); 
    C = cross(v1,v2,2); 
    D = sqrt(sum(C.^2,2))./Dist(idx); 
    result(i) = sum(2 * sqrt(R^2 - D.^2)); 
end 
toc 

Elapsed time is 0.077025 seconds. 

---

方法2：枚舉v1和v2

tic 
idlen = cellfun(@(x) length(x),id); 
idai = cell2mat(arrayfun(@(ii) repmat(ii,1,idlen(ii)), (1:N), 'UniformOutput', false)); 
idx2 = cell2mat(id'); 
V1 = A(idai,:) - A(idx2,:); 
V2 = A(idai,:) - B(idx2,:); 
C2 = cross(V1,V2,2); 
d = @(c,id) sqrt(sum(c.^2,2))./Dist(id); 
r = @(d) sum(2 * sqrt(R^2 - d.^2)); 
result3 = splitapply(@(c,id) r(d(c,id)), C2,idx2', idai'); 
toc 
isequal(result,result3) 


Elapsed time is 0.471092 seconds. 

ans = 

    1 

之間所有可能的組合最慢的線是splitapply

---

方法3：使用cellfun，這並不能保證矢量

tic 
V1 = cellfun(@(Ai,idx) bsxfun(@minus, Ai, A(idx,:)), num2cell(A,2), id, 'UniformOutput', false); 
V2 = cellfun(@(Ai,idx) bsxfun(@minus, Ai, B(idx,:)), num2cell(A,2), id, 'UniformOutput', false); 
C = cellfun(@(v1,v2) cross(v1,v2,2), V1,V2, 'UniformOutput', false); 
D = cellfun(@(c,idx) sqrt(sum(c.^2,2))./Dist(idx), C, id, 'UniformOutput', false); 
result2 = cellfun(@(d) sum(2 * sqrt(R^2 - d.^2)), D, 'UniformOutput', false); 
result2 = cell2mat(result2); 
toc 
isequal(result,result2) 

Elapsed time is 0.122700 seconds. 

ans = 

    1 

---

方法4：使用parfor

tic 
result = zeros(N,1); 
parfor i = 1:N 
    idx = id{i}'; 
    v1 = bsxfun(@minus, A(i,:), A(idx,:)); 
    v2 = bsxfun(@minus, A(i,:), B(idx,:)); 
    C = cross(v1,v2,2); 
    D = sqrt(sum(C.^2,2))./Dist(idx); 
    result(i) = sum(2 * sqrt(R^2 - D.^2)); 
end 
toc 

Elapsed time is 0.177929 seconds. 

parfor抱怨A，B和Dist不能切片，與減慢計算速度。

---

編輯：上述測試使用N=1000。如果我使用N=10000則方法2的時間減少40％，方法3的消耗時間相同，方法4的時間減少約90％。所以，如果你使用的是多核計算機，那麼你可以去parfor。