6

我编写了一个程序来构造一个 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;
4

3 回答 3

11

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

首先，我将只讨论对第一个双循环进行矢量化，您可以对第二个循环遵循相同的逻辑。

我试图从头开始展示一个思考过程，所以虽然最终答案只有 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

上面的值是找到值的矩阵的线性索引。每列表示in的特定值的线性索引。例如，索引都包含值等...每一行都包含完整，因此，索引 [29 38 47 56 65 74] 包含。您可以注意到，对于每一行，索引之间的差异为 9，或者，每个索引以步长 分隔，其中有 6 个，这只是向量的长度（也是由 获得的）。对于列，相邻元素之间的差异为 28，或者步长为。 BvvB[1 29 57 ...]v(1)vv=[v(1) v(2) ... v(6)]Nvnumel(v)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;

但这对于 big 是不切实际的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不那么重要）

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

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

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

于 2013-07-11T05:32:27.050 回答
8

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

通过构造稀疏矩阵的坐标形式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我们得到以下结果（为了比较，我还包括了bsxfun建议的基于方法的结果natan）：

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)
于 2013-07-12T05:03:33.540 回答
2

您可以使用一种巧妙的方法来创建一个块对角矩阵，如下所示：

>> v=[-0.117377016134830 0.54433105395181 -0.0187057473531300 ...
          -0.699119564792890 -0.136082763487960 0.426954037816980]；
>>lendiff=长度（v）-3；
>> B=repmat([v zeros(1,3^n-lendiff)],3^(n-1),1);
>> B=重塑(B',3^(n),3^(n-1)+1);
>> B(:,end-1)=B(:,end-1)+B(:,end);
>> B=B(:,1:end-1)';

在这里，lendiff用于创建行的 3^{n-1} 个副本，v后跟零，长度为 3^n+3，因此矩阵大小为 [3^{n-1} 3^n+3] .

该矩阵被重塑为大小 [3^n 3^{n-1}+1] 以创建移位。额外的列需要添加到最后，B 需要转置。

不过应该快得多。

编辑

看到 Darren 的解决方案并意识到它也reshape适用于稀疏矩阵，让我想出了这个 - 没有 for 循环（未编码原始解决方案）。

首先是开始的值：

>> v=[-0.117377016134830  ...
       0.54433105395181   ...
      -0.0187057473531300 ...
      -0.699119564792890  ...
      -0.136082763487960  ...
       0.426954037816980];    
>> rows = 3^(n-1);                  % same number of rows
>> cols = 3^(n)+3;                  % add 3 cols to implement the shifts

然后使矩阵每行增加 3 列

>> row=(1:rows)'*ones(1,length(v)); % row number where each copy of v is stored'
>> col=ones(rows,1)*(1:length(v));  % place v at the start columns of each row
>> val=ones(rows,1)*v;              % fill in the values of v at those positions
>> B=sparse(row,col,val,rows,cols); % make the matrix B[rows cols+3], but now sparse

然后重塑以实现移位（额外的行，正确的列数）

>> B=reshape(B',3^(n),rows+1);      % reshape into B[3^n rows+1], shifted v per row'
>> B(1:3,end-1)=B(1:3,end);         % the extra column contains last 3 values of v
>> B=B(:,1:end-1)';                 % delete extra column after copying, transpose

对于 n=4,5,6,7，这会导致s中的 cpu 时间：

n    original    new version
4    0.033       0.000
5    0.206       0.000
6    1.906       0.000
7    16.311      0.000

由 Profiler 测量。对于原始版本，我无法运行 n>7，但新版本给出了

n    new version
8    0.002
9    0.009
10   0.022
11   0.062
12   0.187
13   0.540
14   1.529
15   4.210

这就是我的 RAM 能走多远 :)

于 2013-07-17T16:14:05.440 回答