提高matlab运算速度的几种方法

由于matlab是一种解释性语言，所以在matlab程序中最忌讳直接使用循环语句，如果不得已要使用for循环，可以采用以下方法提高速度。

1、使用6.5以上版本，对循环已作优化；
2、尽可能转化为矩阵运算；
3、转化为二进制执行文件运算，如使用matlab内带的编译系统或matcom以及com组件技术。

其中com组件技术最方便的就是利用com builder来实现，这里重点介绍。

com builder是matlab6.5才有的，也是mathworks公司推荐使用于混合编程的，这些日子进行了全方位的摸索，感觉是爽呆了，下面我们一起来揭开它的神秘面纱。

此系列分为以下几块：
1.matlab下做com组件
2.vb,c#.net实现调用
3.vc实现调用
4.打包
5.优缺点评注
其中2,3部分可以选择一个看
matlab下做com组件
com是component object module的简称，它是一种通用的对象接口，任何语言只要按照这种接口标准，就可以实现调用它。

matlab6.5新推出来的combuilder就是把matlab下的程序做成com组件，供其他语言调用。

我们先准备两个测试文件，并copy一个图片到c盘下，起名叫1.jpg(这些你都可以改的，我这儿是为了程序方便)：
第一个叫im_test.m如下:
function im_test %这个文件不带输入与输出
I=imread('c:\1.jpg'); %因为以前带有imshow的程序用mcc编成dll后用不
%了，所以想试combuilder是否
imshow(I); %能支持这些函数
第二个叫split2rgb.m，就是前些日子Zosco发给我的那个程序，因为它用mcc编成dll后有问题，所以我在这儿继续将它进行测试，而且它也带有多个输入及输出参数，所以也正好拿来测试
在matlab的workspace下打comtool，就打开了matlab com builder，点击file-new project，新建一个工程，在component name里填上comtest，Class name里填上一个sgltest(并将自动生成classes里的comtest remove掉),complie
code in选c或c++都无所谓，将Complier options里的Use Handle Graphics library
的复选框画上，点击ok就行了。

然后点击project--Add files,将im_test.m和split2rgb.m 添加入工程，然后点Build-Com Object，就会在comtest\distrib\文件夹下生成一个comtest_1_0.dll(它就是做好的com组件)，Build时matlab已经自动将此dll在注册表中注册,为了下面能用其他编译器调用，我们还需做一个准备工作，开一个dos窗口，进入<matlabroot>/bin/win32目录下(matlabroot为你机器上matlab安装的路径),打regsvr32 mwcomutil.dll，即对mwcomutil.dll进行注册(这个dll是matlab下作的任何com组件都要用到的dll),下面我们在其他编译器下调用时就可以用了。

是不是觉得做起com组件来很简单呢，matlab下可以给com组件中的类添加成员、事件、方法等，我这儿其实是给类sgltest添加了两个方法，怎么添加成员和方法可以参看matlab 的com builder的帮助。

附录：split2rgb.m的源代码
%%// 测试文件
function [m_nHeight,m_nWidth,mOrigR,mOrigG,mOrigB]=Split2RGB(FileName)
%%// 读入一个Jpg文件，
mOrigData=imread(FileName);
%mDestData=imresize(mOrigData,0.5);
imwrite(mOrigData,'c:\2.jpg');
%%// 用三个变量保存其R,G,B分量
mOrigR=mOrigData(:,:,1);
mOrigG=mOrigData(:,:,2);
mOrigB=mOrigData(:,:,3);
%%// 获得图象的高度，宽度
[m_nHeight,m_nWidth]=size(mOrigR);
figure(1);
set(gcf,'MenuBar','none')
imshow(mOrigData);
title(['Orginal Image:',FileName],'Color','b','FontSize',14);
xlabel(['Height: ',num2str(m_nHeight),' Width :',num2str(m_nWidth)],'Color'
,'b','FontSize',12);
%%// 写param文件
paraFName=[FileName(1:length(FileName)-4),'.param'];
fid=fopen(paraFName,'w');
fwrite(fid,m_nHeight,'uint32');
fwrite(fid,m_nWidth,'uint32');
fclose(fid);
%%// 写 R 分量文件
RFName=[FileName(1:length(FileName)-4),'_R.rot'];
fid=fopen(RFName,'w');
fwrite(fid,mOrigR,'uint8');
fclose(fid);
%%// 写 G 分量文件
GFName=[FileName(1:length(FileName)-4),'_G.rot'];
fid=fopen(GFName,'w');
fwrite(fid,mOrigG,'uint8');
fclose(fid);
%%// 写 B 分量文件
BFName=[FileName(1:length(FileName)-4),'_B.rot'];
fid=fopen(BFName,'w');
fwrite(fid,mOrigB,'uint8');
fclose(fid);
这一部分讲vb,c#.net下怎么实现调用上一部分我们生成的comtest_1_0.dll(matlab下做的com组件),记得一定先要对mwcomutil.dll进行注册(怎么注册参看上一部分)
1.vb下实现调用
打开或新建一个vb工程，点project-Reference,在弹出来的窗口中找到comtest 1.0 Type Library,将前面的复选框选上，点击ok，此时便将此com组件添加到工程里面去了，此时你可以用对象浏览器看到comtest下有个sgltest类，这个类里面有两个方法im_test,split2rgb,还有个MWFlags成员(这个成员是自动生成的),vb下测试代码如下:
测试im_test方法的代码:
Dim st As sgltest
Set st = New sgltest
Call st.im_test
测试split2rgb方法的代码:
Dim st As sgltest
Set st = New sgltest
Dim h As Variant, w As Variant, r As Variant, g As Variant, b As Variant, filename As Variant
filename = "c:\\1.jpg"
Call st.split2rgb(5, h, w, r, g, b, filename)
可见matlab下函数的输入输出参数在com组件里全是variant型的变量，测试大获成功，结果就跟matlab下运行的一摸一样，爽
2.c#.net下实现调用
打开或新建一个c#项目(我采用的是编辑器)，选中右边的解决方案资源管理器中的引用，点鼠标右键，选添加引用，在弹出来的窗口中选com，然后也找到comtest_1_0.dll，点选择，然后确定就可，此时此com组件也添加到工程里面去了，同样我们可以选择对象浏览器来看comtest及下面的sgltest类，c#测试项目如下:
测试im_test方法的代码:
comtest.sgltestClass st=new comtest.sgltestClass();
st.im_test();
测试split2rgb方法的代码:
comtest.sgltestClass st=new comtest.sgltestClass();
object h=null,w=null,r=null,g=null,b=null;
object filename="c:\\1.jpg";
st.split2rgb(5,ref h,ref w,ref r,ref g,ref b,filename);
可见输入参数是ref object型的，而输出参数是object型的，测试同样大获成功，与matlab 下运行的结果一摸一样，爽呆了。

这一部分讲vc下实现调用第一部分我们生成的comtest_1_0.dll，同样要记得先对mwcomutil.dll进行注册(怎么注册参看第一部分)，vb和.net下实现对com组件的调用很简单，而vc下实现这一步我可是摸索了好几天(主要是vc用的不好)
1.先做一些准备工作，用ole viewer工具开始--程序--Microsoft visual studio6.0--Microsoft visual studio6.0 Tools--OLE viewer(这个工具也可以通过在vc 下点Tools--OLE/COM Object Viewer来打开，在Oleviewer工具里，在右边选择Type libraries,将他展开，找到comtest 1.0 Type Library，选中它，点鼠标右键，选view，便又弹出一窗口，点工具栏上的save按钮，分别将他保存为comtest_1_0.h,comtest_1_0.c(也可以存为comtest_1_0.Idl接口文件),我们就可以通过这两个文件实现对comtest_1_0.dll调用
2.vc下调用
新建或打开一个vc工程，注意，此时不用对编译器进行任何设置(而用mcc生成的dll我们么设置一大堆,我这儿只设置了Precomplied headers,选Automatic use of precompliedheaders,写上 stdafx.h),将comtest_1_0.h和comtest_1_0.c加入工程，并复制一个comtest_1_0.dll到工程目录下，由于comtest_1_0.dll还要用到mwcomutil.dll，所以将<matlabroot>/extern/include/下的mwcomutil.h和mwcomtypes.h也加入工程(将这两个文件拷贝入工程路径下，如果设置了library path,可以不加)此时可以通过classView 看到多出了_IID、IMWUtil,Isgltest类，Isgltest就是我们在matlab下建起来的sgltest 类
vc下代码如下
//这几个是引入dll和头文件
#import "mwcomutil.dll"
#import "comtest_1_0.DLL"
#include "mwcomutil.h"
#include "comtest_1_0.h"
#include "comutil.h" //此文件是用来处理由char *向VARIANT类型的转换
测试im_test方法的代码:
if(FAILED(CoInitialize(NULL))) //初始化调用com
{
AfxMessageBox("unable to initialize COM");
}
Isgltest *st=NULL;
//创建一个com组件，CLSID_sgltest和IID_Isgltest可以从comtest_1_0.h和comtest_ 1_0.c里找到
HRESULT hr=CoCreateInstance(CLSID_sgltest,NULL,CLSCTX_ALL,IID_Isgltest,(voi
d **)&st);
if(SUCCEEDED(hr))
{
st->im_test();
AfxMessageBox("succeed");
st->Release();
}
else
{
AfxMessageBox("unsucceed");
}
如果你的vc工程是console project的话，上述的AfxMessageBox可改成printf或cout，测试split2rgb方法的代码(类型转换我参照visual c的精华区也转换成功了)
if(FAILED(CoInitialize(NULL)))
{
AfxMessageBox("unable to initialize COM");
}
Isgltest *st=NULL;
HRESULT hr=CoCreateInstance(CLSID_sgltest,NULL,CLSCTX_ALL,IID_Isgltest,(voi
d **)&st);
VARIANT m,n,r,g,b,filename;
VariantInit(&m);
VariantInit(&n);
VariantInit(&r);
VariantInit(&g);
VariantInit(&b);
VariantInit(&filename);
filename.vt=VT_BSTR;
filename.bstrVal=_com_util::ConvertStringToBSTR("C:\\1.jpg");
if(SUCCEEDED(hr))
{
st->split2rgb(5,&m,&n,&r,&g,&b,filename);
st->Release();
AfxMessageBox("succeed");
}
else
{
AfxMessageBox("unsucceed");
}
同样，运行结果与matlab下的结果一摸一样，记得我们的im_test里面可是使用了imshow 阿，以前用mcc生成的程序中用它可是有错哦，爽呆了。

关于VC下用com组件及其类型的转变请参看msdn及其Visual C的精华区
combuilder系列可以结尾了
一.打包:
在matlab下的workspace里打comtool，点file-open project将我们先前建好的comtest.cbl工程文件打开，再点component--package component就实现了打包，此时到comtest\distrib文件夹里看，生成的comtest.exe就是打包后的解压程序，双击它会解压出一些文件，再点击解压出来的_install.bat就可以实现安装(按他的步骤去做就行了)
二.优缺点评注
这几天用这个combuilder可把我给爽死了，特别是在vc下调用成功时，记得精华区里曾讲combuilder没有什么实质性的突破，我可不这么认为，它的突破可大了
1.做出来的是com，通用的，任何编译器只要支持com,就可以实现调用，想c++ builder,Delphi等的，我想只要按照调用com组件去做，也能成功的
2.支持imshow等一些原来混编用不了的函数，对图形库的支持也比以前强(这些还需各位大侠共同测试)
3.实现方法简单，没有像以前混编还要设置一大堆东东
4.能够在matlab下写自己的类，并为自己的类编写成员、方法和事件，管理工程也方便(这个有点像vc、vb下管理工程一样)用的太爽了，一下子还不知道怎么写缺点了，^0^，我想缺点还需大家一起用来找出我这儿说一个缺点，感觉它的参数全是VARIANT型的，不怎么方便。

一、基本技术
-----------------------------------------------------
1)MATLAB索引或引用(MATLAB Indexing or Referencing)
在MATLAB中有三种基本方法可以选取一个矩阵的子阵。

它们分别是
下标法，线性法和逻辑法(subscripted, linear, and logical)。

如果你已经熟悉这个内容，请跳过本节
1.1)下标法
非常简单，看几个例子就好。

A = 6:12;
A([3,5])
ans =
8 10
A([3:2:end])
ans =
8 10 12
A = [11 14 17; ...
12 15 18; ...
13 16 19];
A(2:3,2)
ans =
15
16
1.2)线性法
二维矩阵以列优先顺序可以线性展开，可以通过现行展开后的元素序号来访问元素。

A = [11 14 17; ...
12 15 18; ...
13 16 19];
A(6)
ans =
16
A([3,1,8])
ans =
13 11 18
A([3;1;8])
ans =
13
11
18
1.3)逻辑法
用一个和原矩阵具有相同尺寸的0-1矩阵，可以索引元素。

在某个
位置上为1表示选取元素，否则不选。

得到的结果是一个向量。

A = 6:10;
A(logical([0 0 1 0 1]))
ans =
8 10
A = [1 2
3 4];
B = [1 0 0 1];
A(logical(B))
ans =
1 4
-----------------------------------------------------
2)数组操作和矩阵操作(Array Operations vs. Matrix Operations)
对矩阵的元素一个一个孤立进行的操作称作数组操作；而把矩阵视为一个整体进行的运算则成为矩阵操作。

MATLAB运算符*,/,\,^都是矩阵运算，而相应的数组操作则是.*, ./, .\, .^
A=[1 0 ;0 1];
B=[0 1 ;1 0];
A*B % 矩阵乘法
ans =
0 1
1 0
A.*B % A和B对应项相乘
ans =
0 0
0 0
------------------------------------------------------
3)布朗数组操作(Boolean Array Operations)
对矩阵的比较运算是数组操作，也就是说，是对每个元素孤立进行的。

因此其结果就不是一个“真”或者“假”，而是一堆“真假”。

这个
结果就是布朗数组。

D = [-0.2 1.0 1.5 3.0 -1.0 4.2 3.14];
D >= 0
ans =
0 1 1 1 0 1 1
如果想选出D中的正元素：
D = D(D>0)
D =
1.0000 1.5000 3.0000 4.2000 3.1400
除此之外，MATLAB运算中会出现NaN,Inf,-Inf。

对它们的比较参见下例Inf==Inf返回真
Inf<1返回假
NaN==NaN返回假
同时，可以用isinf，isnan判断，用法可以顾名思义。

在比较两个矩阵大小时，矩阵必须具有相同的尺寸，否则会报错。

这是你用的上size和isequal,isequalwithequalnans(R13及以后)。

------------------------------------------------------
4)从向量构建矩阵(Constructing Matrices from Vectors)
在MATLAB中创建常数矩阵非常简单，大家经常使用的是：
A = ones(5,5)*10
但你是否知道，这个乘法是不必要的？
A = 10;
A = A(ones(5,5))
A =
10 10 10 10 10
10 10 10 10 10
10 10 10 10 10
10 10 10 10 10
10 10 10 10 10
类似的例子还有：
v = (1:5)';
n = 3;
M = v(:,ones(n,1))
M =
1 1 1
2 2 2
3 3 3
4 4 4
5 5 5
事实上，上述过程还有一种更加容易理解的实现方法：
A = repmat(10,[5 5]);
M = repmat([1:5]', [1,3]);
其中repmat的含义是把一个矩阵重复平铺，生成较大矩阵。

更多详细情况，参见函数repmat和meshgrid。

-----------------------------------------------------
5)相关函数列表(Utility Functions)
ones 全1矩阵
zeros 全0矩阵
reshape 修改矩阵形状
repmat 矩阵平铺
meshgrid 3维plot需要用到的X－Y网格矩阵
ndgrid n维plot需要用到的X－Y－Z...网格矩阵
filter 一维数字滤波器，当数组元素前后相关时特别有用。

cumsum 数组元素的逐步累计
cumprod 数组元素的逐步累计
eye 单位矩阵
diag 生成对角矩阵或者求矩阵对角线
spdiags 稀疏对角矩阵
gallery 不同类型矩阵库
pascal Pascal 矩阵
hankel Hankel 矩阵
toeplitz Toeplitz 矩阵
=================================================== =======
二、扩充的例子
------------------------------------------------------
6)作用于两个向量的矩阵函数
假设我们要计算两个变量的函数F
F(x,y) = x*exp(-x^2 - y^2)
我们有一系列x值，保存在x向量中，同时我们还有一系列y值。

我们要对向量x上的每个点和向量y上的每个点计算F值。

换句话
说，我们要计算对于给定向量x和y的所确定的网格上的F值。

使用meshgrid，我们可以复制x和y来建立合适的输入向量。

然后可以使用第2节中的方法来计算这个函数。

x = (-2:.2:2);
y = (-1.5:.2:1.5)';
[X,Y] = meshgrid(x, y);
F = X .* exp(-X.^2 - Y.^2);
如果函数F具有某些性质，你甚至可以不用meshgrid，比如
F(x,y) = x*y ，则可以直接用向量外积
x = (-2:2);
y = (-1.5:.5:1.5);
x'*y
在用两个向量建立矩阵时，在有些情况下，稀疏矩阵可以更加有
效地利用存储空间，并实现有效的算法。

我们将在第8节中以一个实例来进行更详细地讨论.
--------------------------------------------------------
7)排序、设置和计数(Ordering, Setting, and Counting Operations) 在迄今为止讨论过的例子中，对向量中一个元素的计算都是独立
于同一向量的其他元素的。

但是，在许多应用中，你要做的计算
则可能与其它元素密切相关。

例如，假设你用一个向量x来表示一个集合。

不观察向量的其他元素，你并不知道某个元素是不是一
个冗余元素，并应该被去掉。

如何在不使用循环语句的情况下删除冗余元素，至少在现在，并不是一个明显可以解决的问题。

解决这类问题需要相当的智巧。

以下介绍一些可用的基本工具
max 最大元素
min 最小元素
sort 递增排序
unique 寻找集合中互异元素（去掉相同元素）
diff 差分运算符[X(2) - X(1), X(3) - X(2), ... X(n) - X(n-1)]
find 查找非零、非NaN元素的索引值
union 集合并
intersect 集合交
setdiff 集合差
setxor 集合异或
继续我们的实例，消除向量中的多余元素。

注意：一旦向量排序后，任何多余的元素就是相邻的了。

同时，在任何相等的相邻元素在向量diff运算时变为零。

这是我们能够应用以下策略达到目的。

我们现在在已排序向量中，选取那些差分非零的元素。

% 初次尝试。

不太正确!
x = sort(x();
difference = diff(x);
y = x(difference~=0);
这离正确结果很近了，但是我们忘了diff函数返回向量的元素个数比
输入向量少1。

在我们的初次尝试中，没有考虑到最后一个元素也可能
是相异的。

为了解决这个问题，我们可以在进行差分之前给向量x加入
一个元素，并且使得它与以前的元素一定不同。

一种实现的方法是增
加一个NaN。

% 最终的版本。

x = sort(x();
difference = diff([x;NaN]);
y = x(difference~=0);
我们使用(运算来保证x是一个向量。

我们使用~=0运算，而不用find 函数，因为find函数不返回NaN元素的索引值，而我们操作中差分的最后元素一定是NaN。

这一实例还有另一种实现方式：
y=unique(x);
后者当然很简单，但是前者作为一个练习并非无用，它是为了练习使用
矢量化技术，并示范如何编写你自己的高效代码。

此外，前者还有一个
作用：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
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
/company/digest/sept01/matrix.shtml
下面的说明链接将指导你如何使用MATLAB中的矩阵操作。

含括矩阵创建
、索引、操作、数组运算、矩阵运算以及其它主题。

MATLAB Documentation: Getting Started
/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?
/support/solutions/data/26623.shtml
下面的技术手册是关于MATLAB存储管理方面的一个包罗广泛的指南：
The Technical Support Guide to Memory Management
/support/tech-notes/1100/1106.shtml
----------------------------------------------------------------
12)发挥MATLAB的最高性能(Maximizing MATLAB Performance)
这个技术手册对代码矢量化进行一般的讨论，而许多时候，我们的目的是加快代码的速度。

因此，本节提供一些附加的资源，它们涉及如何是代码达到最高性能。

加速代码的常用技巧：
3257: How do I increase the speed or performance of MATLAB?
/support/solutions/data/3257.shtml
编译你的代码并从MATLAB中调用它
22512: What are the advantages of using the MATLAB Compiler to
compile my M-files?
/support/solutions/data/22512.shtml
=================================================== ===================
结论
正如本文所介绍的，针对代码矢量化和代码提速，有一些常规方法，同时，还有
性能。

加速代码的常用技巧：
3257: How do I increase the speed or performance of MATLAB?
/support/solutions/data/3257.shtml
编译你的代码并从MATLAB中调用它
22512: What are the advantages of using the MATLAB Compiler to
compile my M-files?
/support/solutions/data/22512.shtml
=================================================== ===================
结论
正如本文所介绍的，针对代码矢量化和代码提速，有一些常规方法，同时，还有
一些非常规方法，正如第9节最后那个例子(对sparse函数进行独特利用)那样。

以上的讨论远没有穷尽MATLAB中所有矢量化技术。

它仅仅是代码矢量化的思想和理论的一个导引。

根据我的一些经验和相关的参考资料，要想提高其运算的速度，以下的方法可以试一下：
（1）尽可能地用向量化的数组运算代替循环，尽可能地减少使用户for或while 循环，这是因为matlab执行循环运算效率很低而数组运算效率较高。

举个最经典的例子，下面的循环：
k=0;
for t=0:pi/20:2*pi;
k=k+1;
y(k)=cos(t);end;
就可以用t=0:pi/20:2*pi;y=cos(t)代替，这样即减少代码量又能提高运算速度(2)如果非要用循环的话尽可能进行循环内数组的预配置而不是让程序在循环中不断地动态配置，好比像C语言中要想使用数组必须先定义数组的长度一样，比如说如果想要做以下的循环：
k=0;
for t=1:1:20;
k=k+1;
y(k)=t^2;end;
可以事先定义y=linspace(0,0,20)即定义y为一个拥有20个0的一维数组，这样就相当于是事先对循环内的数组进行了预配置；
（3）尽可能采用MATLAB自带的函数指令，这些函数都是一些很经典的算法构成的，比如说要想找最大值就应用函数max而不是自己去编程序，这样只会事倍功半；
（4）尽量采用M函数文件替代M的脚本文件，因为函数文件运行时是采用P 码方式驻留在内存中而不是像脚本文件一样每运行一次都要经历把程序装入内存的过程，因而比较省时；
（5）在循环比较费时的时候可以考虑采用非解释执行的MEX文件来对此进行表达；
（6）尽可能找出导致程序运行缓慢的瓶颈，可以在MATLAB的View:profiler 中打开程序剖析器找出程序运行中的瓶颈，具体使用的方法请参考MATLAB的帮助或者在网上寻找；
（7）从硬件来看电脑的CPU或者内存（不是大小而是读写速度）可能会是瓶颈，在做大型运算时可以考虑较好的配置。

希望对你有所帮助！。