3

在我的 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 的大小。对于不同的 NI，得到以下输出：

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

4

2 回答 2

2

您可以使用 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;
于 2016-08-31T09:52:20.987 回答
1

您应该尝试使用Tony 的技巧在 Matlab 中以最快的方式进行矢量堆叠/平铺。我在这里回答了一个类似的问题。这是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代码分解为零碎的代码，使其不再享受完全矢量化的好处。

于 2016-08-31T12:26:08.027 回答