MatLab 向量化技巧.docx

MatLab 向量化技巧.docx

《MatLab 向量化技巧.docx》由会员分享，可在线阅读，更多相关《MatLab 向量化技巧.docx（11页珍藏版）》请在冰豆网上搜索。

MatLab向量化技巧

MatLab向量化技巧

家家@2007-05-1713:

53

大家都知道MatLab是一种解释性语言，它的长处在矩阵运算。

因此需要将问题尽量用矩阵表示，并且要避免对单个矩阵元素操作，

而用整体矩阵运算。

一个重要过程就是所谓将问题矢量化（vectorizing）。

一个问题是如果需要对矩阵的不同元素作不同操作，那么怎么办？

举个例子

对矩阵X，作A=5*sin（X）/X，显然当X某元素->0时，A对应元素->5。

按照

通常在其他语言思路编程为：

Nx=2000;

X=2.0*pi*（rand（Nx,Nx）-0.5）;

X（1:

Nx/2,:

）=0.0;        %生成Nx X Nx矩阵X，一半为随机分布，一半为0。

A=zeros（Nx,Nx）;

for m=1:

Nx

for n=1:

Nx

    if（abs（X（m,n））>=eps）

        A（m,n）=5.0.*sin（X（m,n））./X（m,n）;

    else

        A（m,n）=5.0;

    end

end

end                                                     %（ProgramNo.1）

以上代码通过循环语句对每一个矩阵元素进行单独操作，这段代码在MatLab里效

率特别低。

为避免循环，考虑如下方法：

方法I：

用所谓logical subscripting。

Nx=2000;

X=2.0*pi*（rand（Nx,Nx）-0.5）;

X（1:

Nx/2,:

）=0.0;

A=zeros（Nx,Nx）;

I=（abs（X）>=eps）;

A=I.*5.0.*sin（X）./（X+（1-I））+（1-I）.*5.0;                 %（ProgramNo.2）

或者上面最后一行用如下代码代替：

A（I）=5.0.*sin（X（I））./X（I）;

A（~I）=5.0;                                              %（ProgramNo.3）

一个看起来更有效的办法是先生成一个所有元素都是5的A矩阵。

代码如下：

Nx=2000;

X=2.0*pi*（rand（Nx,Nx）-0.5）;

X（1:

Nx/2,:

）=0.0;

A=5.0;

A=A（ones（Nx,Nx））;  % 用Tony's trick生成矩阵A，所有元素=5

I=（abs（X）>=eps）;

A（I）=5.0.*sin（X（I））./X（I）;                             %（ProgramNo.4）

方法II：

用find函数。

Nx=2000;

X=2.0*pi*（rand（Nx,Nx）-0.5）;

X（1:

Nx/2,:

）=0.0;

A=5.0;

A=A（ones（Nx,Nx））;

I=find（abs（X）>=eps）;

A（I）=5.0.*sin（X（I））./X（I）;                             %（ProgramNo.5）

我在P4 1.6G，512M DDR266 上作了个测试，每个程序执行5遍，选最小执行时间：

ProgramNo.     time（s）

1              71.20

2              3.31

3              2.64

4              2.80

5              3.22

显然ProgramNo.3是最快的。

如果X矩阵中0元素很少，那么ProgramNo.2、

ProgramNo.3以及ProgramNo.4的执行时间会非常接近。

大家可以用profiler分析每一个语句所用时间情况。

以上是我用过的两种方法，有没有更有效的方法呢？

还有个问题就是，如何将主对角元素的Index加入到变量I中？

例如需要将满足abs（X）

phi=0:

pi/20:

pi/2;

phi=phi';

theta=0:

pi/30:

pi;

A=repmat（phi,1,31）;

B=repmat（theta,11,1）;

x1=10*sin（A）.*cos（300*pi*20*B/1000-pi）-20*cos（B）;

浅析MATLAB程序设计语言的效率（2007-07-0617:

18:

20）

标签：

学术科学统计学习计算机视觉matlab编程语言程序设计分类：

学术科学

在统计学习和计算机视觉领域，MATLAB是进行实验的最主要的工具之一。

如果你觉得MATLAB很慢，可是这往往不是MATLAB的问题，而是因为自己程序没有写好。

MATLAB是我非常欣赏的一种语言，堪称灵活易用和高效运算的结合典范。

作为一种解释语言，我认为MATLAB进行高效运算的秘诀有几个方面：

（1）JIT即时编译技术。

Matlab在第一次加载运行一个函数的时候，会在后台对它进行某种程度的编译和优化，使得后续运行更快。

因此，除了第一次运行需要进行语句解释之外，后面运行的其实是已经放在内存中的经过优化的中间码。

所以，很多时候我们可能会看到第一次运行一个新函数，比它后面几次的运行明显慢一些。

不过目前的JIT技术还不是非常成熟，和标准的编译语言相比还有相当差距，仅凭这个MATLAB还不能称为高效。

（2）向量化（Vectorization）。

这是MATLAB最著名的特色，尤其对于计算机视觉和图像分析等多个研究领域，利用MATLAB的向量化程序设计有其天然的优势。

向量化配合经过高度优化的数值运算引擎，是其高效运算的最重要的基石——很多MATLAB的使用者都明白这一点。

能够转化成矩阵操作的规则运算过程，使用MATLAB计算远比自己手工用C/C++实现高效。

相比之下，MATLAB在关键的核心运算上的实现往往可以比自己用C/C++按照标准的routine进行实现快几十倍。

不过，这也不能完全归功与MATLAB程序设计平台，其实MATLAB是建基于BLAS和LAPACK等的基础数学运算包的基础上的。

Intel和AMD都发布了这些东西的vendor-implementation，并且针对各自的CPU进行了大量的优化。

因此，在一般情况下，我们都不会选择使用C/C++去实现数值运算的关键代码（学习数值分析课者除外）。

但是，哲学的思想告诉我们世界上没有十全十美的东西，任何事情都是有利有弊的。

因此，也不能把MATLAB向量化的优势无限制的扩大，这种优势在某些case下往往变成shortcoming。

下面举些例子说明有时for-loop比vectorization更适合：

*（a）粗粒度的算法流程控制。

比如，你要实现一个迭代算法，循环做一个事情直到收敛。

只要循环几次或者几十次，但是每个循环内部要进行大量的复杂运算，那么你就没有必要费心把这层循环给vectorize掉了。

除非收敛结果有某种close-form的解析解。

*（b）如果向量化可能导致产生巨型矩阵，则使用前要三思！

这是很头疼的问题，相信大家与我有同感。

很多情况下，向量化是一种用空间换时间的行为，就是通过把数据组织成某种方式，从而使得内建的高效引擎能得以应用。

但是，有些时候要处理大量的数据（尤其是图像数据），可能导致其组织过程需要耗费额外的数百兆乃至千兆内存空间，那么这有可能造成效率的不升反降。

原因有几个方面：

第一，数据组织过程也是需要时间的，最起码它也需要大量的操作进行密集的内存读写和用于定位的偏移量计算。

这方面增加的时间vs.向量化节省的运算时间，孰轻孰重，需要衡量。

第二，分配额外的大块内存是一件非常耗时的事情，它可能导致虚拟内存的使用，那么当对这块矩阵进行读写和计算时可能涉及频繁的内存与外存交换区的I/O，势必造成效率的严重下降。

这时应该做的是对程序进行重新思考和设计，降低其对内存的耗用。

第三，vectorization有些时候还可能增加实际运算次数，这往往出现在不适合的向量化过程中。

这样，即使你通过生搬硬套的向量化过程让操作变成矩阵运算，但是增加的无用计算使得即使是更高效的引擎也无法挽回你的损失了。

说了这么多，并不是要把Vectorization贬的一钱不值，而是希望提醒大家在做一些事情的时候，要考虑得周全一些。

（3）InplaceOperation。

这是一个很有趣的事情，MATLAB的对象管理策略是Copy-on-Write，就是说你把一个矩阵传给一个函数的时候，会先传引用，而不产生副本，而当函数要对这个矩阵修改的时候，才会制造出它的副本出来，让函数去修改。

这样听上去很完美，既保护了输入矩阵不被改变，又避免了大量无谓的复制。

不过，这个策略真的没有缺陷么？

可以看看下面的例子

A=rand（2000,2000）;

fori=1:

1000

A=f（A）;

end

functionA=f（A）

A

（1）=A

（1）+1;%atemporarycopyofAisproducedformodification

return;%themodifiedtemporarycopyisassignedtoA

在上面的代码中，只是想对A的第一个元素做1000次加法，结果却导致了对整个有四百万元素的大矩阵做了1000次副本复制！

而且这些副本其实没有用，只要f（A）直接对A的原矩阵进行修改，这些巨大的浪费就能避免。

你可能跟我argue说，这里完全可以写成A

（1）=A

（1）+1000，不用这么搞。

我想说明的是，我只是想举一个简单例子说明，Copy-on-write的策略为什么可能造成巨大的效率浪费。

其实，很多小修改没法像上面那么容易合并。

以前，为了解决这个问题，f函数的作者需要自己写一个mexfunction来避开MATLAB的保护机制，直接改写A。

后来matlab自己也意识到这种策略在实际中的问题，于是他们改写了解释器，采用了一种办法：

A=f（A,...）当它发现这样的定义和调用的时候，它会很智能地了解到你其实是想改写A这个输入，于是它就把操作改成inplace，直接对A进行，避免拷贝。

这种修改的策略，在2006b后才比较正式的运用，在旧版Matlab中仍旧不同程度地存在上述问题。

我在2007a中，测试了两种f的写法：

functiony=f0（x）y=x+1;end

functionx=f1（x）x=x+1;end

确实发现大量调用f1比f0快了好多倍。

所以，如果确实是想改变input的话，函数定义中要让inputparameter和outputparameter同名，以触发这种智能的解释机制。

并且建议升级到最新的2007a版本，这个版本在这方面进行了重要改善。

最近发现的一个技巧，台湾去年就有人用了，在simwe、研学、振动上搜了搜都没有这个函数的介绍。

因此特转过来。

我们知道，MATLAB编程核心思想之一就是向量化。

MATLAB的很多built-in函数，向sin,cos,find等等都支持向量运算。

但是，很多时候我们编写的函数的输入变量是标量，而我们又要对很多组参数进行函数调用，MATLAB7以前的版本中我们只能通过循环来实现，如果参数的维数增加，就会出现循环套循环的现象，效率恐怖。

从MATLAB7.1开始，MATLAB新增arrayfun这个built-in函数来实现将任意函数应用到数组内包括结构在内的所有元素。

这样很多以前不可避免的循环现在可以向量化了。

举例如下：

例1：

生成一个这样的n*n矩阵a：

a（i,j）=dblquad（@（u,v）sin（u）*sqrt（v）,0,i,0,j）。

以n=10为例，

以前我们可能这样做：

代码:

a=zeros（10）;

forii=1:

10

forjj=1:

10

a（ii,jj）=dblquad（@（u,v）sin（u）*sqrt（v）,0,ii,0,jj）;

end

end

现在我们只需这样：

代码:

[J,I]=meshgrid（1:

10）;

a1=arrayfun（@（ii,jj）dblquad（@（u,v）sin（u）*sqrt（v）,0,ii,0,jj）,I,J）;

例2：

验证角谷猜想，一个正整数n,如果是偶数除以2，如果是奇数乘以3加1，得到的新数继续按上述规则运算，最后结果都为1。

验证1到100000内的正整数。

先编写单个数的验证函数

代码:

functionf=jiaogu（n）

while（n>1）

if（mod（n,2）==1）

n=n*3+1;

elseif（mod（n,2）==0）

n=n/2;

elsebreak;

end

end

f=n;

验证1:

100000,以前我们可能这样做：

代码:

a=zeros（1,100000）;

fork=1:

100000

a（k）=jiaogu（k）;

end

all（a）

现在我们只需这样：

代码:

all（arrayfun（@jiaogu,1:

100000））

举这两个例子目的是为了让大家知道很多以前不能优化的程序利用这个函数还可以进一步优化。

具体了解请docarrayfun。

另外cellfun和structfun也提供了类似的功能，关于其用法，帮助文档里写的很详细。

向量化操作的又一重要函数accumarray的用法总结

向量化编程是MATLAB编程区别于其他语言的最重要特征之一，MATLAB不断增强其向量化编程的能力。

版本7以后出现的accumarray函数就是一个很好的例子。

accumarray函数最早出现于7.0版（R14），在随后的7.1（R14SP3），7.2（R2006a）版里又对其功能进行增强。

大家可以在帮助文档里了解其基本信息，这里结合自己的使用经验给出其用法的举例，大家可以看看其效果。

例1：

生成一个100000*1的向量，其元素服从[1，100000]之间离散均匀分布，找出都有哪些元素出现了

a=unidrnd（100000,100000,1）;

方法1：

tic;b=union（a,a）;toc;

Elapsedtimeis0.077001seconds.

方法2：

tic;c=unique（a）;toc;

Elapsedtimeis0.035060seconds.

方法3

tic,d=accumarray（a,1,[100000,1]）;e=find（d）;toc

Elapsedtimeis0.016507seconds.

d给出了每个元素出现的个数，如果还需要这方面的信息，无疑方法3占绝对优势。

例2：

在1000*1000的正方形区域内随机生成100000个点（坐标值是整数），统计每个坐标点上生成的点的个数。

在这个例子下，像例1一样简单应用union和unique就不行了。

通常我们考虑用循环：

p=unidrnd（1000,100000,2）;

tic;

A=zeros（1000）;

fork=1:

100000

A（p（k,1）,p（k,2））=A（p（k,1）,p（k,2））+1;

end

toc

Elapsedtimeis0.399575seconds.

而用accumarray,可以这样：

>>tic;a=accumarray（p,1,[1000,1000]）;toc

Elapsedtimeis0.038169seconds.

>>isequal（A,a）

ans=

1

速度整整提高了10倍。

例3:

1000人，身高分布在170-180，体重在110-200斤，年龄分布在20-50岁，计算身高体重都相等的人的年龄平均值。

结果用矩阵表示：

行数表示身高，列数表示体重，元素表示年龄的平均值。

首先生成数据：

rand（'state',0）

height=unidrnd（10,1000,1）+170;

rand（'state',0）

weight=unidrnd（90,1000,1）+110;

rand（'state',0）

old=unidrnd（30,1000,1）+20;

利用accumarray计算的语句如下：

tic;mo=accumarray（[height,weight],old,[],@mean）;toc

Elapsedtimeis0.005170seconds.

这个矩阵比较稀疏也可以结果用稀疏矩阵来表示

tic;mo=accumarray（[height,weight],old,[],@mean,0,true）;toc

Elapsedtimeis0.005312seconds.

维数大后，稀疏矩阵会有优势，这个例子还不明显。

大家有兴趣可以试试传统方法效果。

以上仅举了三个例子，实际上，accumaary的应用方法非常灵活，尤其是对于很多要操作大矩阵的情况下。

大家可以仔细看看帮助文档