📄 matlab vector guid代码矢量化指南.txt
字号:
矢量化技术,并示范如何编写你自己的高效代码。此外,前者还有一个
作用:Unique函数提供了一些超出我们要求的额外功能,这可能降低代
码的执行速度。
假设我们不只是要返回集合x,而且要知道在原始的矩阵里每个相异元素
出现了多少个“复本”。一旦我们对x排序并进行了差分,我们可以用
find来确定差分变化的位置。再将这个变化位置进行差分,就可以得到
复本的数目。这就是"diff of find of diff"的技巧。基于以上的讨论,
我们有:
% Find the redundancy in a vector x
x = sort(x(:));
difference = diff([x;max(x)+1]);
count = diff(find([1;difference]));
y = x(find(difference));
plot(y,count)
这个图画出了x中每个相异元素出现的复本数。注意,在这里我们避开了
NaN,因为find不返回NaN元素的索引值。但是,作为特例,NaN和Inf
的复本数可以容易地计算出来:
count_nans = sum(isnan(x(:)));
count_infs = sum(isinf(x(:)));
另一个用于求和或者计数运算的矢量化技巧是用类似建立稀疏矩阵的方
法实现的。这还将在第9节中作更加详细的讨论.
-------------------------------------------------------
8)稀疏矩阵结构(Sparse Matrix Structures)
在某些情况下,你可以使用稀疏矩阵来增加计算的效率。如果你构造一
个大的中间矩阵,通常矢量化更加容易。在某些情况下,你可以充分利
用稀疏矩阵结构来矢量化代码,而对于这个中间矩阵不需要大的存储空
间。
假设在上一个例子中,你事先知道集合y的域是整数的子集,
{k+1,k+2,...k+n};即,
y = (1:n) + k
例如,这样的数据可能代表一个调色板的索引值。然后,你就可以对集
合中每个元素的出现进行计数(构建色彩直方图?译者)。这是对上一
节中"diff of find of diff"技巧的一种变形。
现在让我们来构造一个大的m x n矩阵A,这里m是原始x向量中的元素数,
n是集合y中的元素数。
A(i,j) = 1 if x(i) = y(j)
0 otherwise
回想一下第3节和第4节,你可能认为我们需要从x和y来构造矩阵A。如果
当然可以,但要消耗许多存储空间。我们可以做得更好,因为我们知道,
矩阵A中的多数元素为0,x中的每个元素对应的行上只有一个值为1。
以下就是构造矩阵的方法(注意到y(j) = k+j,根据以上的公式):
x = sort(x(:));
A = sparse(1:length(x), x+k, 1, length(x), n);
现在我们对A的列进行求和,得到出现次数。
count = sum(A);
在这种情况下,我们不必明确地形成排序向量y,因为我们事先知道
y = 1:n + k.
这里的关键是使用数据,(也就是说,用x控制矩阵A的结构)。由于x在
一个已知范围内取整数值,我们可以更加有效地构造矩阵。
假设你要给一个很大矩阵的每一列乘以相同的向量。使用稀疏矩阵,不仅
可以节省空间,并且要比在第5节介绍的方法更加快速. 下面是它的工作
方式:
F = rand(1024,1024);
x = rand(1024,1);
% 对F的所有行进行点型乘法.
Y = F * diag(sparse(x));
% 对F的所有列进行点型乘法.
Y = diag(sparse(x)) * F;
我们充分利用矩阵乘法算符来执行大规模运算,并使用稀疏矩阵以防止临
时变量变得太大。
--------------------------------------------------------
9)附加的例子(Additional Examples)
下面的例子使用一些在本技术手册中讨论过的方法,以及其它一些相关方
法。请尝试使用tic 和toc(或t=cputime和cputime-t),看一下速度加快
的效果。
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
用于计算数组的双重for循环。
使用的工具:数组乘法
优化前:
A = magic(100);
B = pascal(100);
for j = 1:100
for k = 1:100;
X(j,k) = sqrt(A(j,k)) * (B(j,k) - 1);
end
end
优化后:
A = magic(100);
B = pascal(100);
X = sqrt(A).*(B-1);
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
用一个循环建立一个向量,其元素依赖于前一个元素
使用的工具:FILTER, CUMSUM, CUMPROD
优化前:
A = 1;
L = 1000;
for i = 1:L
A(i+1) = 2*A(i)+1;
end
优化后:
L = 1000;
A = filter([1],[1 -2],ones(1,L+1));
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
如果你的向量构造只使用加法或乘法,你可使用cumsum或cumprod函数。
优化前:
n=10000;
V_B=100*ones(1,n);
V_B2=100*ones(1,n);
ScaleFactor=rand(1,n-1);
for i = 2:n
V_B(i) = V_B(i-1)*(1+ScaleFactor(i-1));
end
for i=2:n
V_B2(i) = V_B2(i-1)+3;
end
优化后:
n=10000;
V_A=100*ones(1,n);
V_A2 = 100*ones(1,n);
ScaleFactor=rand(1,n-1);
V_A=cumprod([100 1+ScaleFactor]);
V_A2=cumsum([100 3*ones(1,n-1)]);
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
向量累加,每5个元素进行一次:
工具:CUMSUM , 向量索引
优化前:
% Use an arbitrary vector, x
x = 1:10000;
y = [];
for n = 5:5:length(x)
y = [y sum(x(1:n))];
end
优化后(使用预分配):
x = 1:10000;
ylength = (length(x) - mod(length(x),5))/5;
% Avoid using ZEROS command during preallocation
y(1:ylength) = 0;
for n = 5:5:length(x)
y(n/5) = sum(x(1:n));
end
优化后(使用矢量化,不再需要预分配):
x = 1:10000;
cums = cumsum(x);
y = cums(5:5:length(x));
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
操作一个向量,当某个元素的后继元素为0时,重复这个元素:
工具:FIND, CUMSUM, DIFF
任务:我们要操作一个由非零数值和零组成的向量,要求把零替换成为
它前面的非零数值。例如,我们要转换下面的向量:
a=2; b=3; c=5; d=15; e=11;
x = [a 0 0 0 b 0 0 c 0 0 0 0 d 0 e 0 0 0 0 0];
为:
x = [a a a a b b b c c c c c d d e e e e e e];
解(diff和cumsum是反函数):
valind = find(x);
x(valind(2:end)) = diff(x(valind));
x = cumsum(x);
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
将向量的元素累加到特定位置上
工具:SPARSE
优化前:
% The values we are summing at designated indices
values = [20 15 45 50 75 10 15 15 35 40 10];
% The indices associated with the values are summed cumulatively
indices = [2 4 4 1 3 4 2 1 3 3 1];
totals = zeros(max(indices),1);
for n = 1:length(indices)
totals(indices(n)) = totals(indices(n)) + values(n);
end
优化后:
indices = [2 4 4 1 3 4 2 1 3 3 1];
totals = full(sparse(indices,1,values));
注意:这一方法开辟了稀疏矩阵的新用途。在使用sparse命令创建稀疏矩阵
时,它是对分配到同一个索引的所有值求和,而不是替代已有的数值。这称
为"向量累加",是MATLAB处理稀疏矩阵的方式。
========================================================================
=
三、更多资源
--------------------------------------------------------------
10)矩阵索引和运算
下面的MATLAB文摘讨论矩阵索引。它比本技术手册的第1节提供了更加
详细的信息
MATLAB Digest: Matrix Indexing in MATLAB
http://www.mathworks.com/company/digest/sept01/matrix.shtml
下面的说明链接将指导你如何使用MATLAB中的矩阵操作。含括矩阵创建
、索引、操作、数组运算、矩阵运算以及其它主题。
MATLAB Documentation: Getting Started
http://www.mathworks.com/access/helpdesk/help/techdoc/
learn_MATLAB/learn_MATLAB.shtml
Peter Acklam是我们的一个power user,他建立了一个网页提供MATLAB
数组处理中的一些技巧。
Peter Acklam's MATLAB array manipulation tips and tricks
http://home.online.no/~pjacklam/MATLAB/doc/mtt/index.html
---------------------------------------------------------------
11)矩阵存储管理(Matrix Memory Management)
有关预分配技术如何加快计算速度的更多的信息,参见下面的联机解决方案:
26623: How do I pre-allocate memory when using MATLAB?
http://www.mathworks.com/support/solutions/data/26623.shtml
下面的技术手册是关于MATLAB存储管理方面的一个包罗广泛的指南:
The Technical Support Guide to Memory Management
http://www.mathworks.com/support/tech-notes/1100/1106.shtml
----------------------------------------------------------------
12)发挥MATLAB的最高性能(Maximizing MATLAB Performance)
这个技术手册对代码矢量化进行一般的讨论,而许多时候,我们的目的是加快
代码的速度。因此,本节提供一些附加的资源,它们涉及如何是代码达到最高
性能。
加速代码的常用技巧:
3257: How do I increase the speed or performance of MATLAB?
http://www.mathworks.com/support/solutions/data/3257.shtml
编译你的代码并从MATLAB中调用它
22512: What are the advantages of using the MATLAB Compiler to
compile my M-files?
http://www.mathworks.com/support/solutions/data/22512.shtml
======================================================================
结论
正如本文所介绍的,针对代码矢量化和代码提速,有一些常规方法,同时,还有
一些非常规方法,正如第9节最后那个例子(对sparse函数进行独特利用)那样。
以上的讨论远没有穷尽MATLAB中所有矢量化技术。它仅仅是代码矢量化的思想和
理论的一个导引。
--
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -