【转载】让你的MATLAB代码飞起来

原文地址:http://developer.51cto.com/art/201104/255128_all.htm

MATLAB语言是一种被称为是“演算纸”式的语言,因此追求的是方便性、灵活性以及交互性。在快速性上要比C语言这种性能强劲著称的稍逊一筹。然而,通过一些手段,我们也能让MATLAB语言快起来,甚至和C差不多了!

MATLAB语言是一种被称为是“演算纸”式的语言,因此追求的是方便性、灵活性以及交互性。在快速性上要比C语言这种性能强劲著称的稍逊一筹。然而,通过一些手段,我们也能让MATLAB语言快起来,甚至和C差不多了!

首先声明:本文是一个初级教程,因此很多知识是假定你已经很熟悉了的;虽然我在讨论让代码飞起来,但从来不会说最快有多快,究竟有多快你要自己感觉;作者水平不是很高,难免误导你,小心甄别。



在正式讨论之前,先看看这些好习惯你有没有?

1. 使用 M-Lint

M-Lint是一个代码分析检查工具,它在你写代码的过程中实时交互,发现你代码的问题,按照最佳性能和最可维护性给出修改建议。

注意:我可没说是最正确!

如果没有激活这个功能,依次使用File > Preferences > M-Lint,勾选Enable integrated M-Lint warning and error messages 。同时,还可以设定你的偏好。

激活后,在你写代码时就会实时交互了,错误的或者不推荐的部分会以红色下划线标出,鼠标经过红色下划线的语句或单词,M-Lint给出提示信息。想一下子看遍全部提示信息。使用Tools >M-Lint > (Save and) Show M-Lint Report2。

注:首次“观看”先提示先保存一下。

2. 组织

给每一个项目(project)建立一个单独的文件夹。同属于一个项目的文件保存在哪儿的都有,你找的时候就不费劲吗!

写头部注释,尤其是H1。第一行就是H1。MATLAB中的内置函数的 help的内容其实就是读取的这个函数的头部注释。怎么写,参照MATLAB内置函数。

将经常用到的控制台命令存储为脚本(script)。如果有些命令反复使用,还是存为脚本吧,没别的意思,你要少敲多少次键盘啊!

3. 避免数据丢失

不要在脚本中使用 clear all。不幸的是这是一个大家常用的命令,有些书上还作为一条规则确立起来,建议必须使用!要知道这个命令一执行,工作空间的数据可就不可逆转的全没了啊!

警告:注意呦! ?

小心同名覆盖。如果你在一个文件中,本来你的意思是两个变量,你却给他们起了相同的名字,那么第一次的数据可就没了。比如:

  1. result=max(a,b); %想求 a和 b之较大者
  2. result=max(c,d); %想求 c 和d之较大者

result结果是什么?恐怕不是你想要的。不妨将其改为result1和 result2。类似的,也要小心文件重名的覆盖,这个后果貌似更严重些。

下条内容请重视!

如何让 MATLAB崩溃。

尽管 MATLAB是很稳定的,但是我们仍然可以让它崩溃!使用第三方的MEX函数或者耗内存的操作比如视频处理或者超大规模矩阵都可能会造成MATLAB崩溃。



如果你已经有这些好习惯,那么恭喜,你要是还有其他好习惯麻烦也告诉我一声!如果没有,相信你看完之后总该有了吧?好了,我们开始!

1.使用profile

profile,Longman 给出的解释是:a short description that gives important details about a person, a group of people, or a place。

MATLAB中内置了一个叫做profile的工具,来协助评估程序,也就是对程序运行过程的一个short description吧。主要命令有:

profile on 开启

profile off 关闭

profile clear 清空数据

profile viewer 在profiler中看结果

下面我们评估一下下面这个函数:

  1. function result =example1(count)
  2. for k = 1:count
  3. result(k) = sin(k/50);
  4. if result(k)<-0.9
  5. result(k) = gammaln(k);
  6. end
  7. end

为了分析这个函数的效率,首先开启并清空 profiler,然后运行这个函数,接下来看结果报告。即依次输入:

  1. >> profile on, profile clear
  2. >> example1(5000);
  3. >> profile viewer

这就是 profile 的基本语法。也有使用鼠标操作的方法,这里就不介绍了,那样虽然直观单远不及使用,命令方便。

由于系统的不同,报告的结果一般是不一样的。以下是我的系统得出的结果。

1.先看profile summary:

2.点击example1链接,进入具体各小项的评估。

(1)调用函数(children)、被调用函数(parents)。本例中都没有。如果被 profile 的对象有调用函数或者被调用函数的话,会给出相应的数据。

(2)时间在哪些行被消耗(Lines where the most time was spent):

从数据中我们可以看出哪些行消耗了多少时间(总时间和相对时间),被调用了多少次,以及直观的柱形图。

(3)另一个有用的项目是 M-Lint 结果,给出了错误(警告、提示)所在的行,以及对应的建议修改信息,这些建议对代码的改进是很有价值的信息:

(4)最下面还有一个函数列表,是(2)的另一种形式。看图:

最右侧是函数代码,前有行号、每一行调用的次数和小号的时间。消耗时间最多的行被标示了出来。最红的消耗时间最多。

profiler工具的时间分辨率不是很高,因此,如果你的代码运行的时间很短,有时候恐怕不能感知到。这时候不妨人为的加入几个循环,让程序所运行几次,然后进行分析。

必须指出,profile工具的作用主要是分析程序,获得程序运行的信息。如果想要知道程序运行的精确时间,使用计时器 tic/toc。以上面程序为例,在命令行输入:

  1. >> tic;example1(5000);toc

输出是:

  1. Elapsed time is 0.058522 seconds.

为了获得更为精准的结果,你最好把浏览器、杀毒软件、防火墙等等占用CPU时间片的程序先关了,只剩下不能关掉的系统进程。

注意:profile在新版本中不断被加强,可使用的参数也越来越多,不过大多数根本用不着,如果你觉得那些参数很有用,我相信你根本用不找看我这个小册子了,要真是这样,麻烦您不吝赐教,分享一些经验。更详细的内容,您还是去看文档去吧!

2. 预分配矩阵

MATLAB中的矩阵变量可以动态增长行和列。比如:

  1. >>x=2
  2. x=
  3. 2
  4. >>x(2,3)=1
  5. x=
  6. 2   0   0
  7. 0   0   1

看到没?MATLAB自动调整了矩阵的大小!从内部实现上看,矩阵数据存储单元被重新分配了更大的单元。如果矩阵的大小被反复的调整(比如在循环中),重新分配存储空间带来的额外开销会是很显著的。为了避免反复的矩阵存储重新分配,预分配矩阵的存储单元是一个不错的选择。一个推荐的方法是使用 zeros 函数命令。看下面的代码:

  1. a(1) = 1;
  2. b(1) = 0;
  3. for k = 2:8000
  4. a(k) = 0.99803 * a(k-1)-0.06279 * b(k-1);
  5. b(k) = 0.06279 * a(k-1) + 0.99803 * b(k-1);
  6. end
  7. tic/toc计时运行得到:
  8. Elapsed time is 0.013306 seconds.

简单分析上面的代码,知道,每一次 for,矩阵 a 和 b 的大小都要被重新分配,最终的大小事 8000 的列向量。如果我们提前就给它们分配好大小为 8000的存储空间,看看结果怎么样:

  1. a=zeros(1,8000);    %预分配矩阵存储单元
  2. b=zeros(1,8000);
  3. a(1) = 1;
  4. b(1) = 0;
  5. for k = 2:8000
  6. a(k) = 0.99803 * a(k-1)-0.06279 * b(k-1);
  7. b(k) = 0.06279 * a(k-1) + 0.99803 * b(k-1);
  8. end
  9. 及时运行得到:
  10. Elapsed time is 0.000753 seconds.

看出来没?速度提高了近 18 倍!像这种只需添加几行代码就能做到的情况是很多的。这个例子也有特殊性,就是最后的结果大小已知,如果结果的大小可变、未知呢?没关系,我们可以估计一下,最终结果最大能是多少?比估计到的最大再留出一些余量就成了!如果你估计的还是不够大,那超出的部分还要反复重新分配,不过这样节省下来的时间也是很可观的,毕竟可以少分配很多次了! 最后呢,还要处理一下后事,比如你分配给变量 a 有 1000 个单元,但最终它只占了300个,那你还要将那700个给收回来。看下面的代码:

  1. a = zeros(1,10000);     %  预分配
  2. count = 0;
  3. for k = 1:10000
  4. v = exp(rand*rand);
  5. if v > 0.5        %  增长结果不确定的来源
  6. count = count + 1; a(count) = v;
  7. end
  8. end
  9. a = a(1:count);    %调整矩阵大小
  10. 未预分配时:Elapsed time is 0.052395 seconds.
  11. 预分配后:Elapsed time is 0.008935 seconds.

感慨:些微时间的意义在哪呢?背后是你对 MATLAB 的理解深度。哥玩的不是时间,是技术。

3. 向量化

很多情况下,程序中的某些代码可以被向量化,向量化前后的速度往往在10 倍以上!向量化是最基本和最有效的让代码快起来的技巧,我都不愿意在后面叫“之一”了。

(1)向量化的计算

很多常规函数都是向量化的,它们作用于数组时,就好像是作用于数组中的每一个元素。例如:

  1. >> sqrt([1,4,9,16])
  2. ans =
  3. 1     2     3     4
  4. 考虑下面的函数:
  5. function d = minDistance(x,y,z)    %寻找点集中距离远点最近点
  6. nPoints = length(x);
  7. d = zeros(nPoints,1);    %  预分配
  8. for k = 1:nPoints    %  计算每一个点的距离
  9. d(k) = sqrt(x(k)^2 + y(k)^2 + z(k)^2);
  10. end
  11. d = min(d);    %  得到最小距离
  12. 取  x=[1 2 3 4 5 6]; y=[2 3 5 2 1 4];z=[9 2 3 2 1 5];
  13. 计时运行:Elapsed time is 0.008006 seconds.

如果你写出上面类似的代码,说明你认真看了前面的内容。为d预分配空间确实为本例节省了不少时间。如果采用向量化计算,我们可以去掉for循环,直接计算向量。这里要隆重推出“.”运算符,它表示的是对应元素进行运算。有.*和./和.\和.‘和.^等。分别表示不带.运算的对应元素运算。假设A是方阵,A^2是矩阵的 2 次乘幂,而 A.^2 表示矩阵 A 中的元素各自求平方组成新的矩阵。考虑下面的代码:

  1. function d = minDistance(x,y,z)
  2. d = sqrt(x.^2 + y.^2 + z.^2);    %  计算每一点的距离
  3. d = min(d);
  4. 计时运行:
  5. Elapsed time is 0.005326 seconds.

貌似差别不大?这就对了,别忘了,咱可就计算了6个值啊!这么几个值就有了这样的差距,那x、y、z向量要是大一点,结果的差异就可想而知了!

更进一步的,我们可以使用d = sqrt(min(x.^2 + y.^2 + z.^2))取代后两行语句,让程序更加简洁。

一下函数使用向量化的计算会更为节省时间:min, max, repmat, meshgrid,sum, diff, prod等等。

(2)向量化逻辑

上面讨论了计算的向量化,其实MATLAB的逻辑运算也是向量化的。比如:

  1. >> [1 4 2]>[2 3 1]
  2. ans =
  3. 0     1     1

两个数组“按元素”进行比较。向量的逻辑操作返回二进制的逻辑结果向量,即用0代表假,用1代表真。这为什么有用呢?因为MATLAB中有一些强劲的针对逻辑向量的函数。例如:

  1. >> find([1,5,3] < [2,2,4])
  2. ans =
  3. 1      3
  4. >> any([1,5,3] < [2,2,4])
  5. ans =
  6. 1
  7. >> all([1,5,3] < [2,2,4])
  8. ans =
  9. 0

其实,对一般向量(非逻辑向量)也是适用的!

  1. >> find(eye(4)==1)
  2. ans =
  3. 1
  4. 6
  5. 11
  6. 16

以上函数的用法请自己查阅函数说明。

4. 示例

(1)向量归一标准化

将一个向量v归一标准化,我们可是使用v = v/norm(v),norm函数的作用是求模(范数)。

如果对一组向量 v(:,1), v(:,2),…进行归一标准化,可以使用一个循环计算v(:,k)/norm(v(:,k))。更好的策略是向量化计算:

  1. vMag = sqrt(sum(v.?2));
  2. v = v./vMag(ones(1,size(v,1)),:);

(2)剔除元素

有时候,我们需要将矩阵中的符合某些条件的元素剔除,当然可以使用条件判断加循环。我们使用向量化剔除矩阵中的NaN和无穷两类数:

  1. i = find(isnan(x) | isinf(x));   %在x中找到符合条件的数的位置
  2. x(i) = [];    %剔除它
  3. 或者,同样的功能:
  4. i = find(?isnan(x) & ?isinf(x));     %找到不符合的数
  5. x = x(i);     %保留它
  6. 进一步的,我们可以更加简化,省略中间变量:
  7. x(isnan(x) | isinf(x)) = [];
  8. 以及
  9. x = x(?isnan(x) & ?isinf(x));

(3)分段函数

信号分析中十分重要的 sinc(x)函数是分段的:x=0 时的值是 1,x!=0 时,sinc(x)=sin(x)/x。下面的代码使用向量化方法处理分段:

  1. function y = sinc(x)
  2. y = ones(size(x));          %  先设所有的y都是1
  3. i = find(x ?= 0);          %  找到非零x值
  4. y(i) = sin(x(i)) ./ x(i);      %  计算非零值处的函数值
  5. 更简洁的,可以写成:
  6. y = (sin(x) + (x == 0))./(x + (x == 0))

能看出来吗?里面用到了逻辑运算,实在是巧妙的很!

(4)其他

还有些不常用的,算了,知道也八辈子用不着,珍惜脑细胞吧!

感慨:向量、矢量、相量、复数、数组、矩阵,这些名词能分清楚么?能分清楚知道内涵也就是为什么要这样规定么?不会也别问我!

5. 内嵌简单函数

内嵌函数的意思就是将函数调用的函数的代码直接写到这个函数里面来。由于函数调用要做保护现场以及恢复现场等工作,也会额外增加一些时间消耗。如果调用的次数不是很多,这些时间是可以忽略的,但是当调用次数很多的时候(比如500次),这个时间就很可观了!

什么样的被调用函数适合内嵌呢?正如标题所说,是简单的函数,特征呢就是这个函数只有几行代码。如果这个函数很复杂,代码很长,还是死了这个心吧,内嵌是内嵌了,可是你看不懂代码了,得不偿失。程序的可读性是非常重要的!

注意:必须是 M-File 实现的函数才能内嵌!

下面的代码演示一个反复调用median函数的内嵌方法。原代码:

  1. y = zeros(size(x));      %  预分配
  2. for k = 3:length(x)-2
  3. y(k) = median(x(k-2:k+2));
  4. end
  5. 取  x=rand(1,2500);
  6. 计时运行:Elapsed time is 0.030949 seconds.

下面我们试试内嵌。首先,要研究一下你要内嵌的函数,本例中就是median。在命令行中输入:edit median,发现它是使用sort进行工作的。将核心代码内嵌:

  1. y = zeros(size(x));
  2. for k = 3:length(x)-2
  3. tmp = sort(x(k-2:k+2));
  4. y(k) = tmp(3); ;
  5. end
  6. 仍取x=rand(1,2500);
  7. 计时运行:Elapsed time is 0.011379 seconds.

以上就是一个演示,可见时间确实省去了不少。为了确认你想内嵌的函数是否是用M-File实现的,你可以使用“edit 函数名”命令试试看。

时间: 2024-08-25 18:45:13

【转载】让你的MATLAB代码飞起来的相关文章

Latex 中插入 Matlab 代码

这篇文章将介绍如何在 Latex 排版过程中添加 Matlab 代码 功能效果 主要有如下排版功能: 语法高亮 自动添加边框 自动添加行号 先上图,大家感受一下效果 listings 包 首先确保你能使用使用 listings 包 简单快捷的使用方法如下 \usepackage{listings} \lstset{language=Matlab} \begin{lstlisting} % Plot function f(x) = 2*x^3 - x - 2 ezplot('2*x^3-x-2',

调试和运行matlab代码(源程序)的技巧和教程

转载请标明出处:专注matlab代码下载的网站http://www.downma.com/ 本文主要给大家分享使用matlab编写代码,完成课程设计.毕业设计或者研究项目时,matlab调试程序的技巧和方法. 快速完成一个项目,最简单的方法就是利用前人的开源代码,然后根据自己项目的具体需求和参数,对已有代码进行调试,并增加或删减部分功能,最终实现自己项目的全部功能.所谓"站在前人的肩膀上". 闲话少叙,我们先从最基本的matlab软件安装开始,已经成功安装matlab软件的可以直接跳过

SAR成像学习(四)距离方向成像matlab代码解析 2

如果发射信号是线性调频信号,上一次讲的距离成像算法流程(匹配滤波方法)依然可以用,但那个流程要求T x =4X 0 c >T p  .如果T x <T p  ,即幅宽相对较小的情况,上一讲中的流程会带来一个问题,解决这个问题的办法是pulse compression.本文将会讨论这个puse compression的原理和实现. 1 what is pulse compression 对于线性调频信号:p(t)=a(t)exp(jβt+jαt 2 ) ,信号持续时间为T p  ,瞬时频率为β+

Latex中Matlab代码的环境

需要用到listings宏包 使用方法: 导言区\usepackage{listings}\lstset{language=Matlab}      %代码语言使用的是matlab\lstset{breaklines}                   %自动将长的代码行换行排版\lstset{extendedchars=false}  %解决代码跨页时,章节标题,页眉等汉字不显示的问题 正文区:\begin{lstlisting}此处写下Matlab代码 \end{lstlisting}

多分类问题中,实现不同分类区域颜色填充的MATLAB代码(demo:Random Forest)

之前建立了一个SVM-based Ordinal regression模型,一种特殊的多分类模型,就想通过可视化的方式展示模型分类的效果,对各个分类区域用不同颜色表示.可是,也看了很多代码,但基本都是展示二分类,当扩展成多分类时就会出现问题,所以我的论文最后就只好画了boundary的图了.今天在研究Random Forest时,找到了下面的demo的MATLAB代码,该代码很好的实现了各分类区域的颜色填充,效果非常漂亮. 下面是一个Demo代码:Demo.m %% generate data

Poisson image editing算法实现的Matlab代码解析

之前我发了数篇系列博文来仔细研究Poisson Image Editing算法,每次重新审视和深入,仿佛都能有更为深刻的认识和很大的收获.这应该算是我这个系列的完结篇,会用用Matlab代码一点一点的演示,原文作者到底是如何设计和实现他那个强大且影响深远的算法的.希望你在看本文之前务必参考一下文章来了解算法原理,本文将主要讲解编程实现的问题,对于前面讲过的内容,我不会深究.但我个人总体的感觉是,现在图像处理算法对数学的要求是越来越高了,像泊松融合.泊松抠图这样的算法如果没有偏微分方程(本算法中涉

使用ecilpse(Java)调用Matlab代码

1 安装java环境: http://www.oracle.com/technetwork/java/javase/downloads/index.html 下载JDK最新版本并安装,CloudSim需要运行在jdk1.6以上版本. 以jdk1.6.0_24为例,默认的安装目录为C:\Program Files\Java\jdk1.6.0_24. 设置环境变量: 新建系统变量JAVA_HOME,变量值设为JDK安装目录,即C:\Program Files\Java\jdk1.6.0_24: 在P

hog特征原理详解及matlab代码学习笔记

1.HOG特征: 方向梯度直方图(Histogram of Oriented Gradient, HOG)特征是一种在计算机视觉和图像处理中用来进行物体检测的特征描述子.它通过计算和统计图像局部区域的梯度方向直方图来构成特征.Hog特征结合SVM分类器已经被广泛应用于图像识别中,尤其在行人检测中获得了极大的成功.需要提醒的是,HOG+SVM进行行人检测的方法是法国研究人员Dalal在2005的CVPR上提出的,而如今虽然有很多行人检测算法不断提出,但基本都是以HOG+SVM的思路为主. (1)主

在C#的Web项目中调用Matlab代码的方法

为了毕设的图形检索方向的研究,本人需要在信科的师兄师姐们已经完成的C#界面中,调用现在研究的算法的Matlab代码,以便看到实验的效果.前段时间已经拖延了1个多月,一方面因为实习越来越多事情,时间减少了很多:但更重要在于C#调用Matlab的方法真心麻烦,C#的Web项目中进行这个操作貌似会碰到更多细节上的问题.而且总是很不稳定,操作系统.Matlab或VS的版本.遗漏一些文件或步骤都会造成失败!这个问题本人已经搞了很长时间,直至前几天,在同学的帮助下,自己再弄一遍,总算成功了!下面我及时把这个