正确的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速度更快，但它也是内存昂贵，所以在内存问题的情况下，一个循环更喜欢。此外，并行化时不需要全局变量。如果你需要这个，你可能需要仔细查看你的设计（但是如果代码很短，没有什么关键的，那么你就不需要这么麻烦）。