关于matlab代码矢量化的理解

forcal

代码矢量化是matlab的精髓，其基本特点是运行速度快和代码简洁，它是如何实现的？

按我的理解，代码矢量化的本质就是设计专门的函数对数组元素集中运算，这样可提高运行速度，同时兼有代码简洁的特点。
, `2 N1 @. s7 Q" b) h) u! r
对matlab的理解比较肤浅，但也确实看不出有更深意义的东西，望解惑。
- M- x( i- B  Z( D) n0 K
大家有什么看法，愿畅所欲言。

qbist

我 正在  学  Matlab  软件！~~

zhwqqiangge

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

wznzy0822

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

master_math

,,,顶顶更健康！

forcal

; |! w0 r2 N9 X- R  G6 a' A
我正在练手设计的FcMath库也打算以矩阵运算为基础，设计一些专门的函数对数组元素集中运算，运行效率确实有所提高（甚至有些涉及矩阵的算法比matlab还快），代码也简洁了，但不知这是不是矢量化？
, l: J' A' L# c. S. A1 u  |
脚本运行效率应该取决于函数调度效率、对象管理效率和函数内部算法的实现。
! |7 j5 s0 N) c3 v( R6 z
4 {# n5 b' [  @* R- c8 o我感觉，matlab的函数调度效率较低，对象管理效率这个不好说，但一些函数内部的设计比较优秀。故有些Forcal代码比matlab快，而有些慢。

以下例子体现了Forcal和matlab的效率差别所在。

这个matlab程序段是网友lin2009 给出的，理论结果是每个元素均为275000。

1. clear all
   , s2 A# l' ^, ^$ E
2. clc
   * [  y, I* r; T0 b
3. tic\" ?: I+ {3 U9 j( y; Y' k7 R\" t. m
   
4. k = zeros(5,5); % //生成5×5全0矩阵
   ) d0 c9 o0 H5 e
5. % 循环计算以下程序段1000 00次：  z; J% F8 b3 ^8 X& M$ y- g
   
6. for m = 1:1000 00
   8 G; l- D+ L. s4 k. E3 T; H' J
7.     a = rand(5,7);
   2 k+ R9 t3 a6 l+ y$ b' L
8.     b = rand(7,5);%//生成5×7矩阵a,7×5矩阵b,用0～1之间的随机数初始化+ M& z7 T/ f2 J  N2 M
   
9.     k = k + a * b + a(1:5, 2:6) * b(2:6, 1:5) - a(:, 7) * b(3, :);# k& L+ k: L& w. v
   
10. end\" a: F\" j1 R2 f) _- h
    
11. k
    3 b1 D% E; A' O* l\" d
12. toc

复制代码

Forcal代码1：运行稍快的代码，比matlab约快10％吧？
: p0 [9 v2 r, u' R" d

2. !using["math","sys"];9 p! i1 T, {* X: H0 U- a+ H
   
4. mvar:
   
5. 6 W- c* _' D- k
6. t0=clock(),
   
7. $ z( [\\" S4 o8 V1 w* K
8. oo{k=zeros[5,5]},2 ^1 r9 ?) e2 \$ w- V2 I$ Q
   
10. i=0,((i++)<100000).while{' q  K4 u8 u( _
    
12.   oo{
    
13.   F( [: ]. v! {0 v
14.     a=rand[5,7], b=rand[7,5],3 \: ?* P, L\\" n0 Y$ v
    
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)]+ E1 f* M5 b% R& r3 h\\" ?
    
18.   }1 X% B% E, n. q! f
    
20. },) ?8 Y- J' D( p\\" P) c! f
    
22. k.outm(),$ D$ ^( s3 U: T2 ?  v6 r
    
24. [clock()-t0]/1000;

!using["math","sys"];<font class=\\"jammer\\">9 p! i1 T, {* X: H0 U- a+ H</font><br /> mvar:<br /> <span style=\\"display:none\\">6 W- c* _' D- k</span> t0=clock(),<br /> <span style=\\"display:none\\">$ z( [\\" S4 o8 V1 w* K</span> oo{k=zeros[5,5]},<font class=\\"jammer\\">2 ^1 r9 ?) e2 \$ w- V2 I$ Q</font><br /> i=0,((i++)<100000).while{<font class=\\"jammer\\">' q  K4 u8 u( _</font><br />   oo{<br /> <span style=\\"display:none\\">  F( [: ]. v! {0 v</span>     a=rand[5,7], b=rand[7,5],<font class=\\"jammer\\">3 \: ?* P, L\\" n0 Y$ v</font><br />     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)]<font class=\\"jammer\\">+ E1 f* M5 b% R& r3 h\\" ?</font><br />   }<font class=\\"jammer\\">1 X% B% E, n. q! f</font><br /> },<font class=\\"jammer\\">) ?8 Y- J' D( p\\" P) c! f</font><br /> k.outm(),<font class=\\"jammer\\">$ D$ ^( s3 U: T2 ?  v6 r</font><br /> [clock()-t0]/1000;

在我的电脑上运行时间为3.344秒。
7 l, X* @) S* E0 p: R, v! ^5 C$ J
" h7 p% D* M0 r1 l/ X. I2 x0 AForcal代码2：比较好看些的代码，似乎也比matlab稍快吧？

2. !using["math","sys"];) l& {6 t4 r9 E7 W
   
4. (:t0,k,i,a,b)=
   
5. - r* e7 n: C: @; y\\" V
6. {# Y. ?* H1 p# R8 m8 Y
   
8.   t0=clock(),
   
9. / y. b, N7 X8 X, c
10.   oo{k=zeros[5,5]},- J, M% t+ H+ {\\" }% \. @
    
12.   i=0,((i++)<100000).while{
    
13. / w\\" y/ d8 ~: s, A  q& X
14.     oo{
    
15. ' m' \+ L5 |; c# X, f$ e; f
16.       a=rand[5,7], b=rand[7,5],* l  t) F: ^6 ]. p9 n, X. F
    
18.       k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)-a(neg:6)*b(3:neg)7 @& f8 ~\\" g7 z' x* Z2 {
    
20.     }
    
21. ! d3 m% K# Z0 W0 r7 c
22.   },
    
23. 6 H4 F8 |$ Q3 ?, Y0 x  |! l
24.   k.outm(),3 H1 ]9 Q: n8 h* u( _/ h
    
26.   [clock()-t0]/10007 @0 |  o\\" U9 p2 \' s* p# x
    
28. };

!using["math","sys"];<font class=\\"jammer\\">) l& {6 t4 r9 E7 W</font><br /> (:t0,k,i,a,b)=<br /> <span style=\\"display:none\\">- r* e7 n: C: @; y\\" V</span> {<font class=\\"jammer\\"># Y. ?* H1 p# R8 m8 Y</font><br />   t0=clock(),<br /> <span style=\\"display:none\\">/ y. b, N7 X8 X, c</span>   oo{k=zeros[5,5]},<font class=\\"jammer\\">- J, M% t+ H+ {\\" }% \. @</font><br />   i=0,((i++)<100000).while{<br /> <span style=\\"display:none\\">/ w\\" y/ d8 ~: s, A  q& X</span>     oo{<br /> <span style=\\"display:none\\">' m' \+ L5 |; c# X, f$ e; f</span>       a=rand[5,7], b=rand[7,5],<font class=\\"jammer\\">* l  t) F: ^6 ]. p9 n, X. F</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 @& f8 ~\\" g7 z' x* Z2 {</font><br />     }<br /> <span style=\\"display:none\\">! d3 m% K# Z0 W0 r7 c</span>   },<br /> <span style=\\"display:none\\">6 H4 F8 |$ Q3 ?, Y0 x  |! l</span>   k.outm(),<font class=\\"jammer\\">3 H1 ]9 Q: n8 h* u( _/ h</font><br />   [clock()-t0]/1000<font class=\\"jammer\\">7 @0 |  o\\" U9 p2 \' s* p# x</font><br /> };

在我的电脑上运行时间为3.579秒。
, D, q. S: g+ E8 O8 i1 }
例子2：
" l1 q+ x0 q( O. l0 A8 q4 C一段程序的Forcal实现：
) g( D  `$ D8 ^9 P5 d3 v0 e% ]0 o

1. //用C++代码描述为：
   - I0 P* M0 b$ W$ Q9 k& K
2. s=0.0;  / m. f8 q3 r8 T, T3 u9 |
   
3. for(x=0.0;x<=1.0;x=x+0.0011)  + J1 S' y( S* G; w1 }# h
   
4. {1 ?0 p/ Q2 B0 r1 J- z- c* W
   
5.   for(y=1.0;y<=2.0;y=y+0.0009)
   2 o' {& a4 g9 I7 {/ I: Q
6.   {; x+ z/ e) L* L2 p+ ^1 ]6 N3 ?
   
7.   s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   2 N# F# R2 P' n; ~: d7 g
8.   }
   5 M9 ^. n& Q; C5 W
9. }  

复制代码

结果：
9 K/ Z0 F' x0 J" L! W, K$ e1008606.64947441
0.609 //时间

: l, x6 Z  V: R! h这个matlab程序段是网友yycs001给出的。

1. %file speedtest.m9 W% k% A+ e  S! R
   
2. function speedtest# A6 w8 S+ U+ H- `: H2 z
   
3. format long
   6 O& J$ w: _- Z3 n( S
4. tic4 V3 o7 Z* w; N, b4 e
   
5. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   ( W' y. L7 Y7 X6 V* r
6. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))9 S\" m* l5 f: M1 U, O1 N8 v: H
   
7. toc

复制代码

! L* m- z' H; e# S$ @Forcal代码1：**数组求和函数Sum，完全矢量化的代码
9 F' O! u  h+ i! y8 B( D

2. !using["math","sys"];1 Z4 M- m0 r6 W2 ]7 }
   
4. mvar:* X% Y; i2 {4 @. y: @9 |1 e
   
6. t=clock(),
   
7. 9 r, C, a! z1 y5 K* y
8. oo{
   
9. % ], b; k2 N& d. x3 R  d. h
10.   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],! D0 Z* a9 K' H
    
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. ' V, L' R; x- n+ ]) s5 }0 O( |9 m5 K' l
14. };
    
15. 4 c% P* x! S3 @/ s, ^) w3 [
16. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">1 Z4 M- m0 r6 W2 ]7 }</font><br /> mvar:<font class=\\"jammer\\">* X% Y; i2 {4 @. y: @9 |1 e</font><br /> t=clock(),<br /> <span style=\\"display:none\\">9 r, C, a! z1 y5 K* y</span> oo{<br /> <span style=\\"display:none\\">% ], b; k2 N& d. x3 R  d. h</span>   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">! D0 Z* a9 K' H</font><br />   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\\">' V, L' R; x- n+ ]) s5 }0 O( |9 m5 K' l</span> };<br /> <span style=\\"display:none\\">4 c% P* x! S3 @/ s, ^) w3 [</span> [clock()-t]/1000;

结果：
0 @5 r+ K4 h0 x6 \, ~8 ^: o1008606.64947441
7 o1 q# Z9 G% J% H0.625 //时间
* \- a) Z9 d; C" }7 x
或者这个，与上面效率差别不大：
. e% E+ I. E$ k0 u. g" h

2. !using["math","sys"];0 m+ q: [2 z$ m
   
4. mvar:
   
5. 3 q( b5 [- p* b  w, ~
6. t=clock(),
   
7. 3 m$ z7 v2 J+ A0 D
8. oo{, |  s9 y# m/ C5 L& o
   
10.   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
11. 7 i, J/ n4 {9 Q! Z
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))))))]]% x1 p- E- W0 q2 J( z$ V6 f\\" e
    
14. };- `4 w( X5 y; N; i9 J) |; [
    
16. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">0 m+ q: [2 z$ m</font><br /> mvar:<br /> <span style=\\"display:none\\">3 q( b5 [- p* b  w, ~</span> t=clock(),<br /> <span style=\\"display:none\\">3 m$ z7 v2 J+ A0 D</span> oo{<font class=\\"jammer\\">, |  s9 y# m/ C5 L& o</font><br />   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">7 i, J/ n4 {9 Q! Z</span>   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\\">% x1 p- E- W0 q2 J( z$ V6 f\\" e</font><br /> };<font class=\\"jammer\\">- `4 w( X5 y; N; i9 J) |; [</font><br /> [clock()-t]/1000;

Forcal代码2：求和函数sum，非矢量化代码
! L6 u# v; Q  L

1. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   . m( Y0 m1 O2 c' g6 k% K
2. sum["f",0,1,0.0011,1,2,0.0009];

复制代码

结果：
1008606.64947441
# z* g7 m5 f- k/ x5 [0.719 //时间

Forcal代码3：while循环

1. mvar:: Y3 `! V8 s: [& l8 }5 \& ]$ p
   
2. t=sys::clock();8 ^* B/ W  Z$ Q4 C
   
3. s=0,x=0,
   . M/ k) [8 @1 M1 P3 u: a4 w
4. while{x<=1,  //while循环算法； ) k/ N+ _  a5 ?* r' M3 Y; J- B
   
5.    y=1,
   1 t- K4 c/ L+ k% W0 z
6.    while{y<=2,
   / \, s2 n2 K, S% B+ X9 z6 w% ]
7.        s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
   7 d% y2 C3 c6 q3 L
8.        y=y+0.0009 2 ~$ i% ^# l; u, j: g# o( S
   
9.       }, + H) h# }4 M+ v- b0 B* b
   
10.    x=x+0.0011 * K& E8 _. I! d: W  l
    
11. },
    : `, Y3 s5 L/ U+ K9 L
12. s;1 R7 W# a+ x9 D* R
    
13. [sys::clock()-t]/1000;

复制代码

结果：
8 u6 l* r6 ~* n. c1008606.64947441
5 d& p- r* ~* n8 k" |: ^0.734 //时间
/ J& o- O$ C% E. ?# L) N' M
大家可下载OpenFC进行测试：http://www.forcal.net/xiazai/forcal9/openfc32w.rar

# M; c  ]# T  B: C1 h& h- U. W注意Forcal的矢量化代码第一次运行有时效率较低。
, H. P5 \8 G# }$ H
4 |0 x% y/ b9 a% w7 _. @- |2 ~6 f例子1中Forcal和matlab都是矢量化代码，但matlab跑不过Forcal。该例子的特点是函数调用频繁，临时变量生成多，但矩阵很小，矩阵的各种函数运行时耗时较少。故说明Forcal函数调用＋变量管理效率优于matlab。

例子2中Forcal的矢量化代码是最快的，但与matlab的矢量化代码相比仍有差距。该例子的特点是函数调用少，临时变量也少，但矩阵大。故说明Forcal的各种矩阵函数Sin、Cos及矩阵的加减运算等函数的内部设计不及matlab。
5 w$ `9 J- Q, b; f
+ W" U3 n2 z. m如能在函数内部设计上下点功夫，例子2超越matlab也是可能的。在这方面，期待高手们的指点。

' o. L. i% r! }0 a如果例子2速度也超越了matlab ，matlab矢量化的神秘面纱就揭开了。

顺便说一下，例子1如果用C++的运算符重载来实现，速度将比Forcal慢一些，也就是说，在涉及运算符重载时，脚本的效率有时比C++还要高些。

forcal

讨论有益！以下是在其他论坛的讨论帖子：
4 O0 g+ }# X2 P5 d! Psimwe：http://forum.simwe.com/thread-952532-1-3.html
csdn：http://topic.csdn.net/u/20101006/21/2ed9e5c1-cc9f-4623-b1c0-ebd5d1b5a98a.html
/ s0 H! g( ^7 ~6 Tcadn：http://topic.csdn.net/u/20101010/15/3bcf2fe0-0abd-4c29-b854-b1d007b16863.html

forcal

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

) j% ~4 U' }4 G3 R3 P1 |1 M$ X8 R# F我在多个帖子中有不同说明，但将这些说明再集中到一个帖子中比较麻烦，大家相互参考一下，看能否把这个问题解决了。

forcal

参考：http://bbs.emath.ac.cn/thread-2727-1-1.html
& E4 z$ \  I9 T, ~
关于最快速的矩阵乘实现的讨论。
  e' S' k5 }, ~# R- }