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)。此访问是一个 write 操作 (*=)，但 order 对于此操作无关紧要。

因为嵌套 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

这样，我可以将向量 xx 转换为另一个向量 r。

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，...）？

更新

为了尝试简化任务并将其减少到问题的最小工作示例，我省略了 i~=j 的检查以使其更容易，尽管导致所有 NaN 结果。此外，代码的性质在添加此检查时会导致 all 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

这道题不需要任何并行for循环来执行。一个问题是 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 更快，但它的内存也很昂贵，因此在内存问题的情况下最好使用循环。此外，在并行化的情况下不需要全局变量。如果你需要这个，你可能需要检查你的设计（但如果代码很短，没有一些关键的东西，那么你就不需要那么麻烦了）。

Matlab parfor 切片正确

Matlab parfor slice correctly

arrays

optimization

matlab

slice

parfor