关于matlab代码矢量化的理解

forcal

代码矢量化是matlab的精髓，其基本特点是运行速度快和代码简洁，它是如何实现的？
: `( |4 t) a. z0 E6 X
按我的理解，代码矢量化的本质就是设计专门的函数对数组元素集中运算，这样可提高运行速度，同时兼有代码简洁的特点。
. y8 h. [( J2 s
对matlab的理解比较肤浅，但也确实看不出有更深意义的东西，望解惑。

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

forcal

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

关于最快速的矩阵乘实现的讨论。
! b0 k% F% g$ w8 g4 u

forcal

参考：http://bbs.emath.ac.cn/thread-2709-1-1.html
0 j, I% L9 u* f. L3 ^
我在多个帖子中有不同说明，但将这些说明再集中到一个帖子中比较麻烦，大家相互参考一下，看能否把这个问题解决了。

forcal

讨论有益！以下是在其他论坛的讨论帖子：
simwe：http://forum.simwe.com/thread-952532-1-3.html
csdn：http://topic.csdn.net/u/20101006/21/2ed9e5c1-cc9f-4623-b1c0-ebd5d1b5a98a.html
$ S% ^9 G+ N' W' E8 h( o* h! D8 a( Ocadn：http://topic.csdn.net/u/20101010/15/3bcf2fe0-0abd-4c29-b854-b1d007b16863.html

forcal

" n! y1 K4 P6 U/ ?  f1 k
: v* S; w: M" h& S3 V% g我正在练手设计的FcMath库也打算以矩阵运算为基础，设计一些专门的函数对数组元素集中运算，运行效率确实有所提高（甚至有些涉及矩阵的算法比matlab还快），代码也简洁了，但不知这是不是矢量化？
) V4 [0 G* c% I8 l  l: X. X
; z* o6 @, T# D( Z3 T* w* ~, p. I- U2 D脚本运行效率应该取决于函数调度效率、对象管理效率和函数内部算法的实现。
& I: K/ G8 R# j* Q7 R; Q
6 }/ h/ D1 D5 E6 z7 w" K! R我感觉，matlab的函数调度效率较低，对象管理效率这个不好说，但一些函数内部的设计比较优秀。故有些Forcal代码比matlab快，而有些慢。

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

( C* m6 s% b& W9 j4 @4 a这个matlab程序段是网友lin2009 给出的，理论结果是每个元素均为275000。

1. clear all
   ; e) m( K* `/ e, h
2. clc
   ) M! c0 b1 {; E3 [+ _; P' d. K
3. tic% F3 W3 J- C8 I2 L7 H3 d( y* g
   
4. k = zeros(5,5); % //生成5×5全0矩阵3 `6 ?2 ?# }. \# h\" \
   
5. % 循环计算以下程序段1000 00次：
   / X$ o: Z( a% l( @1 w4 b/ t
6. for m = 1:1000 00
   \" d5 ]0 B! l& y2 I
7.     a = rand(5,7);
   : q& V3 H- s/ n. q
8.     b = rand(7,5);%//生成5×7矩阵a,7×5矩阵b,用0～1之间的随机数初始化, [\" w) u2 [6 F8 e9 R; j
   
9.     k = k + a * b + a(1:5, 2:6) * b(2:6, 1:5) - a(:, 7) * b(3, :);+ M+ H- m! x; D) c\" S
   
10. end  ~7 a1 Q* U- k7 R# t2 k
    
11. k
    6 b0 h# X- O\" X2 T2 }
12. toc

复制代码

& Z  ~4 y( _$ P7 A) VForcal代码1：运行稍快的代码，比matlab约快10％吧？

2. !using["math","sys"];
   
3. - ?; ?- M! G/ B) v8 z) I
4. mvar:
   
5. \\" \/ _, x! K9 z5 w, w
6. t0=clock(),5 u7 B/ H- B# Y5 c$ m; ?6 b1 m
   
8. oo{k=zeros[5,5]},% A6 Z* @1 Z5 b0 C) c
   
10. i=0,((i++)<100000).while{
    
11. , d1 v- r; {5 c- ]& f4 j! N
12.   oo{& K5 Q5 m. [3 P9 k' |
    
14.     a=rand[5,7], b=rand[7,5],
    
15. / T3 _* y5 A# N' e9 F- m# J
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. : R: v8 j' `2 b\\" U1 I
18.   }9 g5 C5 a\\" {7 _, R  T! u2 t
    
20. },+ E& c/ A3 m3 l  X! e9 T1 W6 K' v
    
22. k.outm(),2 n7 c( G! I0 e. A\\" _# K  i
    
24. [clock()-t0]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">- ?; ?- M! G/ B) v8 z) I</span> mvar:<br /> <span style=\\"display:none\\">\\" \/ _, x! K9 z5 w, w</span> t0=clock(),<font class=\\"jammer\\">5 u7 B/ H- B# Y5 c$ m; ?6 b1 m</font><br /> oo{k=zeros[5,5]},<font class=\\"jammer\\">% A6 Z* @1 Z5 b0 C) c</font><br /> i=0,((i++)<100000).while{<br /> <span style=\\"display:none\\">, d1 v- r; {5 c- ]& f4 j! N</span>   oo{<font class=\\"jammer\\">& K5 Q5 m. [3 P9 k' |</font><br />     a=rand[5,7], b=rand[7,5],<br /> <span style=\\"display:none\\">/ T3 _* y5 A# N' e9 F- m# J</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\\">: R: v8 j' `2 b\\" U1 I</span>   }<font class=\\"jammer\\">9 g5 C5 a\\" {7 _, R  T! u2 t</font><br /> },<font class=\\"jammer\\">+ E& c/ A3 m3 l  X! e9 T1 W6 K' v</font><br /> k.outm(),<font class=\\"jammer\\">2 n7 c( G! I0 e. A\\" _# K  i</font><br /> [clock()-t0]/1000;

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

6 T5 ^, W, j/ _' ~9 ^6 e5 _Forcal代码2：比较好看些的代码，似乎也比matlab稍快吧？

2. !using["math","sys"];
   
3. 9 [6 R& S0 V; q& o
4. (:t0,k,i,a,b)=
   
5.   l6 M& D9 K8 k* T: K, H6 a( h
6. {
   
7. 0 ~* j7 O2 E9 V' Y$ L9 U' s
8.   t0=clock(),
   
9. & w( P' @! ~6 k3 K
10.   oo{k=zeros[5,5]},
    
11. 9 @6 s\\" o: y' b$ e% u& d
12.   i=0,((i++)<100000).while{
    
13. : B- i, z1 e7 J
14.     oo{: u3 ~; M6 y3 @6 `! {' C
    
16.       a=rand[5,7], b=rand[7,5],
    
17. & a* e8 Z- x! c6 W5 K1 c
18.       k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)-a(neg:6)*b(3:neg)
    
19. ! `( v; C* w# _, o7 P# u' p. U
20.     }
    
21. 5 W/ X/ b4 }\\" P) F' H% U
22.   },' f5 F9 l% e. Q& _. _& ^
    
24.   k.outm(),
    
25. 1 t) k$ A  R- ?; f7 P8 v
26.   [clock()-t0]/1000
    
27. ) M# |$ U% a* M, R1 o. V
28. };

!using["math","sys"];<br /> <span style=\\"display:none\\">9 [6 R& S0 V; q& o</span> (:t0,k,i,a,b)=<br /> <span style=\\"display:none\\">  l6 M& D9 K8 k* T: K, H6 a( h</span> {<br /> <span style=\\"display:none\\">0 ~* j7 O2 E9 V' Y$ L9 U' s</span>   t0=clock(),<br /> <span style=\\"display:none\\">& w( P' @! ~6 k3 K</span>   oo{k=zeros[5,5]},<br /> <span style=\\"display:none\\">9 @6 s\\" o: y' b$ e% u& d</span>   i=0,((i++)<100000).while{<br /> <span style=\\"display:none\\">: B- i, z1 e7 J</span>     oo{<font class=\\"jammer\\">: u3 ~; M6 y3 @6 `! {' C</font><br />       a=rand[5,7], b=rand[7,5],<br /> <span style=\\"display:none\\">& a* e8 Z- x! c6 W5 K1 c</span>       k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)-a(neg:6)*b(3:neg)<br /> <span style=\\"display:none\\">! `( v; C* w# _, o7 P# u' p. U</span>     }<br /> <span style=\\"display:none\\">5 W/ X/ b4 }\\" P) F' H% U</span>   },<font class=\\"jammer\\">' f5 F9 l% e. Q& _. _& ^</font><br />   k.outm(),<br /> <span style=\\"display:none\\">1 t) k$ A  R- ?; f7 P8 v</span>   [clock()-t0]/1000<br /> <span style=\\"display:none\\">) M# |$ U% a* M, R1 o. V</span> };

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

例子2：
一段程序的Forcal实现：
6 `: m0 y8 `9 x

1. //用C++代码描述为：
   * |# q4 J. s1 Z4 P* z2 |% O
2. s=0.0;  ! G3 @0 u\" F\" B
   
3. for(x=0.0;x<=1.0;x=x+0.0011)  5 a. L\" m  \; d, W5 }
   
4. {
   9 I8 t9 y8 r1 ^9 R2 I' \
5.   for(y=1.0;y<=2.0;y=y+0.0009)
     t/ x% y0 a6 A* ]
6.   {5 J\" k$ l1 {9 q1 u5 _
   
7.   s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   : S5 W9 P* U( }9 S7 G
8.   }8 I! l# s7 y4 {9 Z0 ]1 W, h- z+ ^
   
9. }  

复制代码

结果：
$ W) m/ R% ?8 X. N/ |1008606.64947441
* v9 x$ o4 A2 l! |0.609 //时间

2 d# j) O, L5 J: K9 F/ S这个matlab程序段是网友yycs001给出的。
0 P; E* D! a0 s3 P& t/ Y1 L

1. %file speedtest.m
   8 @& o+ c: `# P2 @* P9 G
2. function speedtest
   ! N4 v3 T/ _! }. K\" a
3. format long
   ! y  |) \% [8 o% F
4. tic8 g& {4 W0 D6 h( T
   
5. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);2 A( V2 B4 a/ V5 N' v3 O7 S
   
6. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))
   9 ^) O# }9 W1 g4 E. P% H* K  w
7. toc

复制代码

( o% U5 d1 ^: t( eForcal代码1：**数组求和函数Sum，完全矢量化的代码
& Z5 p9 p/ ?7 e* T! j; J9 \& M/ g

2. !using["math","sys"];3 g+ L' d$ w- q1 i
   
4. mvar:
   
5. ! e; s. T8 U$ X* H\\" a% `, A) A
6. t=clock(),
   
7. 7 `3 f2 D6 n+ ~3 T+ f3 A; B
8. oo{
   
9. 6 m( L/ |; h( C\\" [; N5 h' l
10.   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
11. 9 U! i2 a& o3 y& V0 M# z% k' f
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]$ B3 r3 V/ p) T0 N% E
    
14. };
    
15. ' t$ _, n! C2 R) v6 {
16. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">3 g+ L' d$ w- q1 i</font><br /> mvar:<br /> <span style=\\"display:none\\">! e; s. T8 U$ X* H\\" a% `, A) A</span> t=clock(),<br /> <span style=\\"display:none\\">7 `3 f2 D6 n+ ~3 T+ f3 A; B</span> oo{<br /> <span style=\\"display:none\\">6 m( L/ |; h( C\\" [; N5 h' l</span>   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">9 U! i2 a& o3 y& V0 M# z% k' f</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]<font class=\\"jammer\\">$ B3 r3 V/ p) T0 N% E</font><br /> };<br /> <span style=\\"display:none\\">' t$ _, n! C2 R) v6 {</span> [clock()-t]/1000;

结果：
" w* ]0 R$ L% d+ k1008606.64947441
3 u! N, w" \1 q) u2 w9 o0.625 //时间

或者这个，与上面效率差别不大：

2. !using["math","sys"];
   
3. & }+ j; H. ~' P/ J( B4 M7 L
4. mvar:
   
5. % }\\" f0 {. y3 X* }5 |5 F, g
6. t=clock(),7 z4 N2 o& Y\\" F/ C0 m+ Q$ {
   
8. oo{! v3 }/ u. q! `
   
10.   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
11. % e, ^1 U* i. U* i& Y/ f/ c
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))))))]]
    
13. 1 ]) y1 R\\" a. O6 \+ a
14. };- k9 K\\" ~7 n: d6 r/ f6 T
    
16. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">& }+ j; H. ~' P/ J( B4 M7 L</span> mvar:<br /> <span style=\\"display:none\\">% }\\" f0 {. y3 X* }5 |5 F, g</span> t=clock(),<font class=\\"jammer\\">7 z4 N2 o& Y\\" F/ C0 m+ Q$ {</font><br /> oo{<font class=\\"jammer\\">! v3 }/ u. q! `</font><br />   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">% e, ^1 U* i. U* i& Y/ f/ c</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))))))]]<br /> <span style=\\"display:none\\">1 ]) y1 R\\" a. O6 \+ a</span> };<font class=\\"jammer\\">- k9 K\\" ~7 n: d6 r/ f6 T</font><br /> [clock()-t]/1000;

Forcal代码2：求和函数sum，非矢量化代码
$ f# K  a' l% Z4 V

1. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   ; z; @  E; w' Y$ R* |0 W
2. sum["f",0,1,0.0011,1,2,0.0009];

复制代码

结果：
1008606.64947441
0.719 //时间

Forcal代码3：while循环
. V+ _2 a  [, B1 i" ^  J! I% A* I: ]

1. mvar:! h& h* F) l! A, e2 @2 W1 C& Y4 x
   
2. t=sys::clock();- R- y! b) X+ d2 I& l
   
3. s=0,x=0,
   ( q5 o8 j4 c& p* S# t/ \) g
4. while{x<=1,  //while循环算法； ! X6 R\" T$ Z6 p7 L( ?9 B- J
   
5.    y=1, % b) F/ _. s- |* a/ e. r& O
   
6.    while{y<=2,
   , ?- u; t9 r2 E
7.        s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))), : `  X& P4 t% ]\" f% {; A
   
8.        y=y+0.0009
   0 {+ i0 Q/ i4 _* l\" t- B7 x) [
9.       }, - x+ e; M8 X3 [$ C. R
   
10.    x=x+0.0011
    9 v! F/ p  y  O) q! b. [
11. }, * P7 A1 i0 @) y# F
    
12. s;
    ! m( ?4 z3 S6 }3 `5 O$ K
13. [sys::clock()-t]/1000;

复制代码

结果：
4 }+ e+ B. x! R: ^- u' d# q1008606.64947441
0.734 //时间
5 W5 f! z+ z& d0 r% A8 `
大家可下载OpenFC进行测试：http://www.forcal.net/xiazai/forcal9/openfc32w.rar

0 M" ^' {( f% t, N" l0 d注意Forcal的矢量化代码第一次运行有时效率较低。

例子1中Forcal和matlab都是矢量化代码，但matlab跑不过Forcal。该例子的特点是函数调用频繁，临时变量生成多，但矩阵很小，矩阵的各种函数运行时耗时较少。故说明Forcal函数调用＋变量管理效率优于matlab。

例子2中Forcal的矢量化代码是最快的，但与matlab的矢量化代码相比仍有差距。该例子的特点是函数调用少，临时变量也少，但矩阵大。故说明Forcal的各种矩阵函数Sin、Cos及矩阵的加减运算等函数的内部设计不及matlab。

如能在函数内部设计上下点功夫，例子2超越matlab也是可能的。在这方面，期待高手们的指点。
( H4 N( [" U) H. `* b8 Q% W, I
- E6 e; l/ b% m# D如果例子2速度也超越了matlab ，matlab矢量化的神秘面纱就揭开了。

顺便说一下，例子1如果用C++的运算符重载来实现，速度将比Forcal慢一些，也就是说，在涉及运算符重载时，脚本的效率有时比C++还要高些。
$ |! f. v- R8 b1 _$ [) @4 ^

master_math

,,,顶顶更健康！

wznzy0822

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

zhwqqiangge

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

qbist

我 正在  学  Matlab  软件！~~