如何在Matlab中高效地构建一个依赖于索引的矩阵。

Question

在我的Matlab程序中，有几个实例需要创建矩阵，其条目依赖于其索引，并对其执行矩阵向量运算。我想知道如何最有效地实现这一点。

例如，我需要加速以下操作：

N = 1e4;
x = rand(N,1);

% Option 1
tic
I = 1:N;
J = 1:N;
S = zeros(N,N);
for i = 1:N
    for j = 1:N
        S(i,j) = (i+j)/(abs(i-j)+1);
    end
end
a = x'*S*x
fprintf('Option 1 takes %.4f sec\n',toc)
clearvars -except x N

我尝试加快速度，因此我尝试了以下选项：

% Option 2
tic
I = 1:N;
J = 1:N;
Sx = zeros(N,1);
for i = 1:N
    Srow_i = (i+J)./(abs(i-J)+1);
    Sx(i)= Srow_i*x;
end
a = x'*Sx
fprintf('Option 2 takes %.4f sec\n',toc)
clearvars -except x N

并且

% Option 3
tic
I = 1:N;
J = 1:N;
S = bsxfun(@plus,I',J)./(abs(bsxfun(@minus,I',J))+1);
a = x'*S*x
fprintf('Option 3 takes %.4f sec\n',toc)
clearvars -except x N

同时（感谢其中一个答案）

% options 4
tic
[I , J] = meshgrid(1:N,1:N);
S = (I+J) ./ (abs(I-J) + 1);
a = x' * S * x;
fprintf('Option 4 takes %.4f sec\n',toc)
clearvars -except x N

选项2是最有效的。有没有更快的选项来执行这个操作？

更新：

我也尝试了Abhinav的选项。

% Option 5 using Tony's Trick
tic
i = 1:N;
j = (1:N)';
I = i(ones(N,1),:);
J = j(:,ones(N,1));
S = (I+J)./(abs(I-J)+1);
a = x'*S*x;
fprintf('Option 5 takes %.4f sec\n',toc)
clearvars -except x N

看起来最有效的程序取决于N的大小。对于不同的N，我得到以下输出：

N = 100：

Option 1 takes 0.00233 sec
Option 2 takes 0.00276 sec
Option 3 takes 0.00183 sec
Option 4 takes 0.00145 sec
Option 5 takes 0.00185 sec

N = 10000:

Option 1 takes 3.29824 sec
Option 2 takes 0.41597 sec
Option 3 takes 0.72224 sec
Option 4 takes 1.23450 sec
Option 5 takes 1.27717 sec

因此，对于较小的N，选项2是最慢的，但对于更大的N，它变得最有效。可能是因为内存？能有人解释一下吗？

- Optimist

---

你只对 a 的值感兴趣吗？ - obchardon

我包括了计算a以进行公平比较，因为在选项2中它略有不同。 - Optimist

Answer 1

您可以使用meshgrid创建索引，无需循环：

N = 1e4;
[I , J] = meshgrid(1:N,1:N);
x = rand(N,1);
S = (I+J) ./ (abs(I-J) + 1);
a = x' * S * x;

更新：

由于 @Optimist 表明此代码的性能不如 Option2 和 Option3，因此我决定稍微改进 Option2：

N = 1e4;
x = rand(N,1);
Sx = zeros(N,1);
for i = 1:N
    Srow_i = (i+1:i+N)./[i:-1:2,1:N-i+1] ;
    Sx(i)= Srow_i*x;
end
a = x'*Sx;

- rahnema1

---

谢谢你的回答。这确实避免了for循环，但不幸的是比选项2和3都要低效。 - Optimist

@Optimist 通常我认为循环会导致代码效率降低，但是你的例子表明如果正确使用循环可以提高效率。我更新了我的答案并稍微修改了你的option2。 - rahnema1

Answer 2

你应该尝试使用 Tony的技巧 在Matlab中以最快的方式进行向量叠加/平铺。我在这里回答了一个类似的问题。这是Tony的技巧选项。

% Option using Tony's Trick
tic
i = 1:N;
j = (1:N)';
I = i(ones(N,1),:);
J = j(:,ones(N,1));
S = (I+J)./(abs(I-J)+1);

a = x'*S*x
fprintf('Option 1 takes %.4f sec\n',toc)

Edit 1: 我进行了几项测试并发现以下情况。在N=1000之前，Tony's trick选项略快于Option 2。超过此值后，Option 2再次追上并变得更快。

可能的原因: 这是因为，在数组大小适合缓存的情况下，完全向量化的代码（Tony's Trick option）更快，但是一旦数组大小增加（N>1000），它会溢出到远离CPU的内存缓存中，然后Matlab使用某些内部优化将Tony's Trick代码拆分成分段代码，因此不再享受完全向量化的好处。

- Abhinav

---

谢谢您的解释。我已经尝试了您的选项，确实在较小的N（例如N=200）时比选项2更快。在这种情况下，选项4是最快的。对于较大的N，选项2优于其他选项。您的解释听起来很合理。 - Optimist