改进MATLAB矩阵构造代码:或者，面向初学者的代码矢量化 [英] Improving MATLAB Matrix Construction Code : Or, code Vectorization for beginners

问题描述

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

I wrote a program in order to construct a portion of a 3-Band Wavelet Transform Matrix. However, given that the size of the matrix is 3^9 X 3^10, it takes a while for MATLAB to finish constructing it. Therefore, I was wondering if there was a way to improve the code I am using to make it run faster. I am using n=10 when running the code.

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;

推荐答案

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

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

How to vectorize without knowing how to vectorize:

First, I'll only discuss vectorizing the first double loop, you can follow the same logic for the second loop.

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

I tried to show a thought process from scratch, so though the final answer is just 2 lines long, it is worthwhile to see how a beginner can try to obtain it.

首先，我建议在简单的情况下按摩"代码，以便对此有所了解.例如，我使用n=3和v=1:6并运行了第一个循环一次，这就是B的样子:

First, I recommend "massaging" the code in simple cases, to get a feel for it. For example, I've used n=3 and v=1:6 and run the first loop once, this is how B looks like:

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

imagesc(B); 

因此，您可以看到我们得到了一个像矩阵这样的阶梯，这很正常！我们只需要为v的正确值分配正确的矩阵索引即可.

So you can see we get a staircase like matrix, which is pretty regular! All we need is to assign the right matrix index to the right value of v and that's it.

有很多方法可以实现这一目标，有些方法比其他方法更优雅.最简单的方法之一是使用功能 find :

There are many ways to achieve that, some more elegant than others. One of the simplest ways is to use the function 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的特定值.例如，索引[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.

The values above are the linear indices of matrix B where values of v were found. Each column represents the linear index of a specific value of v in B. For example, indices [1 29 57 ...] all contain the value v(1), etc... Each row contains v entirely, so, indices [29 38 47 56 65 74] contain v=[v(1) v(2) ... v(6)]. You can notice that for each row the difference between indices is 9, or, each index is separated with a step size N, and there are 6 of them, which is just the length of vector v (also obtained by numel(v)). For the columns, the difference between adjacent elements is 28, or, the step size is M+1.

我们只需要根据此逻辑在适当的索引中分配v的值即可.一种方法是写每个行":

We shall simply need to assign the values of v in the proper indices according to this logic. one way is to write each "row":

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循环中进行操作:

But that is impractical for big N-2, so you can do it in a for loop if you really want:

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

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

Matlab offers a more efficient way to get all the indices at once using bsxfun (this replaces the for loop), for example:

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

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

so now we can use ind to assign v into the matrix N-1 times. For that we need to "flatten" ind into a row vector:

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

并将v连接到自身N-1次(或再增加N-1个v副本):

and concatenate v to itself N-1 times (or make N-1 more copies of v):

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

最终我们得到了答案:

B(ind)=vec;

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

Making a long story short, and writing it all compactly, we get a 2 line solution (given size B is already known: [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重要性会降低)

---

For n=9 the vectorized code is ~850 faster in my machine. (small n will be less significant)

由于获得的矩阵主要由零组成，因此您无需将它们分配给完整的矩阵，而是使用稀疏矩阵，这是完整的该代码(非常相似):

Since the matrix obtained is mostly made up from zeros, you don't need to assign these to a full matrix, instead use a sparse matrix, Here's the full code for that (pretty similar):

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秒钟.

For n=10 I could only run the sparse matrix code (out of memory otherwise), and in my machine it took ~6 seconds.

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

Now try to apply this to the second loop...

这篇关于改进MATLAB矩阵构造代码:或者，面向初学者的代码矢量化的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！