关于matlab代码矢量化的理解

forcal

代码矢量化是matlab的精髓，其基本特点是运行速度快和代码简洁，它是如何实现的？
' i/ M2 P- m( Z3 M3 [  s0 U
按我的理解，代码矢量化的本质就是设计专门的函数对数组元素集中运算，这样可提高运行速度，同时兼有代码简洁的特点。
5 e6 u& a, [, o+ R2 w( [$ ]- W
6 o: k7 d$ u5 ~对matlab的理解比较肤浅，但也确实看不出有更深意义的东西，望解惑。

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

forcal

参考：http://bbs.emath.ac.cn/thread-2727-1-1.html

1 E3 Y% t" L0 p( [" U* T关于最快速的矩阵乘实现的讨论。
1 V& v/ q  G' [/ C

forcal

参考：http://bbs.emath.ac.cn/thread-2709-1-1.html
  R4 l+ ~1 O+ l4 B
6 ]' n; Q% I% j5 I" [; U我在多个帖子中有不同说明，但将这些说明再集中到一个帖子中比较麻烦，大家相互参考一下，看能否把这个问题解决了。

forcal

讨论有益！以下是在其他论坛的讨论帖子：
7 @0 c( q' f+ Q# G1 v* S: ssimwe：http://forum.simwe.com/thread-952532-1-3.html
csdn：http://topic.csdn.net/u/20101006/21/2ed9e5c1-cc9f-4623-b1c0-ebd5d1b5a98a.html
  [: O$ U/ M! x" Z3 }cadn：http://topic.csdn.net/u/20101010/15/3bcf2fe0-0abd-4c29-b854-b1d007b16863.html

forcal

  v; f5 X7 t( c) h# s
$ e1 x- e( u% P6 i) ?我正在练手设计的FcMath库也打算以矩阵运算为基础，设计一些专门的函数对数组元素集中运算，运行效率确实有所提高（甚至有些涉及矩阵的算法比matlab还快），代码也简洁了，但不知这是不是矢量化？

$ W5 L( u% _$ {2 c5 [8 Z1 F脚本运行效率应该取决于函数调度效率、对象管理效率和函数内部算法的实现。
5 [& y4 t/ e& J! Y9 @" |, n  s
我感觉，matlab的函数调度效率较低，对象管理效率这个不好说，但一些函数内部的设计比较优秀。故有些Forcal代码比matlab快，而有些慢。

! @/ b+ ]: f$ P* @. J以下例子体现了Forcal和matlab的效率差别所在。
0 l# c3 O! z& W& f0 ?' q( o
! b, K4 ^# i4 C+ `这个matlab程序段是网友lin2009 给出的，理论结果是每个元素均为275000。

1. clear all- B- l6 R& ]% p# T1 X
   
2. clc7 [- b, \! s7 n+ |0 G
   
3. tic. N8 Y; T' x- O9 @/ |' K
   
4. k = zeros(5,5); % //生成5×5全0矩阵
   5 q1 ?$ n6 f. w- Y  k( ]\" ?' e
5. % 循环计算以下程序段1000 00次：$ S2 ^1 G- k. @/ I
   
6. for m = 1:1000 003 P) z5 {2 j* b5 I0 a9 u
   
7.     a = rand(5,7);+ I2 x) i  R\" t\" Z9 z
   
8.     b = rand(7,5);%//生成5×7矩阵a,7×5矩阵b,用0～1之间的随机数初始化) Z$ t: Z/ v. I4 t
   
9.     k = k + a * b + a(1:5, 2:6) * b(2:6, 1:5) - a(:, 7) * b(3, :);$ o/ E) i  W6 d- |( ]
   
10. end
    : ?1 |$ w( L, K) A
11. k7 V; s! I* v0 |
    
12. toc

复制代码

( M! R6 x/ Z% ^9 ^: jForcal代码1：运行稍快的代码，比matlab约快10％吧？
/ V- S. _, m. r5 A

2. !using["math","sys"];9 Z* n9 X! P0 l$ X! ~# Y
   
4. mvar:
   
5. 0 P( Z6 g* K* I) A3 c$ @8 T
6. t0=clock(),+ `1 J/ Y& g, L7 ?; {2 K
   
8. oo{k=zeros[5,5]},: g) r3 I1 l8 F7 T7 e4 t/ C8 X
   
10. i=0,((i++)<100000).while{! I9 o: R% W' i
    
12.   oo{; b8 [8 _4 o: g( z\\" p
    
14.     a=rand[5,7], b=rand[7,5],
    
15. + m5 r7 U9 t3 b' W
16.     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)]
    
17. ( P& O! A  ~) P( H
18.   }
    
19. ) ~  l) e* C. @! Z. J/ X/ U
20. },$ l. O4 w3 D\\" K+ B; @) m- x3 H5 \
    
22. k.outm(),( ~+ |+ t. u$ a\\" h  F
    
24. [clock()-t0]/1000;

!using["math","sys"];<font class=\\"jammer\\">9 Z* n9 X! P0 l$ X! ~# Y</font><br /> mvar:<br /> <span style=\\"display:none\\">0 P( Z6 g* K* I) A3 c$ @8 T</span> t0=clock(),<font class=\\"jammer\\">+ `1 J/ Y& g, L7 ?; {2 K</font><br /> oo{k=zeros[5,5]},<font class=\\"jammer\\">: g) r3 I1 l8 F7 T7 e4 t/ C8 X</font><br /> i=0,((i++)<100000).while{<font class=\\"jammer\\">! I9 o: R% W' i</font><br />   oo{<font class=\\"jammer\\">; b8 [8 _4 o: g( z\\" p</font><br />     a=rand[5,7], b=rand[7,5],<br /> <span style=\\"display:none\\">+ m5 r7 U9 t3 b' W</span>     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)]<br /> <span style=\\"display:none\\">( P& O! A  ~) P( H</span>   }<br /> <span style=\\"display:none\\">) ~  l) e* C. @! Z. J/ X/ U</span> },<font class=\\"jammer\\">$ l. O4 w3 D\\" K+ B; @) m- x3 H5 \</font><br /> k.outm(),<font class=\\"jammer\\">( ~+ |+ t. u$ a\\" h  F</font><br /> [clock()-t0]/1000;

在我的电脑上运行时间为3.344秒。
4 H+ I( ^2 K) b
Forcal代码2：比较好看些的代码，似乎也比matlab稍快吧？

2. !using["math","sys"];1 S' e7 @6 Z; r) d( W2 d
   
4. (:t0,k,i,a,b)=
   
5. 4 U& e! ?, S3 u$ l7 F
6. {$ ~4 `2 _+ `# U; H5 X6 n
   
8.   t0=clock(),  Q# `/ Q7 e( ?7 p4 J
   
10.   oo{k=zeros[5,5]},
    
11. 4 m. ]  ~1 ^\\" F- w( M9 f) F
12.   i=0,((i++)<100000).while{
    
13. 2 S& n* T) h6 W7 _* k5 B
14.     oo{5 w8 T4 r1 l$ N
    
16.       a=rand[5,7], b=rand[7,5],* u( W) g% `% g1 ^+ L* L3 y
    
18.       k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)-a(neg:6)*b(3:neg)! {7 W- n3 L3 q2 N8 H8 Z8 C) b2 W! S
    
20.     }* b$ ~# M8 w6 j! ?
    
22.   },! ?% k; x0 E6 m8 q6 m+ o
    
24.   k.outm(),
    
25. . q\\" ]8 B# r- ~1 b
26.   [clock()-t0]/1000
    
27. 4 d$ E, X' }6 z, o3 y\\" S4 e
28. };

!using["math","sys"];<font class=\\"jammer\\">1 S' e7 @6 Z; r) d( W2 d</font><br /> (:t0,k,i,a,b)=<br /> <span style=\\"display:none\\">4 U& e! ?, S3 u$ l7 F</span> {<font class=\\"jammer\\">$ ~4 `2 _+ `# U; H5 X6 n</font><br />   t0=clock(),<font class=\\"jammer\\">  Q# `/ Q7 e( ?7 p4 J</font><br />   oo{k=zeros[5,5]},<br /> <span style=\\"display:none\\">4 m. ]  ~1 ^\\" F- w( M9 f) F</span>   i=0,((i++)<100000).while{<br /> <span style=\\"display:none\\">2 S& n* T) h6 W7 _* k5 B</span>     oo{<font class=\\"jammer\\">5 w8 T4 r1 l$ N</font><br />       a=rand[5,7], b=rand[7,5],<font class=\\"jammer\\">* u( W) g% `% g1 ^+ L* L3 y</font><br />       k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)-a(neg:6)*b(3:neg)<font class=\\"jammer\\">! {7 W- n3 L3 q2 N8 H8 Z8 C) b2 W! S</font><br />     }<font class=\\"jammer\\">* b$ ~# M8 w6 j! ?</font><br />   },<font class=\\"jammer\\">! ?% k; x0 E6 m8 q6 m+ o</font><br />   k.outm(),<br /> <span style=\\"display:none\\">. q\\" ]8 B# r- ~1 b</span>   [clock()-t0]/1000<br /> <span style=\\"display:none\\">4 d$ E, X' }6 z, o3 y\\" S4 e</span> };

在我的电脑上运行时间为3.579秒。

例子2：
* ?. {/ ?$ |# o$ C. b一段程序的Forcal实现：
" Z/ D! o+ p, G/ m

1. //用C++代码描述为：* e8 ~9 m3 l# x: X. T
   
2. s=0.0;  , P5 q/ X1 v$ q6 P# X9 g
   
3. for(x=0.0;x<=1.0;x=x+0.0011)  / ~6 N, c  p5 b6 C$ V
   
4. {
   + Y4 ~/ z3 ]' \+ y& ?* h( X  b8 O
5.   for(y=1.0;y<=2.0;y=y+0.0009)$ X* K0 u. w) C
   
6.   {$ V$ n1 F: j) n# n\" \% A; {9 }
   
7.   s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   5 F' i- @8 ?8 O( s
8.   }0 x' Z3 ~0 s\" H* q: C9 E2 K
   
9. }  

复制代码

结果：
2 P- t( _$ H9 O1008606.64947441
5 U  m! o. m% `: C: N0.609 //时间
) p1 _" Z1 @4 [! ?( D2 F  u
6 R8 F8 Z& _1 B/ f- x这个matlab程序段是网友yycs001给出的。
) ]0 Z1 a, F" g% K. w

1. %file speedtest.m2 p0 U( X- u8 N' p: w* b
   
2. function speedtest8 U; p, a6 q' \7 n\" ]
   
3. format long
   ( c' m: Z6 d0 `+ B
4. tic
   * I! ]9 R, c: t- Q+ Z
5. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);- l$ e: h3 X% i* i2 p$ E
   
6. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))
   / \) h- m- K: P\" `; ]& _9 C4 f/ V
7. toc

复制代码

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

2. !using["math","sys"];
   
3. ; J9 S/ c* H, p+ j7 L
4. mvar:
   
5. $ C) k! F/ V! B# e\\" \  R
6. t=clock(),5 T) S0 m8 ~; b
   
8. oo{
   
9. 8 _& D& r$ I- k$ |& b
10.   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
11. , T\\" I( \& J3 Q7 t) m
12.   Sum[Cos(rn(1)-Sin(rn(1.2)*x^(y/rn(2))+Cos(rn(1)-Sin(rn(1.2)*y^(x/rn(2)))))),0]
    
13. + P) X; R, O$ T# v9 w: H% r
14. };7 g' a8 y6 g. ?: w, t
    
16. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">; J9 S/ c* H, p+ j7 L</span> mvar:<br /> <span style=\\"display:none\\">$ C) k! F/ V! B# e\\" \  R</span> t=clock(),<font class=\\"jammer\\">5 T) S0 m8 ~; b</font><br /> oo{<br /> <span style=\\"display:none\\">8 _& D& r$ I- k$ |& b</span>   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">, T\\" I( \& J3 Q7 t) m</span>   Sum[Cos(rn(1)-Sin(rn(1.2)*x^(y/rn(2))+Cos(rn(1)-Sin(rn(1.2)*y^(x/rn(2)))))),0]<br /> <span style=\\"display:none\\">+ P) X; R, O$ T# v9 w: H% r</span> };<font class=\\"jammer\\">7 g' a8 y6 g. ?: w, t</font><br /> [clock()-t]/1000;

结果：
; a! L" U: \) p) V1008606.64947441
/ ~. |% h1 B2 S# r0 G0.625 //时间
4 H7 Y$ D7 C, j
& g% M1 @0 T( {' y& B5 k( B- d或者这个，与上面效率差别不大：
& {2 [2 S& `/ _- h

2. !using["math","sys"];: H0 C+ I3 E/ v6 |& H, x
   
4. mvar:
   
5. 1 |6 _3 U\\" F% n& U6 j; z+ x% Q2 r
6. t=clock(),& W9 \% S0 o3 m
   
8. oo{! C1 Z! s! x* G. L; O9 |8 H
   
10.   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],: m# F6 |7 u5 n# H4 i\\" n) y: q
    
12.   Sum[Sum[Cos(rn(1)-Sin(rn(1.2)*x^(y/rn(2))+Cos(rn(1)-Sin(rn(1.2)*y^(x/rn(2))))))]]8 o& Q( \; h+ r& r4 S3 S
    
14. };
    
15. ' T% V. L/ S5 S\\" J
16. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">: H0 C+ I3 E/ v6 |& H, x</font><br /> mvar:<br /> <span style=\\"display:none\\">1 |6 _3 U\\" F% n& U6 j; z+ x% Q2 r</span> t=clock(),<font class=\\"jammer\\">& W9 \% S0 o3 m</font><br /> oo{<font class=\\"jammer\\">! C1 Z! s! x* G. L; O9 |8 H</font><br />   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">: m# F6 |7 u5 n# H4 i\\" n) y: q</font><br />   Sum[Sum[Cos(rn(1)-Sin(rn(1.2)*x^(y/rn(2))+Cos(rn(1)-Sin(rn(1.2)*y^(x/rn(2))))))]]<font class=\\"jammer\\">8 o& Q( \; h+ r& r4 S3 S</font><br /> };<br /> <span style=\\"display:none\\">' T% V. L/ S5 S\\" J</span> [clock()-t]/1000;

+ o, D! H' k, r4 R/ IForcal代码2：求和函数sum，非矢量化代码
. ^4 c% z# |+ y1 O: [

1. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   0 H) G- T, y( ]1 M
2. sum["f",0,1,0.0011,1,2,0.0009];

复制代码

结果：
1008606.64947441
0.719 //时间
  S2 k3 p2 q" d& c9 ^+ d
, x' `0 a/ p6 s/ {Forcal代码3：while循环

1. mvar:4 J5 x9 K' a4 U$ j2 i% F8 c
   
2. t=sys::clock();) ]0 V4 t/ @/ X9 y) I
   
3. s=0,x=0,
   ; d8 m+ U4 r6 T2 m0 h# E
4. while{x<=1,  //while循环算法；
   - B  y1 F/ q  J+ w& C* l
5.    y=1,
   \" \( t* H# i& L  ^) S
6.    while{y<=2,
   : `: q2 L3 t6 d
7.        s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))), ) R9 s5 \. Q1 h% r
   
8.        y=y+0.0009
   ' x; k* }: q( \6 ]- T0 X2 V
9.       }, 9 ?( ^) C% f! c8 V2 o\" @
   
10.    x=x+0.0011 . s8 @: m% K; f\" ?6 s& L$ Y& D0 N
    
11. }, & E  h4 m. n% }) q$ J
    
12. s;2 l\" h1 M! ]; d; Q+ j
    
13. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
" R! j. ]. `( p  }9 b" {7 ^8 Z$ u0.734 //时间
5 T6 w9 o# Q7 M4 y3 F$ w9 h
大家可下载OpenFC进行测试：http://www.forcal.net/xiazai/forcal9/openfc32w.rar

4 Z2 `& `) k* r$ y/ E1 O9 q注意Forcal的矢量化代码第一次运行有时效率较低。
. x- l& U: b- _+ X
1 f1 E# [$ N" o例子1中Forcal和matlab都是矢量化代码，但matlab跑不过Forcal。该例子的特点是函数调用频繁，临时变量生成多，但矩阵很小，矩阵的各种函数运行时耗时较少。故说明Forcal函数调用＋变量管理效率优于matlab。
3 @  Y5 U% c& N
例子2中Forcal的矢量化代码是最快的，但与matlab的矢量化代码相比仍有差距。该例子的特点是函数调用少，临时变量也少，但矩阵大。故说明Forcal的各种矩阵函数Sin、Cos及矩阵的加减运算等函数的内部设计不及matlab。
9 C  H7 V5 G8 W
如能在函数内部设计上下点功夫，例子2超越matlab也是可能的。在这方面，期待高手们的指点。
# Z' G, Y) n* @8 Z$ X. \
8 {5 G6 Q% W+ K- F' ]9 u. ?# u/ M如果例子2速度也超越了matlab ，matlab矢量化的神秘面纱就揭开了。

  Y: E/ {( E- z9 a+ J1 ^顺便说一下，例子1如果用C++的运算符重载来实现，速度将比Forcal慢一些，也就是说，在涉及运算符重载时，脚本的效率有时比C++还要高些。
* A% c& ?# ]$ v

master_math

,,,顶顶更健康！

wznzy0822

顶。。。。。。。。。。。。。。。。。。。。。。

zhwqqiangge

??????????????????

qbist

我 正在  学  Matlab  软件！~~