数学建模社区-数学中国

标题: 关于matlab代码矢量化的理解 [打印本页]

---

作者: forcal    时间: 2010-10-5 09:32
标题: 关于matlab代码矢量化的理解
代码矢量化是matlab的精髓，其基本特点是运行速度快和代码简洁，它是如何实现的？

1 ]: O. F/ O1 [" @6 O) g& ?. {3 U按我的理解，代码矢量化的本质就是设计专门的函数对数组元素集中运算，这样可提高运行速度，同时兼有代码简洁的特点。
$ Z" @4 }- N! s" d
1 Q2 z/ d4 d  r2 D: Z8 ~/ _对matlab的理解比较肤浅，但也确实看不出有更深意义的东西，望解惑。

大家有什么看法，愿畅所欲言。

---

作者: qbist    时间: 2010-10-5 09:52
我 正在  学  Matlab  软件！~~

---

作者: zhwqqiangge    时间: 2010-10-5 10:17
??????????????????

---

作者: wznzy0822    时间: 2010-10-5 10:52
顶。。。。。。。。。。。。。。。。。。。。。。

---

作者: master_math    时间: 2010-10-5 18:32
,,,顶顶更健康！

---

作者: forcal    时间: 2010-10-12 21:36


7 n+ J. s- F$ f" _' M# F* n. a我正在练手设计的FcMath库也打算以矩阵运算为基础，设计一些专门的函数对数组元素集中运算，运行效率确实有所提高（甚至有些涉及矩阵的算法比matlab还快），代码也简洁了，但不知这是不是矢量化？
2 L) V- `$ y0 F  F8 v, t0 N& R7 @
% l9 I! l8 K1 [& p0 @' L9 T- H脚本运行效率应该取决于函数调度效率、对象管理效率和函数内部算法的实现。

我感觉，matlab的函数调度效率较低，对象管理效率这个不好说，但一些函数内部的设计比较优秀。故有些Forcal代码比matlab快，而有些慢。

以下例子体现了Forcal和matlab的效率差别所在。
, V5 a, Q3 k7 C) T9 W. T
这个matlab程序段是网友lin2009 给出的，理论结果是每个元素均为275000。

1. clear all
   
2. clc
   
3. tic
   
4. k = zeros(5,5); % //生成5×5全0矩阵
   * @* t% S( c/ p+ f; G- l8 O3 Z
5. % 循环计算以下程序段1000 00次：
   # R# ?% f- [7 o( Q
6. for m = 1:1000 00
   
7.     a = rand(5,7);
   , {( D1 t3 N2 M
8.     b = rand(7,5);%//生成5×7矩阵a,7×5矩阵b,用0～1之间的随机数初始化
   
9.     k = k + a * b + a(1:5, 2:6) * b(2:6, 1:5) - a(:, 7) * b(3, :);
   
10. end
    $ z5 F) Q7 s1 J9 N% a& j% ]
11. k
    
12. toc

复制代码

Forcal代码1：运行稍快的代码，比matlab约快10％吧？

1. !using["math","sys"];
   
2. mvar:
   - V  W! G) z, G3 N- \4 `# _
3. t0=clock(),
   ) r- M7 q* p, M3 R# {- `
4. oo{k=zeros[5,5]},
   
5. i=0,((i++)<100000).while{
   
6.   oo{
   0 a$ p9 z3 ~6 Z
7.     a=rand[5,7], b=rand[7,5],
   
8.     k.oset[k+a*b+a.subg(0,4:1,5)*b.subg(1,5:0,4)-a.subg(neg:6)*b.subg(3:neg)]
   
9.   }
   6 u' V" S) d# ?% K1 f7 o- o
10. },
    4 ?& l! j' x8 ~7 j4 Y
11. k.outm(),
    % U; r+ q$ g# Y" b% t4 _
12. [clock()-t0]/1000;

复制代码

在我的电脑上运行时间为3.344秒。
& ]# Y* D  I/ t# g  O
Forcal代码2：比较好看些的代码，似乎也比matlab稍快吧？
7 B" Q, _& D3 ?, }- ~

1. !using["math","sys"];
   
2. (:t0,k,i,a,b)=
   
3. {
   
4.   t0=clock(),
   + \+ r# j' b$ H8 w( H6 S3 y
5.   oo{k=zeros[5,5]},
   7 h( |% j5 Z) I- \" \8 Y* m
6.   i=0,((i++)<100000).while{
   
7.     oo{
   
8.       a=rand[5,7], b=rand[7,5],
   
9.       k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)-a(neg:6)*b(3:neg)
   
10.     }
    " I% u8 q; J: S: t3 D: F& W6 h
11.   },
    , `- Z# M, r7 M: n8 q
12.   k.outm(),
    
13.   [clock()-t0]/1000
    : k! D1 C( ~/ z$ Y
14. };

复制代码

在我的电脑上运行时间为3.579秒。
  p4 w+ u7 o" Y; T/ _
- h( J6 {; }4 ^* H" K0 q( y& s例子2：
一段程序的Forcal实现：
0 ]0 T  v6 C  c0 J1 Z- e4 h0 [

1. //用C++代码描述为：
   1 b2 l: W- H0 V' F9 a
2. s=0.0;  
   ) N3 Z: b: r0 N, ^1 D* g5 z
3. for(x=0.0;x<=1.0;x=x+0.0011)  
   
4. {
     O# k# Z# n$ [' B& O$ W1 b
5.   for(y=1.0;y<=2.0;y=y+0.0009)
   
6.   {
   7 c+ o5 |" \$ J, P4 K
7.   s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   9 W, n! \3 m: A5 c) s# b0 u
8.   }
   8 [9 U  W- \/ P" {
9. }  

复制代码

结果：
1008606.64947441
) ~0 g, R% S+ Z' p! G0.609 //时间

这个matlab程序段是网友yycs001给出的。
  w3 U7 V. o, w( d2 w

1. %file speedtest.m
   % T6 b$ d) n  K- O4 k- s
2. function speedtest
   # t3 h& X8 `+ Q0 B
3. format long
   " Z* i* H7 i. x7 D5 J& Y1 C
4. tic
   0 x  I- J9 c' p6 o6 _+ t
5. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   
6. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))
   ( P; G7 R% z. B  p
7. toc

复制代码

Forcal代码1：**数组求和函数Sum，完全矢量化的代码

1. !using["math","sys"];
   
2. mvar:
   $ N0 }9 C$ m& |; R2 m
3. t=clock(),
   0 F) s! r1 ^; c+ V2 `
4. oo{
   
5.   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
   
6.   Sum[Cos(rn(1)-Sin(rn(1.2)*x^(y/rn(2))+Cos(rn(1)-Sin(rn(1.2)*y^(x/rn(2)))))),0]
   
7. };
   % M3 u) C2 ~( o  g& d' Q8 l$ \! h: G
8. [clock()-t]/1000;

复制代码

结果：
1 n% n6 \7 ~; G1 q; L+ p& |1008606.64947441
0.625 //时间
# W9 Z+ J6 }* R6 ~% X7 m0 N
或者这个，与上面效率差别不大：
  P* w' j) Y! x

1. !using["math","sys"];
   
2. mvar:
   9 Y, ?3 m% r- N4 ?
3. t=clock(),
   
4. oo{
   3 p: j0 |7 d: Z5 g, m! d1 J# \
5.   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
   5 j8 D) S8 O& d3 R0 |' r! O. g
6.   Sum[Sum[Cos(rn(1)-Sin(rn(1.2)*x^(y/rn(2))+Cos(rn(1)-Sin(rn(1.2)*y^(x/rn(2))))))]]
   
7. };
   
8. [clock()-t]/1000;

复制代码

Forcal代码2：求和函数sum，非矢量化代码

1. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   8 N5 c9 l+ ?% k5 y. m; b& t
2. sum["f",0,1,0.0011,1,2,0.0009];

复制代码

结果：
- H6 r2 i* I9 v# A" C1008606.64947441
0.719 //时间
. q& T) r! }' L2 Z. E; Z
Forcal代码3：while循环
% }6 I* S6 E3 B. y

1. mvar:
   
2. t=sys::clock();
   
3. s=0,x=0,
   
4. while{x<=1,  //while循环算法；
   
5.    y=1,
   
6.    while{y<=2,
   , p3 r6 E/ _( G1 Q
7.        s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
   
8.        y=y+0.0009
   9 z# F  ^' h& i7 l) K$ _' }5 h( i
9.       },
   
10.    x=x+0.0011
    ! i$ {  t' _  U& ~- ~
11. },
    
12. s;
    
13. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
0.734 //时间

5 e* g$ m; i: t  D% d- Q* D# W大家可下载OpenFC进行测试：http://www.forcal.net/xiazai/forcal9/openfc32w.rar

+ N. B8 n- W- [$ l注意Forcal的矢量化代码第一次运行有时效率较低。
& \. ^' y) L- S+ w, U
例子1中Forcal和matlab都是矢量化代码，但matlab跑不过Forcal。该例子的特点是函数调用频繁，临时变量生成多，但矩阵很小，矩阵的各种函数运行时耗时较少。故说明Forcal函数调用＋变量管理效率优于matlab。

; @- }3 C+ B6 m+ W1 l* f! t例子2中Forcal的矢量化代码是最快的，但与matlab的矢量化代码相比仍有差距。该例子的特点是函数调用少，临时变量也少，但矩阵大。故说明Forcal的各种矩阵函数Sin、Cos及矩阵的加减运算等函数的内部设计不及matlab。

如能在函数内部设计上下点功夫，例子2超越matlab也是可能的。在这方面，期待高手们的指点。
/ u) z7 G( `$ {4 O* m3 O, V
" H2 w8 O" ]) \" V" y如果例子2速度也超越了matlab ，matlab矢量化的神秘面纱就揭开了。

$ S) M9 P9 y2 k# t" f6 L( ^2 G顺便说一下，例子1如果用C++的运算符重载来实现，速度将比Forcal慢一些，也就是说，在涉及运算符重载时，脚本的效率有时比C++还要高些。
* r" Y/ p$ A/ j+ V0 M' m( N0 U+ j

---

作者: forcal    时间: 2010-10-16 16:59
讨论有益！以下是在其他论坛的讨论帖子：
simwe：http://forum.simwe.com/thread-952532-1-3.html
csdn：http://topic.csdn.net/u/20101006/21/2ed9e5c1-cc9f-4623-b1c0-ebd5d1b5a98a.html
* D+ L/ J  c$ F  S. b) rcadn：http://topic.csdn.net/u/20101010/15/3bcf2fe0-0abd-4c29-b854-b1d007b16863.html

---

作者: forcal    时间: 2010-10-24 16:59
参考：http://bbs.emath.ac.cn/thread-2709-1-1.html
0 p6 _5 [! v4 J3 }
- f- o; q5 ^' D我在多个帖子中有不同说明，但将这些说明再集中到一个帖子中比较麻烦，大家相互参考一下，看能否把这个问题解决了。

---

作者: forcal    时间: 2010-11-2 18:06
参考：http://bbs.emath.ac.cn/thread-2727-1-1.html

6 v4 R# D" Z& {& n关于最快速的矩阵乘实现的讨论。

---

欢迎光临 数学建模社区-数学中国 (http://www.madio.net/)

Powered by Discuz! X2.5