改进MATLAB矩阵构造代码:面向初学者的代码向量化

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

如何在不知道如何矢量化的情况下进行矢量化：

首先，我将只讨论向量化第一个双循环，你可以遵循第二个循环的相同逻辑。

我试图从头开始展示一个思考过程，所以尽管最终答案只有2行长，但初学者如何才能获得它是值得一看的。

首先，我建议在简单的情况下“按摩”代码，以获得它的感觉。例如，我使用了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

上面的值是矩阵B的，其中找到了v的值。每一列表示B中v的特定值的linear index。例如，索引[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次(或者再制作N-1个v副本)：

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秒。

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

Answer 2

虽然你的矩阵具有巨大的维度，但它也是非常“稀疏”的，这意味着它的大部分元素都是零。为了提高性能，您可以通过确保只对矩阵的非零部分进行操作来利用MATLAB的稀疏矩阵支持。

通过构造稀疏矩阵的coordinate form，可以有效地构造MATLAB中的稀疏矩阵。这意味着必须定义三个数组，分别为矩阵中的每个非零项定义行、列和值。这意味着我们不会通过传统的A(i,j) = x语法赋值，而是将该非零条目附加到我们的稀疏索引结构中：

row(pos+1) = i;
col(pos+1) = j;
val(pos+1) = x;
% pos is the current position within the sparse indexing arrays!

一旦我们在稀疏索引数组中有了完整的非零集，我们就可以使用sparse命令来构建矩阵了。

对于这个问题，我们为每行添加最多六个非零条目，允许我们提前分配稀疏索引数组。变量pos跟踪我们在索引数组中的当前位置。

rows = 3^(n-1);
cols = 3^(n+0);

% setup the sparse indexing arrays for non-
% zero elements of matrix B
row = zeros(rows*6,1);
col = zeros(rows*6,1);
val = zeros(rows*6,1);
pos = +0;

现在，我们可以通过向稀疏索引数组添加任何非零条目来构建矩阵。因为我们只关心非零条目，所以我们也只循环矩阵的非零部分。

我把最后一行的逻辑留给您来填写！

for j = 1 : 3^(n-1)
    if (j < 3^(n-1))

% add entries for a general row
    for k = max(1,3*(j-1)+1) : min(3^n,3*(j-1)+6)             
        pos = pos+1;
        row(pos) = j;
        col(pos) = k;
        val(pos) = v(k-3*(j-1));                
    end

    else

% add entries for final row - todo!!

    end
end

因为我们没有为每一行添加六个非零，所以我们可能过度分配了稀疏索引数组，所以我们将它们削减到实际使用的大小。

% only keep the sparse indexing that we've used
row = row(1:pos);
col = col(1:pos);
val = val(1:pos);

现在可以使用sparse命令构建最终的矩阵。

% build the actual sparse matrix
B = sparse(row,col,val,rows,cols);

可以通过整理上面的代码片段来运行代码。对于n = 9，我们得到了以下结果(为了进行比较，我还包括了natan建议的基于bsxfun的方法的结果)：

Elapsed time is 2.770617 seconds. (original)
Elapsed time is 0.006305 seconds. (sparse indexing)
Elapsed time is 0.261078 seconds. (bsxfun)

原始代码耗尽了n = 10的内存，尽管这两种稀疏方法仍然可以使用：

Elapsed time is 0.019846 seconds. (sparse indexing)
Elapsed time is 2.133946 seconds. (bsxfun)