matlab - 改进 MATLAB 矩阵构造代码 : Or, 代码 Vectorization for beginners

Question

我编写了一个程序来构建 3 波段小波变换矩阵的一部分。但是，鉴于矩阵的大小为 3^9 X 3^10，MATLAB 需要一段时间才能完成构造。因此，我想知道是否有一种方法可以改进我正在使用的代码以使其运行得更快。我在运行代码时使用 n=10。

B=zeros(3^(n-1),3^n);
v=[-0.117377016134830 0.54433105395181 -0.0187057473531300 -0.699119564792890 -0.136082763487960 0.426954037816980 ];

for j=1:3^(n-1)-1 
    for k=1:3^n;
        if k>6+3*(j-1) || k<=3*(j-1)
            B(j,k)=0;
        else 
            B(j,k)=v(k-3*(j-1));
        end                
    end
end
j=3^(n-1);
    for k=1:3^n
        if k<=3
            B(j,k)=v(k+3);
        elseif k<=3^n-3
            B(j,k)=0;
        else 
            B(j,k)=v(k-3*(j-1));
        end
    end

W=B;

Answer 1

不知道如何向量化的情况下如何向量化:

首先，我将只讨论矢量化第一个双循环，您可以按照相同的逻辑处理第二个循环。

我试图展示一个从头开始的思考过程，所以虽然最终答案只有两行，但值得看看初学者如何尝试获得它。

首先，我建议在简单情况下“按摩”代码，以感受一下。例如，我使用了 n=3 和 v=1:6 并运行了第一个循环一次，这就是 B 的样子:

[N M]=size(B)
N =
     9
M =
    27

imagesc(B); 

所以你可以看到我们得到了一个像矩阵一样的楼梯，非常规则!我们只需要将正确的矩阵索引分配给 v 的正确值即可。

有很多方法可以实现这一点，有些方法比其他方法更优雅。最简单的方法之一是使用函数 find :

pos=[find(B==v(1)),find(B==v(2)),find(B==v(3)),...
     find(B==v(4)),find(B==v(5)),find(B==v(6))]

pos =
     1    10    19    28    37    46
    29    38    47    56    65    74
    57    66    75    84    93   102
    85    94   103   112   121   130
   113   122   131   140   149   158
   141   150   159   168   177   186
   169   178   187   196   205   214
   197   206   215   224   233   242

上面的值是 linear indices矩阵 B 的 ，其中找到了 v 的值。每列代表 linear index B 中 v 的特定值。例如，索引 [1 29 57 ...] 都包含值 v(1)，等等......每行包含 v因此，索引 [29 38 47 56 65 74] 完全包含 v=[v(1) v(2) ... v(6)]。你可以注意到对于每一行，索引之间的差异是 9，或者，每个索引以步长 N 分隔，并且有 6 个，这恰好是向量 的长度v(也是通过numel(v)获得的)。对于列，相邻元素之间的差为 28，或者，步长为 M+1。

我们只需要根据这个逻辑在适当的索引中分配 v 的值。一种方法是写每个“行”:

B([1:N:numel(v)*N]+(M+1)*0)=v;
B([1:N:numel(v)*N]+(M+1)*1)=v;
...
B([1:N:numel(v)*N]+(M+1)*(N-2))=v;

但这对于大的 N-2 来说是不切实际的，所以如果你真的想要的话，你可以在 for 循环中完成:

for kk=0:N-2;
     B([1:N:numel(v)*N]+(M+1)*kk)=v;
end

Matlab 提供了一种使用 bsxfun(这取代了 for 循环)一次获取所有索引的更有效方法，例如:

ind=bsxfun(@plus,1:N:N*numel(v),[0:(M+1):M*(N-1)+1]')

所以现在我们可以使用ind 将v 分配给矩阵N-1 次。为此，我们需要将 ind“展平”为一个行向量:

ind=reshape(ind.',1,[]);

并将 v 与其本身连接 N-1 次(或使 v 的 N-1 个副本):

vec=repmat(v,[1 N-1]);

终于得到答案了:

B(ind)=vec;

长话短说，紧凑地编写所有内容，我们得到一个 2 行解决方案(给定大小 B 已知:[N M]=size(B)) :

---

ind=bsxfun(@plus,1:N:N*numel(v),[0:(M+1):M*(N-1)+1]');
B(reshape(ind.',1,[]))=repmat(v,[1 N-1]);

---

对于 n=9，矢量化代码在我的机器上快了 ~850。 (小的 n 将不那么重要)

由于得到的矩阵主要由零组成，所以不需要将它们赋给一个满矩阵，而是使用稀疏矩阵，这里是完整的代码(非常相似):

N=3^(n-1);
M=3^n;
S=sparse([],[],[],N,M);
ind=bsxfun(@plus,1:N:N*numel(v),[0:(M+1):M*(N-1)+1]');
S(reshape(ind.',1,[]))=repmat(v,[1 N-1]);

对于 n=10，我只能运行稀疏矩阵代码(否则内存不足)，而在我的机器上大约需要 6 秒。

现在尝试将其应用于第二个循环...