改进MATLAB矩阵构造代码:或者,为初学者编写代码矢量化
Mat*_*244 6 matlab matrix vectorization wavelet
我写了一个程序,以构建一个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;
bla*_*bla 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
上面的值是矩阵的线性指数,B其中v找到了值.每列代表in 的特定值的线性索引.例如,索引都包含值等等......每行完全包含,因此,索引[29 38 47 56 65 74]包含.您可以注意到,对于每一行,索引之间的差异为9,或者,每个索引用步长分隔,并且其中有6个,这只是向量的长度(也是通过获得).对于列,相邻元素之间的差异为28,或者步长为. vB[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意义不大)
由于获得的矩阵主要由零组成,因此您无需将这些矩阵分配给完整矩阵,而是使用稀疏矩阵,以下是完整代码(非常相似):
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秒钟.
现在尝试将其应用于第二个循环......
虽然你的矩阵具有巨大的尺寸,但它也非常"稀疏",这意味着它的大多数元素都是零.为了提高性能,您可以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)