正確的Matlab parfor切片

Question

我有兩個我想要並行化的嵌套循環。

n=100; 
x=rand(1,n); 
m=5; 
xx=rand(1,m); 

r = zeros(1,m); 
for i=1:n 
    q = ones(1,m); 
    for j=1:n 
     q = q .* (xx-x(j))/(x(i)-x(j)); 
    end 
    r = r + q; 
end 

爲了準備這個功能齶化，我將局部變量更改爲全局變量。

n=100; 
x=rand(1,n); 
m=5; 
xx=rand(1,m); 

r = ones(n,m); 
for i=1:n 
    for j=1:n 
     r(i,:) = r(i,:) .* (xx-x(j))/x(i)-x(j)) 
    end 
end 
r = sum(r,1); 

而不是一次轉化的整體載體，讓我們嘗試它只有一個標量。也使用依賴於i和j的x中最簡單的元素。最後我還刪除了sum。我們可以稍後添加它。

n=100; 
x=rand(1,n); 

r = ones(n,1); 
for i=1:n 
    for j=1:n 
     y = x(i)+x(j); 
     r(i) = r(i) * y; 
    end 
end 

上面的代碼是示例函數，我想並行化。

對於外環i的一次迭代，內循環始終需要訪問相同的向量r(i)。此操作是寫入操作（*=），但命令對此操作無關緊要。

由於嵌套parfor循環不允許在Matlab中，我試圖在一個parfor循環中打包一切。

n=100; 
x=rand(1,n); 

r = ones(n,1); 
parfor k=1:(n*n) 
    %i = floor((k-1)/n)+1; % outer loop 
    %j = mod(k-1,n)+1;  % inner loop 
    [j,i] = ind2sub([n,n],k); 
    y = x(i)+x(j); 
    r(i) = r(i) * y;  % ERROR here 
end 

由於獨立計算，Matlab仍然不知道熱切片它。 因此，我決定將乘法運算移到外面並使用線性索引。

n=100; 
x=rand(1,n); 

r = ones(n,n); 
parfor k=1:(n*n) 
    [j,i] = ind2sub([n,n],k); 
    y = x(i)+x(j); 
    r(k) = y; 
end 
r = prod(r,1); 
r = squeeze(r); % remove singleton dimensions 

雖然這對內部循環中的標量值有效，但它不適用於內部循環中的向量，因爲必須重新計算索引。

n=100; 
x=rand(1,n); 
m=5; 

r = ones(n,n,m); 
parfor k=1:(n*n) 
    [j,i] = ind2sub([n,n],k); 
    y = x(i)+x(j); 
    r((k-1)*m+1:k*m) = y.*(1:m); % ERROR here 
end 
r = prod(r,1); 
r = squeeze(r); % remove singleton dimensions 

儘管它確實有效，但當我重新整形數組時。

n=100; 
x=rand(1,n); 
m=5; 

r = ones(n*n,m); 
parfor k=1:(n*n) 
    [j,i] = ind2sub([n,n],k); 
    y = x(i)+x(j); 
    r(k,:) = y.*(1:m); % ERROR here 
end 
r = reshape(r,n,n,m); 
r = prod(r,2); 
r = squeeze(r); % remove singleton dimensions 

這樣一來，我可以轉換到另一個向量r矢量xx。

n=100; 
x=rand(1,n); 
m=5; 
xx=rand(1,m); 

r = ones(n*n,m); 
parfor k=1:(n*n) 
    [j,i] = ind2sub([n,n],k); 
    y = x(i)+x(j); 
    r(k,:) = y.*xx; % ERROR here 
end 
r = reshape(r,n,n,m); 
r = prod(r,2); 
r = sum(r,1); 
r = reshape(r,size(xx)); % reshape output vector to input vector 

對於我的並行解決方案，我需要一個n*n*m數組，而不是n*m陣列，這似乎非常低效的。 有沒有更好的方式來做我想做的事？ 其他方式的優點是什麼（更漂亮的代碼，更少的CPU，更少的RAM，...）？

UPDATE

在試圖簡化任務，並減少對問題的最低工作示例中的順序，我省略i~=j檢查，使其更容易，雖然導致全面NaN結果。此外，添加此檢查時，代碼的性質會導致所有1結果。爲了使代碼有意義，這些因素僅僅是另一個向量z的權重。

結構複雜的問題如下所示：

n=100; 
x=rand(1,n); 
z=rand(1,n); 
m=5; 
xx=rand(1,m); 

r = zeros(1,m); 
for i=1:n 
    q = ones(1,m); 
    for j=1:n 
     if i~=j 
      q = q .* (xx-x(j))/(x(i)-x(j)); 
     end 
    end 
    r = r + z(i) .* q; 
end 

Answer 1

這個問題不需要任何並行的循環執行。一個問題是x(i)-x(j)被重複計算了很多次。這是低效的。建議的方法精確地計算每個數字一次，並向xx中的每個元素矢量化操作。由於xx是迄今爲止最短的向量，它幾乎完全向量化。如果你想要矢量化最後一個循環，這可能就像隱藏的for循環一樣，它會有更多的內存，代碼會更復雜（如3D矩陣等）。我爲了測試而自由地將分母轉換爲加號。減號會爲所有數字生成NaN。最後一種方法稍微快一點。 n = 10000時約10次。我建議你嘗試一下更精細的基準。

function test() 
% Initiate variables 
n=100; 
x=rand(1,n); 
m=5; 
xx=rand(1,m); 

tic; 
% Alternative 1 
r = zeros(1,m); 
for i=1:n 
    q = ones(1,m); 
    for j=1:n 
     q = q .* (xx-x(j))/(x(i)+x(j)); 
    end 
    r = r + q; 
end 
toc; 

tic; 
% Alternative 2 
xden = bsxfun(@plus, x, x.'); % Calculate denominator 
xnom = repmat(x,n,1); % Calculate nominator 
xfull = (xnom./xden).'; % calculate right term on rhs. 

for (k = 1:m) 
    tmp= prod(xx(k)./xden - xfull); % Split in 2 calculations 
    r2(k) = sum(tmp); % "r = r + xx(k)" 
end 
toc; 

disp(r); 
disp(r2); 

只是在最後的說明。方案2速度更快，但它也是內存昂貴，所以在內存問題的情況下，一個循環更喜歡。此外，並行化時不需要全局變量。如果你需要這個，你可能需要仔細查看你的設計（但是如果代碼很短，沒有什麼關鍵的，那麼你就不需要這麼麻煩）。