精华内容
下载资源
问答
  • 如何提高 matlab 计算速度 运算效率

    万次阅读 2016-10-20 01:10:10
    最近在网上查了一些资料,并结合自己的经验,就如何提高matlab计算效率,总结一下几个原则: 1. 提前给数组分配大小。 2. 尽量用矢量计算,减少 for 循环。 3. 尽量调用 maltab自带的函数来实现一些功能。 4...

    最近在网上查了一些资料,并结合自己的经验,就如何提高matlab计算效率,总结一下几个原则:


    1. 提前给数组分配大小。

    例如:在程序循环时用到数组变量 a ,若知道数组 a 的大小或者知道 a 的大小上界,则可以提前给数组a 分配大小。 一般MATLAB也会给没有提前

    分配大小的变量数组下面标示红线,提醒修改。  

    a=zeros(10,10)
    生成一个10行10列的零矩阵。


    a=100*ones(10,10)
    
    生成一个10行10列的矩阵,并且每个元素都是100。

    a=(1:10)
    生成一个从1到10的1行10列矩阵。


    2. 尽量用矢量计算,减少 for 循环。

    最近我也才意识到这一点,matlab自带的很多函数都支持矢量计算。直接矢量计算,避免大量的for循环。

    用 min 函数举例:

    clc
    
    a=(1:1000000);
    tic
    for i=1:1000000
        a(i)=min(a(i),100);
    end
    toc
    
    tic
    b=max(a,100);
    toc


    结果显示:

    时间已过 0.028604 秒。
    时间已过 0.006296 秒。

    可见后者比前者的 for 循环快。

    3. 尽量调用 maltab自带的函数来实现一些功能。

    这个是显而易见的,matlab自带的函数由大神级的matlab开发人员所编写,他们当然会在计算速度上做不少优化,一定比我们编写的好。


    4. 少用 find 函数,用 logical 替代。

    查找替换矩阵中的元素,过去用 find 函数,现在用 logical 更好。

    下面用一个例子,分别测试 for 循环,find函数,logical函数查找替换:

    a=(1:10000000);
    
    tic
    for i=1:length(a)
        if a(i)>5
            a(i)=1;
        end
    end
    toc
    
    tic
    a(find(a)>5)=1;
    toc
    
    tic
    a(logical(a)>5)=1;
    toc

    结果显示:

    时间已过 0.407840 秒。
    时间已过 0.378585 秒。
    时间已过 0.119918 秒


    可见,find 函数仅仅比 for 循环快一点,而 logical 函数则显著快于它们。 


    5. 必须用到矩阵拼接,可变矩阵时,用 end 拼接矩阵。

    例如:

    function test
    A=cell(3,1);
    tic 
    N=10000;
    for t=1:3
        a=magic(3);
        for i=1:N
            a=[a;magic(3)];
        end
        A{t}=a;
    end
    toc
    
    B=cell(3,1);
    tic
    for t=1:3
        B{t}(end+1:end+3,:)=magic(3);
    end
    toc
    end


    时间已过 1.723104 秒。
    时间已过 0.001621 秒。

    6. 大规模的循环,可以调用 c语言或 C++来 计算。

    这个见不少资料说,自己还没试过,有机会测试一下。

    展开全文
  • 利用MATLAB编写程序,进行二维光栅的衍射级次效率计算
  • 我的数学公式包含很多的未知数,现在这些未知数我要进行逐个求...求出梯度的符号公式后,由于matlab中的sub和eval的效率太低,而我的公式太过于复杂,导致运行速度太慢,matlab能不能解决这个问题,或者用别的语言解决
  • 在使用TreeBagger可能遇到随着数据量以及不同参数设置导致其效率低下的情况,这里将展示如何使用并行计算提升计算速度。 样本数据 样本数据是1985年汽车进口量的数据库,其中有205个样本,25个预测变量和1个因变量。...

    简介

    在使用TreeBagger可能遇到随着数据量以及不同参数设置导致其效率低下的情况,这里将展示如何使用并行计算提升计算速度。

    样本数据

    样本数据是1985年汽车进口量的数据库,其中有205个样本,25个预测变量和1个因变量。前15个变量是数字变量,后10个变量是分类变量。

    代码

    clc
    clear all
    close all
    
    % 加载样本数据,分配到因变量Y和自变量X中。
    load imports-85;
    Y = X(:,1);
    X = X(:,2:end);
    
    % 设置并行环境,使用双核
    mypool = parpool(2)
    % 设置选项以使用并行处理
    paroptions = statset('UseParallel',true);
    
    % 使用计时功能以进行比较
    tic
    b = TreeBagger(5000,X,Y,'Method','r','OOBVarImp','on', ...
        'cat',16:25,'MinLeafSize',1,'Options',paroptions);
    toc
    
    tic
    b = TreeBagger(5000,X,Y,'Method','r','OOBVarImp','on', ...
        'cat',16:25,'MinLeafSize',1);
    toc
    

    效果

    Starting parallel pool (parpool) using the 'local' profile ...
    connected to 2 workers.
    
    mypool = 
    
     Pool - 属性: 
    
                Connected: true
               NumWorkers: 2
                  Cluster: local
            AttachedFiles: {}
        AutoAddClientPath: true
              IdleTimeout: 30 minutes (30 minutes remaining)
              SpmdEnabled: true
    
    时间已过 27.320103 秒。
    时间已过 35.943752 秒。
    
    展开全文
  • 用过Matlab的人都知道,Matlab是一种解释性语言,存在计算速度慢的问题,为了提高程序的运行效率matlab提供了多种实用工具及编码技巧。1. 循环矢量化Matlab是为矢量和矩阵操作而设计的,因此,可以通过矢量化方法...

    用过Matlab的人都知道,Matlab是一种解释性语言,存在计算速度慢的问题,为了提高程序的运行效率,matlab提供了多种实用工具及编码技巧。

    1. 循环矢量化

    Matlab是为矢量和矩阵操作而设计的,因此,可以通过矢量化方法加速M文件的运行。矢量化是指将for循环和while循环转换为等价的矢量或矩阵操作。下面给出一个循环的例子:

    i=0;

    for n = 0:0.1:1000

    i=i+1;

    y(i)=cos(n);

    end

    那么我们可以矢量化为:

    n= 0:0.1:1000;

    y=cos(n);

    我们可以用tic和toc函数来查看上述各代码运行的时间,采用for循环的程序0.39秒(具体时间和计算机配置有关),而矢量化后几乎耗时为0。

    2. 给数组或矩阵预分配内存

    特别是使用大型数组或矩阵时,Matlab进行动态内存分配和取消时,可能会产生内存碎片,这将导致大量闲置内存产生,预分配可通过提前给大型数据结构预约足够空间来避免这个问题。

    3. 用函数代替脚本文件

    因为每次调用MATLAB的脚本文件都需要将不必要的中间变量加载到内存中,每执行一次,就加载一次。函数在调用时被编译成了伪代码,只需要加载到内存一次。当多次调用同一个函数时会运行快一些。因此尽量多使用函数文件而少使用脚本文件,也是提高执行效率的一种方法。

    4. 用Mex文件编写循环代码

    Matlab提供了与C和C++的接口,那么我们可以在用C或C++语言编写耗时的循环代码,然后通过接口程序在Matlab中转换成dll文件,这就是我们所要的Mex文件,通过这种方法可以极大地提高计算速率。

    1.尽量避免使用循环结构

    MATLAB变量的基本类型是矩阵,当对矩阵的每个元素循环处理时,运算速度很慢。因此编程时应尽量把数组和矩阵看作一个整体来进行编程,而不是像其他的程序设计语言那样,使用循环结构对矩阵的元素循环进行处理。利用MATLAB提供的用于矢量化操作的函数,把循环矢量化,这样既可以提高编程效率,也可以提高程序的执行效率。下面给出一个循环的例子:

    i=0;

    for n = 0:0.1:100

    i=i+1;

    y(i)=cos(n)

    end

    上述程序段把数组中的每个元素都进行函数值计算,这样会耗费大量的运算时间,我们可以把数组看作一个整体来处理,计算函数值,可以修改这个程序段如下。

    n = 0:0.1:100;

    y = cos(n)

    通过使用MATLAB专门提供的测试程序运行时间的函数,可以发现,把数组看作一个整体,进行操作后,执行效率提高约300倍。

    另外,在必须使用多重循环的情况下,建议在循环的外环执行循环次数少的,内环执行循环次数多的,这样也可以显著提高程序执行速度。

    2.在使用数组或矩阵之前先定义维数

    MATLAB中的变量在使用之前不需要明确地定义和指定维数。但当未预定义数组或矩阵的维数时,当需赋值的元素下标超出现有的维数时,MATLAB 就为该数组或矩阵扩维一次,这样就会大大降低程序的执行效率。因此,在使用数组或矩阵之前,预定义维数可以提高程序的执行效率。

    3.对矩阵元素使用下标或者索引操作

    在MATLAB中,矩阵元素的引用可用两个下标来表示。例如:A(i,j)

    表示矩阵的第i行第j列的元素;A(1:k,j)表示矩阵A的第j列的前k个元素;A(:,j)

    表示矩阵的第j列的所有元素。求矩阵A的第j列元素的平均值的表达式为mean(A(:,j))。

    4.尽量多使用函数文件少使用脚本文件

    因为每次调用MATLAB的脚本文件都需要将不必要的中间变量加载到内存中,每执行一次,就加载一次。函数在调用时被编译成了伪代码,只需要加载到内存一次。当多次调用同一个函数时会运行快一些。因此尽量多使用函数文件而少使用脚本文件,也是提高执行效率的一种方法。

    5.在必须使用循环时,可以考虑转换为C-MEX

    当必须使用耗时的循环时,可以考虑将循环体中的语句转换为C-MEX。C-MEX是将M文件通过MATLAB的编译器转换为可执行文件,是按照

    MEX

    技术要求的格式编写相应的程序,通过编译连接,生成扩展名为.dll的动态链接库文件,可以在MATLAB环境下直接执行。这样,循环体中的语句在执行时不必每次都解释(interpret)。一般来说,C-MEX

    文件的执行速度是相同功能的M文件执行速率的20~40倍。编写C-MEX不同于M文件,需要了解MATLAB C-MEX规范。幸运的是MATLAB提供了将M文件转换为C-MEX的工具。

    6.内存优化

    MATLAB在进行复杂的运算时需要占用大量的内存。合理使用内存和提高内存的使用效率,可以加快运行速度,减少系统资源的占用。

    7.内存管理函数和命令

    (1)Clear variablename:从内存中删除名称为variablename的变量。

    (2)Clear all:从内存中删除所有的变量。

    (3)Save:将指令的变量存入磁盘。

    (4)Load:将save命令存入的变量载入内存。

    (5)Quit:退出MATLAB,并释放所有分配的内存。

    (6)Pack:把内存中的变量存入磁盘,再用内存中的连续空间载回这些变量。考虑到执行效率问题,不能在循环中使用。

    8.节约内存的方法

    (1)避免生成大的中间变量,并删除不再需要的临时变量。

    (2)当使用大的矩阵变量时,预先指定维数并分配好内存,避免每次临时扩充维数。

    (3)当程序需要生成大量变量数据时,可以考虑定期将变量写到磁盘,然后清除这些变量。

    当需要这些变量时,再重新从磁盘加载。

    (4)当矩阵中数据极少时,将全矩阵转换为稀疏矩阵。

    提高MATLAB程序效率的几点原则,这些都是俺在这两年中参加四次数模编写大量m程序总结的经验,加之网上很多英雄也是所见略同。

    1.“计算向量、矩阵化,尽量减少for循环。”[/B]

    因为MATLAB本来就是矩阵实验室的意思,他提供了极其强大而灵活的矩阵运算能力,你就没必要自己再用自己编写的for循环去实现矩阵运算的功能了。另外由于matlab是一种解释性语言,所以最忌讳直接使用循环语句。但在有些情况下,使用for循环可以提高程序的易读性,在效率提高不是很明显的情况下可以选择使用for循环。

    口说无凭,下面是利用tic与toc命令计算运算所用时间的方法,测试两种编程的效率。需要说明的是没有办法精确计算程序执行时间,matlab帮助这样写到“Keep

    in mind that tic and toc measure overall elapsed time. Make sure

    that no other applications are running in the background on your

    system that could affect the timing of your MATLAB

    programs.”意思是说在程序执行的背后很可能有其他程序在执行,这里涉及到程序进程的的问题,m程序执行的过程中很可能有其他进程中断m程序来利用cup,所以计算出来的时间就不仅是m程序的了,在这种情况下我把那些寄点去掉,进行多次计算求他的平均时间。

    n = 100;

    A(1:1000,1:1000) = 13;

    C(1:1000,1) = 15;

    D(1:1000,1) = 0;

    for k = 1:n

    D(:)

    = 0;

    tic

    for

    i = 1:1000

    for

    j = 1:1000

    D(i)

    = D(i) + A(i,j)*C(j);

    end

    end

    t1(k)

    = toc;

    %------------------

    D(:)

    = 0;

    tic

    D

    = A*C;

    t2(k)

    = toc;

    end

    u = t1./t2;

    u(u<0) = [];

    plot(u)

    p = mean(u)

    t1、t2分别代表用for循环编程和矩阵化编程计算矩阵乘向量所用时间,u代表时间的比值。u(u<0) =

    [];是认为t1不可能小于t2,所以去掉不可能出现的情况。然后画出图形计算平均值。

    经多次试验结果大致相同,其中一次结果如下:

    p =

    9.6196

    ------------t1时间是t2的十倍左右。

    2.“循环内大数组预先定义--预先分配空间”[/U]

    这一点原则是极其重要的,以至于在编写m程序时编辑器会给出提示“'ver' might be growing inside a

    loop.Consider prealloacting for speed.”

    clear

    n = 50;

    m = 1000;

    for k = 1:n

    A

    = [];

    tic

    A(1:m,1:m)

    = 3;

    for

    i = 1:m

    A(i,i)

    = i;

    end

    t1(k)

    = toc;

    %------------

    A

    = [];

    tic

    for

    j = 1:m

    A(j,j)

    = j;

    end

    t2(k)

    = toc;

    end

    t2(t1>10^9) = [];

    t1(t1>10^9) = [];

    plot([t1;t2]')

    t1、t2分别表示预先分配空间和循环中分配空间的时间,下图上面一条t2、下面t1

    3.“尽可能利用matlab内部提供的函数”[/U]

    因为matlab内部提供的函数绝对是各种问题的最优算法,那写程序都是他们大师级人物写出来的,程序应该说相当高效,有现成的为什么不用那! 这个原则就不用实际的程序测试了。

    关于MATLAB程序提速的问题,可以参考网上很多智者的文章,都比较经典。也可以看看我的上一篇文章,和网上大部分帖子有点不同,我是以实际的测试程序作为依据对如何提高MATLAB程序速度进行介绍的。 这里我再补充几点大家需要注意的。下面是我在国内一个比较出名的论坛看到的关于m程序提速的帖子,开始还真以为他们谈论的都应该遵循。(尽信书不如无书)

    帖子的一部分这样说道:“当要预分配一个非double型变量时使用repmat函数以加速,如将以下代码:

    A = int8(zeros(100));

    换成:

    A = repmat(int8(0), 100, 100);”

    凡事不能只凭自己的感觉,没有一点实际的例子,对于权威我们还要有挑战精神那,就不用说现在还不是经典的观点了。下面是我写的测试程序,我本来是想得到这位网友大哥的结果,但是实事不是我们想象的那么简单。

    n = 100;

    m = 1000;

    for k=1:n

    tic

    A

    = int8(ones(m));

    t1(k)

    = toc;

    tic

    B

    = repmat(int8(1),m,m);

    t2(k)

    = toc;

    end

    plot(1:n,t1,'r',1:n,t2)

    isequal(A,B)

    可以看出下面的红线是t1,而且最后的一句返回1,说明两种形式返回的值完全一样。

    由此我想说的是,不管是在我们做论文,还是写博客的时候,别直接从网上或者别人文章那儿找点知识定理之类的补充自己那苍白无力的文章。最好是自己动手编一下,“实践是检验真理的唯一标准”。

    经过这一测试,我感觉有必要,也有责任对这个论坛上的一大批经典谈论加以测试。尽管这个结论是错误的但这还不足以证明论坛上的帖子都不是经典。

    还有一点关于m程序提速的这样说到:“在必须使用多重循环时下,如果两个循环执行的次数不同,则在循环的外环执行循环次数少的,内环执行循环次数多的。这样可以显著提高速度。”

    n=1000;

    A = ones(1000)*13;

    for k=1:n

    tic

    for

    i=1:10

    for

    j=1:1000

    A(i,j)=A(i,j)*15;

    end

    end

    t1(k)=toc;

    tic

    for

    i=1:1000

    for

    j=1:10

    A(i,j)=A(i,j)*16;

    end

    end

    t2(k)=toc;

    end

    t2(t1>10^9)=[];

    t1(t1>10^9)=[];

    t1(t2>10^9)=[];

    t2(t2>10^9)=[];%去除外界因素影响所产生的寄点

    plot(1:size(t1,2),t1,'r',1:size(t1,2),t2)

    由这个时间图可以看出for循环的嵌套顺序对于速度是有影响的,虽然相对来说差别不是很大,但是毕竟论坛上的观点是正确的。至于他所说的“显著”二字就没必要加上了。

    此论坛还有一些提速的观点列举如下:

    “遵守Performance Acceleration的规则

    关于什么是“Performance Acceleration”请参阅matlab的帮助文件。我只简要的将

    其规则总结如下7条:

    1、只有使用以下数据类型,matlab才会对其加速:

    logical,char,int8,uint8,int16,uint16,int32,uint32,double

    而语句中如果使用了非以上的数据类型则不会加速,如:numeric,cell,structure,single,function

    handle,java classes,user classes,int64,uint64

    2、matlab不会对超过三维的数组进行加速。

    3、当使用for循环时,只有遵守以下规则才会被加速:a、for循环的范围只用标量值来表示;

    b、for循环内部的每一条语句都要满足上面的两条规则,即只使用支持加速的数据类型,只使用三维以下的数组;c、循环内只调用了内建函数(build-in

    function)。

    4、当使用if、elseif、while和switch时,其条件测试语句中只使用了标量值时,将加速运行。

    5、不要在一行中写入多条操作,这样会减慢运行速度。即不要有这样的语句:

    x = a.name; for k=1:10000, sin(A(k)), end;

    6、当某条操作改变了原来变量的数据类型或形状(大小,维数)时将会减慢运行速度。

    7、应该这样使用复常量x = 7 + 2i,而不应该这样使用:x = 7 + 2*i,后者会降低运行速度。”

    “尽量用向量化的运算来代替循环操作。如将下面的程序:

    i=0;

    for t = 0:.01:10

    i

    = i+1;

    y(i)

    = sin(t);

    end

    替换为:

    t = 0:.01:10;

    y = sin(t);

    速度将会大大加快。最常用的使用vectorizing技术的函数有:All、diff、ipermute、permute、reshape、ueeze、y、find、logical、prod、shiftdim、sub2ind、cumsum、ind2sub、ndgrid、repmat、sort、sum

    等。”

    “优先使用matlab内建函数,将耗时的循环编写进MEX-File中以获得加速。

    b、使用Functions而不是Scripts 。”

    “ 绝招:你也许觉得下面两条是屁话,但有时候它真的是解决问题的最好方法。

    1、改用更有效的算法

    2、采用Mex技术,或者利用matlab提供的工具将程序转化为C语言、Fortran语言。关于如何将M文件转化为C语言程序运行,可以参阅本版帖子:“总结:m文件转化为c/c++语言文件,VC编译”。

    除了m程序提速的问题,这里还列出了《MATLAB代码矢量化指南(译)》

    一、基本技术 ----------------------------------------------------- 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: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处理稀疏矩阵的方式。

    关于MATLAB的效率问题,很多文章,包括我之前写的一些,主要集中在使用向量化以及相关的问题上。但是,最近我在实验时对代码进行profile的过程中,发现在新版本的MATLAB下,for-loop已经得到了极大优化,而效率的瓶颈更多是在函数调用和索引访问的过程中。

    由于MATLAB特有的解释过程,不同方式的函数调用和元素索引,其效率差别巨大。不恰当的使用方式可能在本来不起眼的地方带来严重的开销,甚至可能使你的代码的运行时间增加上千倍(这就不是多买几台服务器或者增加计算节点能解决的了,呵呵)。

    下面通过一些简单例子说明问题。(实验选在装有Windows Vista的一台普通的台式机(Core2 Duo 2.33GHz +

    4GB Ram)进行,相比于计算集群,

    这可能和大部分朋友的环境更相似一些。实验过程是对某一个过程实施多次的整体进行计时,然后得到每次过程的平均时间,以减少计时误差带来的影响。多次实验,在均值附近正负20%的范围内的置信度高于95%。为了避免算上首次运行时花在预编译上的时间,在开始计时前都进行充分的“热身”运行。) 函数调用的效率

    一个非常简单的例子,把向量中每个元素加1。(当然这个例子根本不需要调函数,但是,用它主要是为了减少函数执行本身的时间,突出函数解析和调用的过程。)

    作为baseline,先看看最直接的实现

    % input u: u is a 1000000 x 1 vectorv = u +

    1; 这个过程平均需要0.00105 sec。

    而使用长期被要求尽量避免的for-loop n = numel(u);% v = zeros(n, 1) has been pre-allocated.for i = 1 :

    n v(i)

    = u(i) + 1;end

    所需的平均时间大概是0.00110

    sec。从统计意义上说,和vectorization已经没有显著差别。无论是for-loop或者vectorization,每秒平均进行约10亿次“索引-读取-加法-写入”的过程,计算资源应该得到了比较充分的利用。

    要是这个过程使用了函数调用呢?MATLAB里面支持很多种函数调用方式,主要的有m-function, function handle,

    anonymous function, inline, 和feval,而feval的主参数可以是字符串名字,function

    handle, anonymous function或者inline。

    用m-function,就是专门定义一个函数

    function y =

    fm(x) y

    = x + 1;

    在调用时 for i = 1 :

    n v(i)

    = fm(u(i));end

    function

    handle就是用@来把一个function赋到一个变量上,类似于C/C++的函数指针,或者C#里面的delegate的作用

    fh = @fm;for i = 1 :

    n v(i)

    = fh(u(i));end

    anonymous function是一种便捷的语法来构造简单的函数,类似于LISP, Python的lambda表达式

    fa = @(x) x + 1;for i = 1 :

    n v(i)

    = fa(u(i));end

    inline function是一种传统的通过表达式字符串构造函数的过程

    fi = inline('x + 1', 'x');for i = 1 :

    n v(i)

    = fi(u(i));end

    feval的好处在于可以以字符串方式指定名字来调用函数,当然它也可以接受别的参数。

    v(i) = feval_r('fm', u(i));v(i) = feval_r(fh, u(i));v(i) = feval_r(fa,

    u(i)); 对于100万次调用(包含for-loop本身的开销,函数解析(resolution),压栈,执行加法,退栈,把返回值赋给接收变量),不同的方式的时间差别很大:

    m-function0.385 secfunction handle0.615 secanonymous function0.635

    secinline function166.00 secfeval_r('fm', u(i))8.328 secfeval_r(fh,

    u(i))0.618 secfeval_r(fa, u(i))0.652 secfeval_r(@fm, u(i))2.788

    secfeval_r(@fa, u(i))34.689 sec

    从这里面,我们可以看到几个有意思的现象:

    · 首先,调用自定义函数的开销远远高于一个简单的计算过程。直接写

    u(i) = v(i) + 1 只需要 0.0011

    秒左右,而写u(i) = fm(v(i))

    则需要0.385秒,时间多了几百倍,而它们干的是同一件事情。这说明了,函数调用的开销远远大于for-loop自己的开销和简单计算过程——在不同情况可能有点差别,一般而言,一次普通的函数调用花费的时间相当于进行了几百次到一两千次双精度浮点加法。

    · 使用function handle和anonymous

    function的开销比使用普通的m-函数要高一些。这可能是由于涉及到句柄解析的时间,而普通的函数在第一次运行前已经在内存中进行预编译。

    · inline

    function的运行时间则令人难以接受了,竟然需要一百多秒(是普通函数调用的四百多倍,是直接计算的十几万倍)。这是因为matlab是在每次运行时临时对字符串表达式(或者它的某种不太高效的内部表达)进行parse。

    · feval_r(fh,

    u(i))和fh(u(i)),feval_r(fa,

    u(i))和fa(u(i))的运行时间很接近,表明feval在接受句柄为主参数时本身开销很小。但是,surprising的是 for

    i = 1 :

    n v(i)

    = feval_r(@fm, u(i));end比起fh = @fm;for i = 1 :

    n v(i)

    = feval_r(fh, u(i));end慢了4.5倍 (前者0.618秒,后者2.788秒)。

    for i = 1 :

    n v(i)

    = feval_r(@(x) x + 1, u(i));end比起fa = @(x) x + 1;for i = 1 :

    n v(i)

    = feval_r(fa, u(i));end竟然慢了53倍(前者0.652秒,后者34.689秒)。

    由于在MATLAB的内部实现中,function

    handle的解析是在赋值过程中进行的,所以预先用一个变量把句柄接下,其效果就是预先完成了句柄解析,而如果直接把@fm或者@(x) x

    +

    1写在参数列上,虽然写法简洁一些,但是解析过程是把参数被赋值到所调函数的局部变量时才进行,每调用一次解析一次,造成了巨大的开销。

    · feval使用字符串作为函数名字时,所耗时间比传入句柄大,因为这涉及到对函数进行搜索的时间(当然这个搜索是在一个索引好的cache里面进行(除了第一次),而不是在所有path指定的路径中搜索。)

    在2007年以后,MATLAB推出了arrayfun函数,上面的for-loop可以写为 v =

    arrayfun(fh, u)

    这平均需要4.48 sec,这比起for-loop(需时0.615

    sec)还慢了7倍多。这个看上去“消除了for-loop"的函数,由于其内部设计的原因,未必能带来效率上的正效果。 元素和域的访问

    除了函数调用,数据的访问方式对于效率也有很大影响。MATLAB主要支持下面一些形式的访问:

    · array-index A(i):

    · cell-index: C{i};

    · struct

    field: S.fieldname

    · struct field

    (dynamic): S.('fieldname')这里主要探索单个元素或者域的访问(当然,MATLAB也支持对于子数组的非常灵活整体索引)。

    对于一百万次访问的平均时间

    A(i) for a numeric array0.0052 secC{i} for a cell array0.2568

    secstruct field0.0045 secstruct field (with dynamic name)1.0394

    sec

    我们可以看到MATLAB对于单个数组元素或者静态的struct

    field的访问,可以达到不错的速度,在主流台式机约每秒2亿次(连同for-loop的开销)。而cell

    array的访问则明显缓慢,约每秒400万次(慢了50倍)。MATLAB还支持灵活的使用字符串来指定要访问域的语法(动态名字),但是,是以巨大的开销为代价的,比起静态的访问慢了200倍以上。

    关于Object-oriented Programming

    MATLAB在新的版本中(尤其是2008版),对于面向对象的编程提供了强大的支持。在2008a中,它对于OO的支持已经不亚于python等的高级脚本语言。但是,我在实验中看到,虽然在语法上提供了全面的支持,但是matlab里面面向对象的效率很低,开销巨大。这里仅举几个例子。

    ·object中的property的访问速度是3500万次,比struct

    field慢了6-8倍。MATLAB提供了一种叫做dependent

    property的属性,非常灵活,但是,效率更低,平均每秒访问速度竟然低至2.6万次(这种速度基本甚至难以用于中小规模的应用中)。

    ·object中method调用的效率也明显低于普通函数的调用,对于instance

    method,每百万次调用,平均需时5.8秒,而对于static

    method,每百万次调用需时25.8秒。这相当于每次调用都需要临时解析的速度,而matlab的类方法解析的效率目前还明显偏低。

    ·MATLAB中可以通过改写subsref和subsasgn的方法,对于对象的索引和域的访问进行非常灵活的改造,可以通过它创建类似于数组的对象。但是,一个符合要求的subsref的行为比较复杂。在一个提供了subsref的对象中,大部分行为都需要subsref进行调度,而默认的较优的调度方式将被掩盖。在一个提供了subsref的类中(使用一种最快捷的实现),object

    property的平均访问速度竟然降到1万5千次每秒。建议[/U]

    根据上面的分析,对于撰写高效MATLAB代码,我有下面一些建议:

    1. 虽然for-loop的速度有了很大改善,vectorization(向量化)仍旧是改善效率的重要途径,尤其是在能把运算改写成矩阵乘法的情况下,改善尤为显著。

    2. 在不少情况下,for-loop本身已经不构成太大问题,尤其是当循环体本身需要较多的计算的时候。这个时候,改善概率的关键在于改善循环体本身而不是去掉for-loop。

    3. MATLAB的函数调用过程(非built-in

    function)有显著开销,因此,在效率要求较高的代码中,应该尽可能采用扁平的调用结构,也就是在保持代码清晰和可维护的情况下,尽量直接写表达式和利用built-in

    function,避免不必要的自定义函数调用过程。在次数很多的循环体内(包括在cellfun,

    arrayfun等实际上蕴含循环的函数)形成长调用链,会带来很大的开销。

    4. 在调用函数时,首选built-in function,然后是普通的m-file函数,然后才是function

    handle或者anonymous function。在使用function handle或者anonymous

    function作为参数传递时,如果该函数被调用多次,最好先用一个变量接住,再传入该变量。这样,可以有效避免重复的解析过程。

    5. 在可能的情况下,使用numeric array或者struct array,它们的效率大幅度高于cell

    array(几十倍甚至更多)。对于struct,尽可能使用普通的域(字段,field)访问方式,在非效率关键,执行次数较少,而灵活性要求较高的代码中,可以考虑使用动态名称的域访问。

    6. 虽然object-oriented从软件工程的角度更为优胜,而且object的使用很多时候很方便,但是MATLAB目前对于OO的实现效率很低,在效率关键的代码中应该慎用objects。

    7. 如果需要设计类,应该尽可能采用普通的property,而避免灵活但是效率很低的dependent

    property。如非确实必要,避免重载subsref和subsasgn函数,因为这会全面接管对于object的接口调用,往往会带来非常巨大的开销(成千上万倍的减慢),甚至使得本来几乎不是问题的代码成为性能瓶颈。

    展开全文
  • MATLAB数值计算

    2016-10-07 20:15:16
    MATLAB数值计算》是关于数值方法、MATLAB软件和工程计算的教材,着重介绍数学软件的熟练使用及其内在的高效率算法。主要内容包括:MATLAB介绍、线性方程组、插值、方程求根、最小二乘法、数值积分、常微分方程、...
  • 用过Matlab的人都知道,Matlab是一种解释性语言,存在计算速度慢的问题,为了提高程序的运行效率matlab提供了多种实用工具及编码技巧。1. 循环矢量化Matlab是为矢量和矩阵操作而设计的,因此,可以通过矢量化方法...

    用过Matlab的人都知道,Matlab是一种解释性语言,存在计算速度慢的问题,为了提高程序的运行效率,matlab提供了多种实用工具及编码技巧。

    1. 循环矢量化

    Matlab是为矢量和矩阵操作而设计的,因此,可以通过矢量化方法加速M文件的运行。矢量化是指将for循环和while循环转换为等价的矢量或矩阵操作。下面给出一个循环的例子:

    i=0;

    for n = 0:0.1:1000

    i=i+1;

    y(i)=cos(n);

    end

    那么我们可以矢量化为:

    n= 0:0.1:1000;

    y=cos(n);

    我们可以用tic和toc函数来查看上述各代码运行的时间,采用for循环的程序0.39秒(具体时间和计算机配置有关),而矢量化后几乎耗时为0。

    2. 给数组或矩阵预分配内存

    特别是使用大型数组或矩阵时,Matlab进行动态内存分配和取消时,可能会产生内存碎片,这将导致大量闲置内存产生,预分配可通过提前给大型数据结构预约足够空间来避免这个问题。

    3. 用函数代替脚本文件

    因为每次调用MATLAB的脚本文件都需要将不必要的中间变量加载到内存中,每执行一次,就加载一次。函数在调用时被编译成了伪代码,只需要加载到内存一次。当多次调用同一个函数时会运行快一些。因此尽量多使用函数文件而少使用脚本文件,也是提高执行效率的一种方法。

    4. 用Mex文件编写循环代码

    Matlab提供了与C和C++的接口,那么我们可以在用C或C++语言编写耗时的循环代码,然后通过接口程序在Matlab中转换成dll文件,这就是我们所要的Mex文件,通过这种方法可以极大地提高计算速率。

    1.尽量避免使用循环结构

    MATLAB变量的基本类型是矩阵,当对矩阵的每个元素循环处理时,运算速度很慢。因此编程时应尽量把数组和矩阵看作一个整体来进行编程,而不是像其他的程序设计语言那样,使用循环结构对矩阵的元素循环进行处理。利用MATLAB提供的用于矢量化操作的函数,把循环矢量化,这样既可以提高编程效率,也可以提高程序的执行效率。下面给出一个循环的例子:

    i=0;

    for n = 0:0.1:100

    i=i+1;

    y(i)=cos(n)

    end

    上述程序段把数组中的每个元素都进行函数值计算,这样会耗费大量的运算时间,我们可以把数组看作一个整体来处理,计算函数值,可以修改这个程序段如下。

    n = 0:0.1:100;

    y = cos(n)

    通过使用MATLAB专门提供的测试程序运行时间的函数,可以发现,把数组看作一个整体,进行操作后,执行效率提高约300倍。

    另外,在必须使用多重循环的情况下,建议在循环的外环执行循环次数少的,内环执行循环次数多的,这样也可以显著提高程序执行速度。

    2.在使用数组或矩阵之前先定义维数

    MATLAB中的变量在使用之前不需要明确地定义和指定维数。但当未预定义数组或矩阵的维数时,当需赋值的元素下标超出现有的维数时,MATLAB 就为该数组或矩阵扩维一次,这样就会大大降低程序的执行效率。因此,在使用数组或矩阵之前,预定义维数可以提高程序的执行效率。

    3.对矩阵元素使用下标或者索引操作

    在MATLAB中,矩阵元素的引用可用两个下标来表示。例如:A(i,j)

    表示矩阵的第i行第j列的元素;A(1:k,j)表示矩阵A的第j列的前k个元素;A(:,j)

    表示矩阵的第j列的所有元素。求矩阵A的第j列元素的平均值的表达式为mean(A(:,j))。

    4.尽量多使用函数文件少使用脚本文件

    因为每次调用MATLAB的脚本文件都需要将不必要的中间变量加载到内存中,每执行一次,就加载一次。函数在调用时被编译成了伪代码,只需要加载到内存一次。当多次调用同一个函数时会运行快一些。因此尽量多使用函数文件而少使用脚本文件,也是提高执行效率的一种方法。

    5.在必须使用循环时,可以考虑转换为C-MEX

    当必须使用耗时的循环时,可以考虑将循环体中的语句转换为C-MEX。C-MEX是将M文件通过MATLAB的编译器转换为可执行文件,是按照

    MEX

    技术要求的格式编写相应的程序,通过编译连接,生成扩展名为.dll的动态链接库文件,可以在MATLAB环境下直接执行。这样,循环体中的语句在执行时不必每次都解释(interpret)。一般来说,C-MEX

    文件的执行速度是相同功能的M文件执行速率的20~40倍。编写C-MEX不同于M文件,需要了解MATLAB C-MEX规范。幸运的是MATLAB提供了将M文件转换为C-MEX的工具。

    6.内存优化

    MATLAB在进行复杂的运算时需要占用大量的内存。合理使用内存和提高内存的使用效率,可以加快运行速度,减少系统资源的占用。

    7.内存管理函数和命令

    (1)Clear

    variablename:从内存中删除名称为variablename的变量。

    (2)Clear

    all:从内存中删除所有的变量。

    (3)Save:将指令的变量存入磁盘。

    (4)Load:将save命令存入的变量载入内存。

    (5)Quit:退出MATLAB,并释放所有分配的内存。

    (6)Pack:把内存中的变量存入磁盘,再用内存中的连续空间载回这些变量。考虑到执行效率问题,不能在循环中使用。

    8.节约内存的方法

    (1)避免生成大的中间变量,并删除不再需要的临时变量。

    (2)当使用大的矩阵变量时,预先指定维数并分配好内存,避免每次临时扩充维数。

    (3)当程序需要生成大量变量数据时,可以考虑定期将变量写到磁盘,然后清除这些变量。

    当需要这些变量时,再重新从磁盘加载。

    (4)当矩阵中数据极少时,将全矩阵转换为稀疏矩阵。

    提高MATLAB程序效率的几点原则,这些都是俺在这两年中参加四次数模编写大量m程序总结的经验,加之网上很多英雄也是所见略同。

    1.“计算向量、矩阵化,尽量减少for循环。”[/B]

    因为MATLAB本来就是矩阵实验室的意思,他提供了极其强大而灵活的矩阵运算能力,你就没必要自己再用自己编写的for循环去实现矩阵运算的功能了。另外由于matlab是一种解释性语言,所以最忌讳直接使用循环语句。但在有些情况下,使用for循环可以提高程序的易读性,在效率提高不是很明显的情况下可以选择使用for循环。

    口说无凭,下面是利用tic与toc命令计算运算所用时间的方法,测试两种编程的效率。需要说明的是没有办法精确计算程序执行时间,matlab帮助这样写到“Keep

    in mind that tic and toc measure overall elapsed time. Make sure

    that no other applications are running in the background on your

    system that could affect the timing of your MATLAB

    programs.”意思是说在程序执行的背后很可能有其他程序在执行,这里涉及到程序进程的的问题,m程序执行的过程中很可能有其他进程中断m程序来利用cup,所以计算出来的时间就不仅是m程序的了,在这种情况下我把那些寄点去掉,进行多次计算求他的平均时间。

    n = 100;

    A(1:1000,1:1000) = 13;

    C(1:1000,1) = 15;

    D(1:1000,1) = 0;

    for k = 1:n

    D(:)

    = 0;

    tic

    for

    i = 1:1000

    for

    j = 1:1000

    D(i)

    = D(i) + A(i,j)*C(j);

    end

    end

    t1(k)

    = toc;

    %------------------

    D(:)

    = 0;

    tic

    D

    = A*C;

    t2(k)

    = toc;

    end

    u = t1./t2;

    u(u<0) = [];

    plot(u)

    p = mean(u)

    t1、t2分别代表用for循环编程和矩阵化编程计算矩阵乘向量所用时间,u代表时间的比值。u(u<0)

    = [];是认为t1不可能小于t2,所以去掉不可能出现的情况。然后画出图形计算平均值。

    经多次试验结果大致相同,其中一次结果如下:

    p =

    9.6196

    ------------t1时间是t2的十倍左右。

    2.“循环内大数组预先定义--预先分配空间”[/U]

    这一点原则是极其重要的,以至于在编写m程序时编辑器会给出提示“'ver'

    might be growing inside a loop.Consider prealloacting for

    speed.”

    clear

    n = 50;

    m = 1000;

    for k = 1:n

    A

    = [];

    tic

    A(1:m,1:m)

    = 3;

    for

    i = 1:m

    A(i,i)

    = i;

    end

    t1(k)

    = toc;

    %------------

    A

    = [];

    tic

    for

    j = 1:m

    A(j,j)

    = j;

    end

    t2(k)

    = toc;

    end

    t2(t1>10^9) = [];

    t1(t1>10^9) = [];

    plot([t1;t2]')

    t1、t2分别表示预先分配空间和循环中分配空间的时间,下图上面一条t2、下面t1

    3.“尽可能利用matlab内部提供的函数”[/U]

    因为matlab内部提供的函数绝对是各种问题的最优算法,那写程序都是他们大师级人物写出来的,程序应该说相当高效,有现成的为什么不用那! 这个原则就不用实际的程序测试了。

    关于MATLAB程序提速的问题,可以参考网上很多智者的文章,都比较经典。也可以看看我的上一篇文章,和网上大部分帖子有点不同,我是以实际的测试程序作为依据对如何提高MATLAB程序速度进行介绍的。 这里我再补充几点大家需要注意的。下面是我在国内一个比较出名的论坛看到的关于m程序提速的帖子,开始还真以为他们谈论的都应该遵循。(尽信书不如无书)

    帖子的一部分这样说道:“当要预分配一个非double型变量时使用repmat函数以加速,如将以下代码:

    A = int8(zeros(100));

    换成:

    A = repmat(int8(0), 100,

    100);”

    凡事不能只凭自己的感觉,没有一点实际的例子,对于权威我们还要有挑战精神那,就不用说现在还不是经典的观点了。下面是我写的测试程序,我本来是想得到这位网友大哥的结果,但是实事不是我们想象的那么简单。

    n = 100;

    m = 1000;

    for k=1:n

    tic

    A

    = int8(ones(m));

    t1(k)

    = toc;

    tic

    B

    = repmat(int8(1),m,m);

    t2(k)

    = toc;

    end

    plot(1:n,t1,'r',1:n,t2)

    isequal(A,B)

    可以看出下面的红线是t1,而且最后的一句返回1,说明两种形式返回的值完全一样。

    由此我想说的是,不管是在我们做论文,还是写博客的时候,别直接从网上或者别人文章那儿找点知识定理之类的补充自己那苍白无力的文章。最好是自己动手编一下,“实践是检验真理的唯一标准”。

    经过这一测试,我感觉有必要,也有责任对这个论坛上的一大批经典谈论加以测试。尽管这个结论是错误的但这还不足以证明论坛上的帖子都不是经典。

    还有一点关于m程序提速的这样说到:“在必须使用多重循环时下,如果两个循环执行的次数不同,则在循环的外环执行循环次数少的,内环执行循环次数多的。这样可以显著提高速度。”

    n=1000;

    A = ones(1000)*13;

    for k=1:n

    tic

    for

    i=1:10

    for

    j=1:1000

    A(i,j)=A(i,j)*15;

    end

    end

    t1(k)=toc;

    tic

    for

    i=1:1000

    for

    j=1:10

    A(i,j)=A(i,j)*16;

    end

    end

    t2(k)=toc;

    end

    t2(t1>10^9)=[];

    t1(t1>10^9)=[];

    t1(t2>10^9)=[];

    t2(t2>10^9)=[];%去除外界因素影响所产生的寄点

    plot(1:size(t1,2),t1,'r',1:size(t1,2),t2)

    由这个时间图可以看出for循环的嵌套顺序对于速度是有影响的,虽然相对来说差别不是很大,但是毕竟论坛上的观点是正确的。至于他所说的“显著”二字就没必要加上了。

    此论坛还有一些提速的观点列举如下:

    “遵守Performance Acceleration的规则

    关于什么是“Performance

    Acceleration”请参阅matlab的帮助文件。我只简要的将

    其规则总结如下7条:

    1、只有使用以下数据类型,matlab才会对其加速:

    logical,char,int8,uint8,int16,uint16,int32,uint32,double

    而语句中如果使用了非以上的数据类型则不会加速,如:numeric,cell,structure,single,function

    handle,java classes,user classes,int64,uint64

    2、matlab不会对超过三维的数组进行加速。

    3、当使用for循环时,只有遵守以下规则才会被加速:a、for循环的范围只用标量值来表示;

    b、for循环内部的每一条语句都要满足上面的两条规则,即只使用支持加速的数据类型,只使用三维以下的数组;c、循环内只调用了内建函数(build-in

    function)。

    4、当使用if、elseif、while和switch时,其条件测试语句中只使用了标量值时,将加速运行。

    5、不要在一行中写入多条操作,这样会减慢运行速度。即不要有这样的语句:

    x = a.name;

    for k=1:10000, sin(A(k)), end;

    6、当某条操作改变了原来变量的数据类型或形状(大小,维数)时将会减慢运行速度。

    7、应该这样使用复常量x = 7 + 2i,而不应该这样使用:x = 7 +

    2*i,后者会降低运行速度。”

    “尽量用向量化的运算来代替循环操作。如将下面的程序:

    i=0;

    for t = 0:.01:10

    i

    = i+1;

    y(i)

    = sin(t);

    end

    替换为:

    t = 0:.01:10;

    y = sin(t);

    速度将会大大加快。最常用的使用vectorizing技术的函数有:All、diff、ipermute、permute、reshape、ueeze、y、find、logical、prod、shiftdim、sub2ind、cumsum、ind2sub、ndgrid、repmat、sort、sum

    等。”

    “优先使用matlab内建函数,将耗时的循环编写进MEX-File中以获得加速。

    b、使用Functions而不是Scripts 。”

    “ 绝招:你也许觉得下面两条是屁话,但有时候它真的是解决问题的最好方法。

    1、改用更有效的算法

    2、采用Mex技术,或者利用matlab提供的工具将程序转化为C语言、Fortran语言。关于如何将M文件转化为C语言程序运行,可以参阅本版帖子:“总结:m文件转化为c/c++语言文件,VC编译”。

    除了m程序提速的问题,这里还列出了《MATLAB代码矢量化指南(译)》

    一、基本技术

    -----------------------------------------------------

    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: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处理稀疏矩阵的方式。

    关于MATLAB的效率问题,很多文章,包括我之前写的一些,主要集中在使用向量化以及相关的问题上。但是,最近我在实验时对代码进行profile的过程中,发现在新版本的MATLAB下,for-loop已经得到了极大优化,而效率的瓶颈更多是在函数调用和索引访问的过程中。

    由于MATLAB特有的解释过程,不同方式的函数调用和元素索引,其效率差别巨大。不恰当的使用方式可能在本来不起眼的地方带来严重的开销,甚至可能使你的代码的运行时间增加上千倍(这就不是多买几台服务器或者增加计算节点能解决的了,呵呵)。

    下面通过一些简单例子说明问题。(实验选在装有Windows

    Vista的一台普通的台式机(Core2 Duo 2.33GHz + 4GB Ram)进行,相比于计算集群,

    这可能和大部分朋友的环境更相似一些。实验过程是对某一个过程实施多次的整体进行计时,然后得到每次过程的平均时间,以减少计时误差带来的影响。多次实验,在均值附近正负20%的范围内的置信度高于95%。为了避免算上首次运行时花在预编译上的时间,在开始计时前都进行充分的“热身”运行。)

    函数调用的效率

    一个非常简单的例子,把向量中每个元素加1。(当然这个例子根本不需要调函数,但是,用它主要是为了减少函数执行本身的时间,突出函数解析和调用的过程。)

    作为baseline,先看看最直接的实现

    % input u: u is a 1000000 x 1

    vectorv = u + 1;

    这个过程平均需要0.00105 sec。

    而使用长期被要求尽量避免的for-loop

    n = numel(u);% v = zeros(n, 1) has

    been pre-allocated.for i = 1 :

    n v(i)

    = u(i) + 1;end

    所需的平均时间大概是0.00110

    sec。从统计意义上说,和vectorization已经没有显著差别。无论是for-loop或者vectorization,每秒平均进行约10亿次“索引-读取-加法-写入”的过程,计算资源应该得到了比较充分的利用。

    要是这个过程使用了函数调用呢?MATLAB里面支持很多种函数调用方式,主要的有m-function,

    function handle, anonymous function, inline,

    和feval,而feval的主参数可以是字符串名字,function handle, anonymous

    function或者inline。

    用m-function,就是专门定义一个函数

    function y =

    fm(x) y

    = x + 1;

    在调用时

    for i = 1 :

    n v(i)

    = fm(u(i));end

    function

    handle就是用@来把一个function赋到一个变量上,类似于C/C++的函数指针,或者C#里面的delegate的作用

    fh = @fm;for i = 1 :

    n v(i)

    = fh(u(i));end

    anonymous

    function是一种便捷的语法来构造简单的函数,类似于LISP, Python的lambda表达式

    fa = @(x) x + 1;for i = 1 :

    n v(i)

    = fa(u(i));end

    inline

    function是一种传统的通过表达式字符串构造函数的过程

    fi = inline('x + 1', 'x');for i = 1

    :

    n v(i)

    = fi(u(i));end

    feval的好处在于可以以字符串方式指定名字来调用函数,当然它也可以接受别的参数。

    v(i) = feval_r('fm', u(i));v(i) =

    feval_r(fh, u(i));v(i) = feval_r(fa,

    u(i));

    对于100万次调用(包含for-loop本身的开销,函数解析(resolution),压栈,执行加法,退栈,把返回值赋给接收变量),不同的方式的时间差别很大:

    m-function0.385 secfunction handle0.615 secanonymous function0.635

    secinline function166.00 secfeval_r('fm', u(i))8.328 secfeval_r(fh,

    u(i))0.618 secfeval_r(fa, u(i))0.652 secfeval_r(@fm, u(i))2.788

    secfeval_r(@fa, u(i))34.689 sec

    从这里面,我们可以看到几个有意思的现象:

    首先,调用自定义函数的开销远远高于一个简单的计算过程。直接写 u(i)

    = v(i) + 1 只需要 0.0011 秒左右,而写u(i)

    = fm(v(i))

    则需要0.385秒,时间多了几百倍,而它们干的是同一件事情。这说明了,函数调用的开销远远大于for-loop自己的开销和简单计算过程——在不同情况可能有点差别,一般而言,一次普通的函数调用花费的时间相当于进行了几百次到一两千次双精度浮点加法。

    使用function handle和anonymous

    function的开销比使用普通的m-函数要高一些。这可能是由于涉及到句柄解析的时间,而普通的函数在第一次运行前已经在内存中进行预编译。

    inline

    function的运行时间则令人难以接受了,竟然需要一百多秒(是普通函数调用的四百多倍,是直接计算的十几万倍)。这是因为matlab是在每次运行时临时对字符串表达式(或者它的某种不太高效的内部表达)进行parse。

    feval_r(fh, u(i))和fh(u(i)),feval_r(fa,

    u(i))和fa(u(i))的运行时间很接近,表明feval在接受句柄为主参数时本身开销很小。但是,surprising的是 for

    i = 1 :

    n v(i)

    = feval_r(@fm, u(i));end比起fh = @fm;for i = 1 :

    n v(i)

    = feval_r(fh, u(i));end慢了4.5倍 (前者0.618秒,后者2.788秒)。

    for i = 1 :

    n v(i)

    = feval_r(@(x) x + 1, u(i));end比起fa = @(x) x + 1;for i = 1 :

    n v(i)

    = feval_r(fa, u(i));end竟然慢了53倍(前者0.652秒,后者34.689秒)。

    由于在MATLAB的内部实现中,function

    handle的解析是在赋值过程中进行的,所以预先用一个变量把句柄接下,其效果就是预先完成了句柄解析,而如果直接把@fm或者@(x) x

    +

    1写在参数列上,虽然写法简洁一些,但是解析过程是把参数被赋值到所调函数的局部变量时才进行,每调用一次解析一次,造成了巨大的开销。

    feval使用字符串作为函数名字时,所耗时间比传入句柄大,因为这涉及到对函数进行搜索的时间(当然这个搜索是在一个索引好的cache里面进行(除了第一次),而不是在所有path指定的路径中搜索。)

    在2007年以后,MATLAB推出了arrayfun函数,上面的for-loop可以写为 v =

    arrayfun(fh, u)

    这平均需要4.48 sec,这比起for-loop(需时0.615

    sec)还慢了7倍多。这个看上去“消除了for-loop"的函数,由于其内部设计的原因,未必能带来效率上的正效果。

    元素和域的访问

    除了函数调用,数据的访问方式对于效率也有很大影响。MATLAB主要支持下面一些形式的访问:

    array-index A(i):

    cell-index: C{i};

    struct field: S.fieldname

    struct field

    (dynamic): S.('fieldname')

    这里主要探索单个元素或者域的访问(当然,MATLAB也支持对于子数组的非常灵活整体索引)。

    对于一百万次访问的平均时间

    A(i) for a numeric array0.0052 secC{i} for a cell array0.2568

    secstruct field0.0045 secstruct field (with dynamic name)1.0394

    sec

    我们可以看到MATLAB对于单个数组元素或者静态的struct

    field的访问,可以达到不错的速度,在主流台式机约每秒2亿次(连同for-loop的开销)。而cell

    array的访问则明显缓慢,约每秒400万次(慢了50倍)。MATLAB还支持灵活的使用字符串来指定要访问域的语法(动态名字),但是,是以巨大的开销为代价的,比起静态的访问慢了200倍以上。

    关于Object-oriented Programming

    MATLAB在新的版本中(尤其是2008版),对于面向对象的编程提供了强大的支持。在2008a中,它对于OO的支持已经不亚于python等的高级脚本语言。但是,我在实验中看到,虽然在语法上提供了全面的支持,但是matlab里面面向对象的效率很低,开销巨大。这里仅举几个例子。

    object中的property的访问速度是3500万次,比struct

    field慢了6-8倍。MATLAB提供了一种叫做dependent

    property的属性,非常灵活,但是,效率更低,平均每秒访问速度竟然低至2.6万次(这种速度基本甚至难以用于中小规模的应用中)。

    object中method调用的效率也明显低于普通函数的调用,对于instance

    method,每百万次调用,平均需时5.8秒,而对于static

    method,每百万次调用需时25.8秒。这相当于每次调用都需要临时解析的速度,而matlab的类方法解析的效率目前还明显偏低。

    MATLAB中可以通过改写subsref和subsasgn的方法,对于对象的索引和域的访问进行非常灵活的改造,可以通过它创建类似于数组的对象。但是,一个符合要求的subsref的行为比较复杂。在一个提供了subsref的对象中,大部分行为都需要subsref进行调度,而默认的较优的调度方式将被掩盖。在一个提供了subsref的类中(使用一种最快捷的实现),object

    property的平均访问速度竟然降到1万5千次每秒。

    建议[/U]

    根据上面的分析,对于撰写高效MATLAB代码,我有下面一些建议:

    虽然for-loop的速度有了很大改善,vectorization(向量化)仍旧是改善效率的重要途径,尤其是在能把运算改写成矩阵乘法的情况下,改善尤为显著。

    在不少情况下,for-loop本身已经不构成太大问题,尤其是当循环体本身需要较多的计算的时候。这个时候,改善概率的关键在于改善循环体本身而不是去掉for-loop。

    MATLAB的函数调用过程(非built-in

    function)有显著开销,因此,在效率要求较高的代码中,应该尽可能采用扁平的调用结构,也就是在保持代码清晰和可维护的情况下,尽量直接写表达式和利用built-in

    function,避免不必要的自定义函数调用过程。在次数很多的循环体内(包括在cellfun,

    arrayfun等实际上蕴含循环的函数)形成长调用链,会带来很大的开销。

    在调用函数时,首选built-in

    function,然后是普通的m-file函数,然后才是function handle或者anonymous

    function。在使用function handle或者anonymous

    function作为参数传递时,如果该函数被调用多次,最好先用一个变量接住,再传入该变量。这样,可以有效避免重复的解析过程。

    在可能的情况下,使用numeric array或者struct

    array,它们的效率大幅度高于cell

    array(几十倍甚至更多)。对于struct,尽可能使用普通的域(字段,field)访问方式,在非效率关键,执行次数较少,而灵活性要求较高的代码中,可以考虑使用动态名称的域访问。

    虽然object-oriented从软件工程的角度更为优胜,而且object的使用很多时候很方便,但是MATLAB目前对于OO的实现效率很低,在效率关键的代码中应该慎用objects。

    如果需要设计类,应该尽可能采用普通的property,而避免灵活但是效率很低的dependent

    property。如非确实必要,避免重载subsref和subsasgn函数,因为这会全面接管对于object的接口调用,往往会带来非常巨大的开销(成千上万倍的减慢),甚至使得本来几乎不是问题的代码成为性能瓶颈。

    展开全文
  • 【转】提高MATLAB运行效率

    万次阅读 多人点赞 2016-05-04 13:35:16
    用过Matlab的人都知道,Matlab是一种解释性语言,存在计算速度慢的问题,为了提高程序的运行效率matlab提供了多种实用工具及编码技巧。   1. 循环矢量化 Matlab是为矢量和矩阵操作而设计的,因此,可以...
  • C++/Python/Matlab执行效率分析

    万次阅读 2017-08-13 02:21:34
    Matlab效率更高。 这为以后选取语言提供了一个很好的参考。 问题起源与对场内期权MC定价时,一步到最后与按天到最后在计算精度上有无差别,镜像问题是FDM定价一步到最后,与按天到最后计算精度上有无差别。我的...
  • 写在前面:Matlab相较于C来说界面更友好,操作起来也更方便,还可以直接做出好看的图,工科学生掌握一下可以避免跟C“苦苦纠缠”;迭代方法用的高斯赛德尔,其实还有其他一些数值计算方法,我还保存着一部分,会在...
  • 使用Matlab运行程序的时候,我们经常需要知道或比较不同程序的具体运行时间,...经常我们需要计算我们程序到底运行多长时间,这样可以比较程序的执行效率。当然这个对于只有几秒钟的小程序没有什么意义,但是对于大程序
  • 为了提高MATALB的滤波效率,采用了spmd并行计算来提高滤波速度。针对大数据分块后滤波的不连续问题,每次分段滤波时,设置每段overlap值等于滤波器中寄存器的个数。下面直接上代码。function作为对比,下面是单进程...
  • 经常我们需要计算我们程序到底运行多长时间,这样可以比较程序的执行效率。当然这个对于只有几秒钟的小程序没有什么意义,但是对于大程序就有很重要的意义了。  下面我们就说说Matlab计算程序运行时间的三种常用...
  • 如果这里概念不清楚,请看这https://zhuanlan.zhihu.com/p/100595803​zhuanlan.zhihu.com如果你能够熟练的完成上述要求,那么今天我们要使用的MATLAB工具箱,用来大幅度的节省你的时间和精力,提高仿真效率。...
  • ROC曲线 ROC-接收器工作特性。 ROC图是组织分类器并可视化其...该函数计算并绘制经典的ROC曲线以及灵敏度,特异性和效率曲线(请参见截图)。 该功能将显示6个截止点: 1)最大灵敏度 2)最大特异性 3)具有成本效..
  • matlab并行计算之parfor

    2020-09-05 18:31:01
    在运行MATLAB程序的时候,当循环较多和计算量较大时,程序运行时间会非常长。经过查阅MATLAB文档发现,使用parfor循环,可以进行多线程并行运算,大大提高计算的...提高matlab代码的执行效率,是很多码农们的迫切愿望.
  • matlab开发-平方效率测定。r square是计算r平方(确定系数)的简单程序。
  • 这里我想告诉大家 MATLAB 很重要的功能 – Profiler – 能够帮助你很快找到程序的问题,然后可以很快的提高程序的效率。当你写完一个程序之后,首先要试着运行。试着运行的时候有两个目的,一个是验证是否正确(这个...
  • 这个程序准确点说并不是效率计算,只是给出了一个波形,没什么实际意义,大家看着玩呗。 function jieguo=xiaolv(san_n,mul) % 输入参数: % san_n:每周期中载波个数; % mul :调制系数,即三角波最大值/正弦...
  • matlab开发-JACCARD效率和理想化。计算JACCard系数和共现矩阵
  • matlab计算程序运行时间

    万次阅读 2017-05-08 13:52:48
    程序遇到tic时Matlab自动开始计时,运行到toc时自动计算此时与最近一次tic之间的时间。实例%test脚本。通过tic,toc命令直接输出程序运行时间。 tic pause(1) t1=toc pause(1) t2=toc执行结果>> test 时间已过 1.000
  • 使用MATLAB将不同抽取因子组合效率进行可视化。可以看出当D1=32,D2=6时,乘法运算量最低,效率最高。 二级抽取乘法运算量 设计二级抽取的滤波器,采用凯泽窗,通过MATLAB编写代码,滤波器频率响应的图像如图所示。 ...
  • Matlab使用:多核并行计算

    万次阅读 2018-11-25 15:17:19
    Matlab使用:多核并行计算 Matlab程序的运行效率,很大程度上决定着科研工作的效率。如果能把循环转变为矩阵运算无疑是最高效的,但实际使用的过程中经常碰到不得不循环的情况。如果循环次数很多,运行速度就会大大...
  • MATLABMATLAB与科学计算与科学计算与科学计算与科学计算 王沫然王沫然 教授教授 清华大学航天学院 清华大学航天学院 mrwang@ 2012/11/20 回顾和思考回顾和思考回顾和思考回顾和思考 回顾非线性方程/组的求解方法 ...
  • 1、GPU与CPU结构上的对比 2、GPU能加速我的应用程序吗? 3、GPU与CPU在计算效率上的对比 ...6、示例Matlab代码——GPU计算与CPU计算效率的对比 1、GPU与CPU结构上的对比 原文: Multicore machines and hyp...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 598
精华内容 239
关键字:

matlab计算效率

matlab 订阅