如何在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

Otion 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，它变得最有效。也许是因为内存？有人可以解释一下吗？

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;

Answer 2

您应该尝试使用Tony's trick 在 Matlab 中以最快的方式进行矢量堆叠/平铺。我已经回答了一个类似的问题here。这是Tony's Trick 选项。

% 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)

编辑 1：我进行了一些测试，发现以下内容。在 N=1000 之前，Tony's trick 选项比Option 2 稍快。除此之外，Option 2 再次赶上并变得更快。

可能的原因： 应该是这样的，因为在数组大小可以放入缓存之前，完全矢量化的代码 (Tony's Trick option) 会更快，但是一旦数组大小增长 (N>1000)，它就会溢出到内存缓存中从 CPU 中提取，然后 Matlab 使用一些内部优化将Tony's Trick 代码分解成零碎的代码，使其不再享受完全向量化的好处。