关于matlab代码矢量化的理解

forcal

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

( J" x3 S. c" J) K3 v1 U! M按我的理解，代码矢量化的本质就是设计专门的函数对数组元素集中运算，这样可提高运行速度，同时兼有代码简洁的特点。

: p6 C- h0 N/ \: {% x8 t对matlab的理解比较肤浅，但也确实看不出有更深意义的东西，望解惑。
' |" P1 o( Z* H- `
大家有什么看法，愿畅所欲言。
; y/ J) \* u/ W

qbist

我 正在  学  Matlab  软件！~~

zhwqqiangge

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

wznzy0822

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

master_math

,,,顶顶更健康！

forcal

我正在练手设计的FcMath库也打算以矩阵运算为基础，设计一些专门的函数对数组元素集中运算，运行效率确实有所提高（甚至有些涉及矩阵的算法比matlab还快），代码也简洁了，但不知这是不是矢量化？
$ V2 B( C# S/ l7 I/ e' d
& ~9 H# I) p( S# T2 ~脚本运行效率应该取决于函数调度效率、对象管理效率和函数内部算法的实现。
3 h- @5 Z6 {2 S3 j, o; ]6 [
- N4 l# r+ f+ t$ G9 ~我感觉，matlab的函数调度效率较低，对象管理效率这个不好说，但一些函数内部的设计比较优秀。故有些Forcal代码比matlab快，而有些慢。
  ]4 `# Y4 V1 a7 Y% ]6 Q* A5 g" S
5 U. d' H8 _( x  D4 I4 P: O以下例子体现了Forcal和matlab的效率差别所在。

# y1 R- n: E* r" }# M- [这个matlab程序段是网友lin2009 给出的，理论结果是每个元素均为275000。

1. clear all
   5 _6 \1 f4 L: X2 w! B6 n+ F4 m
2. clc) h3 v$ Z3 t6 i) B$ w
   
3. tic/ y7 C6 z' F( ^0 T1 \; J
   
4. k = zeros(5,5); % //生成5×5全0矩阵  i0 T- P: t$ w* v% \
   
5. % 循环计算以下程序段1000 00次：9 L6 z2 ]1 y/ m$ o9 x
   
6. for m = 1:1000 006 e, Q% R- x3 q1 T8 u8 D$ @
   
7.     a = rand(5,7);
   ! |, d* w0 G! ?6 V* D\" [+ @7 X+ _
8.     b = rand(7,5);%//生成5×7矩阵a,7×5矩阵b,用0～1之间的随机数初始化4 n2 R* n2 {1 d, D+ Z
   
9.     k = k + a * b + a(1:5, 2:6) * b(2:6, 1:5) - a(:, 7) * b(3, :);
   ) Y: H! H. h1 n5 }2 o
10. end
    , [& o' n0 C' c1 q
11. k% P7 I) z! I3 \, K7 C
    
12. toc

复制代码

$ u8 p$ f' I' N9 U* \$ L" e6 DForcal代码1：运行稍快的代码，比matlab约快10％吧？

2. !using["math","sys"];: Z& E% l0 |. r3 g. B  s; B
   
4. mvar:
   
5. / t1 P3 t. y8 Y
6. t0=clock(),4 `) y9 m- O. }+ p) X# z
   
8. oo{k=zeros[5,5]},
   
9. , r\\" _' `7 `9 `* P2 |1 e5 C( G( s
10. i=0,((i++)<100000).while{
    
11. . m4 F( L6 u! B1 N7 L+ p
12.   oo{2 w! `& h$ L- {3 y* \
    
14.     a=rand[5,7], b=rand[7,5],
    
15. \\" s3 F. Y9 z/ N8 Z$ E5 S+ ~
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. 5 B3 K0 E: {- t5 ^8 P
18.   }. O- E* [- i2 ]0 L9 N* t$ X
    
20. },
    
21. - o% o$ J. \8 D# f. g+ R
22. k.outm(),
    
23. & P# P: |3 ]5 V6 |; W/ a! O2 I
24. [clock()-t0]/1000;

!using["math","sys"];<font class=\\"jammer\\">: Z& E% l0 |. r3 g. B  s; B</font><br /> mvar:<br /> <span style=\\"display:none\\">/ t1 P3 t. y8 Y</span> t0=clock(),<font class=\\"jammer\\">4 `) y9 m- O. }+ p) X# z</font><br /> oo{k=zeros[5,5]},<br /> <span style=\\"display:none\\">, r\\" _' `7 `9 `* P2 |1 e5 C( G( s</span> i=0,((i++)<100000).while{<br /> <span style=\\"display:none\\">. m4 F( L6 u! B1 N7 L+ p</span>   oo{<font class=\\"jammer\\">2 w! `& h$ L- {3 y* \</font><br />     a=rand[5,7], b=rand[7,5],<br /> <span style=\\"display:none\\">\\" s3 F. Y9 z/ N8 Z$ E5 S+ ~</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\\">5 B3 K0 E: {- t5 ^8 P</span>   }<font class=\\"jammer\\">. O- E* [- i2 ]0 L9 N* t$ X</font><br /> },<br /> <span style=\\"display:none\\">- o% o$ J. \8 D# f. g+ R</span> k.outm(),<br /> <span style=\\"display:none\\">& P# P: |3 ]5 V6 |; W/ a! O2 I</span> [clock()-t0]/1000;

在我的电脑上运行时间为3.344秒。
" y, C9 I0 H! I4 n4 K, |- g
  @, J, S  B- n6 }4 }Forcal代码2：比较好看些的代码，似乎也比matlab稍快吧？

2. !using["math","sys"];1 F5 D) h# y# }8 s: v1 E+ F* Z
   
4. (:t0,k,i,a,b)=& @- W& I# }1 T
   
6. {' i5 M) W% I9 C' C
   
8.   t0=clock(),
   
9. ) U' \6 D) m* a7 K
10.   oo{k=zeros[5,5]},
    
11. $ ~5 U! D& [8 o! d' R
12.   i=0,((i++)<100000).while{
    
13. , p8 {9 Y: _3 \, Q/ L
14.     oo{
    
15. 2 ~: ~/ R# D9 [1 X6 ]
16.       a=rand[5,7], b=rand[7,5],
    
17. 1 I- n1 X* d* _0 d; ]
18.       k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)-a(neg:6)*b(3:neg)
    
19. $ C0 a  t* n% i: e
20.     }7 s, c. I$ V. J6 F( X) y8 `
    
22.   },
    
23. 4 R8 M' M* }) P- }, ^2 Z* }. s
24.   k.outm(),9 L2 k5 y. V& j7 d
    
26.   [clock()-t0]/1000
    
27. + `) h/ x  e$ g2 l' T% Q: J
28. };

!using["math","sys"];<font class=\\"jammer\\">1 F5 D) h# y# }8 s: v1 E+ F* Z</font><br /> (:t0,k,i,a,b)=<font class=\\"jammer\\">& @- W& I# }1 T</font><br /> {<font class=\\"jammer\\">' i5 M) W% I9 C' C</font><br />   t0=clock(),<br /> <span style=\\"display:none\\">) U' \6 D) m* a7 K</span>   oo{k=zeros[5,5]},<br /> <span style=\\"display:none\\">$ ~5 U! D& [8 o! d' R</span>   i=0,((i++)<100000).while{<br /> <span style=\\"display:none\\">, p8 {9 Y: _3 \, Q/ L</span>     oo{<br /> <span style=\\"display:none\\">2 ~: ~/ R# D9 [1 X6 ]</span>       a=rand[5,7], b=rand[7,5],<br /> <span style=\\"display:none\\">1 I- n1 X* d* _0 d; ]</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\\">$ C0 a  t* n% i: e</span>     }<font class=\\"jammer\\">7 s, c. I$ V. J6 F( X) y8 `</font><br />   },<br /> <span style=\\"display:none\\">4 R8 M' M* }) P- }, ^2 Z* }. s</span>   k.outm(),<font class=\\"jammer\\">9 L2 k5 y. V& j7 d</font><br />   [clock()-t0]/1000<br /> <span style=\\"display:none\\">+ `) h/ x  e$ g2 l' T% Q: J</span> };

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

例子2：
0 ?2 y$ ]! t8 F) D) K- |* D  N一段程序的Forcal实现：
+ i, R# o* k" I3 {& k. Z% x( E

1. //用C++代码描述为：7 e9 H  O- W# _& L8 t  N
   
2. s=0.0;  \" D, H; {# u5 v' c, {4 F+ f
   
3. for(x=0.0;x<=1.0;x=x+0.0011)  , E4 U3 ~1 O2 \, m+ _! ]' c* f0 z
   
4. {' t3 V9 t8 _' a5 t3 H
   
5.   for(y=1.0;y<=2.0;y=y+0.0009)' m# z6 S5 ^. R1 o
   
6.   {
   & `3 g$ L3 K  X\" m
7.   s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));! E( T0 s( \4 r0 d) n
   
8.   }4 W, m1 l  z0 ^* Q
   
9. }  

复制代码

结果：
& F$ B% X9 h# Z: A* V2 s: ~1008606.64947441
0.609 //时间

这个matlab程序段是网友yycs001给出的。

1. %file speedtest.m; g7 q, D. b) f# b\" Z
   
2. function speedtest
   & F+ G3 f! n) {2 A2 K1 o2 o
3. format long
   7 \- _/ k4 z/ l+ o# m! J( y
4. tic. V. j  {. Q- C. k4 a+ _
   
5. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);; d, b7 L& F! h1 `) i0 O
   
6. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))8 U4 D3 X& P0 w\" g4 H: @' s
   
7. toc

复制代码

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

2. !using["math","sys"];
   
3. 5 t9 P* h! {( p, g- d
4. mvar:3 b2 b* U5 W, Z; l$ B
   
6. t=clock(),
   
7. \\" }9 C* p) S4 r2 k7 C
8. oo{
   
9. 0 E7 b3 a  x' a
10.   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],! t1 x; Q* b  x, s8 z  x
    
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. + S3 a+ V  e4 m; w: s/ w  k
14. };: M$ Z  Z' g, b/ L% l
    
16. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">5 t9 P* h! {( p, g- d</span> mvar:<font class=\\"jammer\\">3 b2 b* U5 W, Z; l$ B</font><br /> t=clock(),<br /> <span style=\\"display:none\\">\\" }9 C* p) S4 r2 k7 C</span> oo{<br /> <span style=\\"display:none\\">0 E7 b3 a  x' a</span>   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">! t1 x; Q* b  x, s8 z  x</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\\">+ S3 a+ V  e4 m; w: s/ w  k</span> };<font class=\\"jammer\\">: M$ Z  Z' g, b/ L% l</font><br /> [clock()-t]/1000;

结果：
1008606.64947441
0.625 //时间

) D9 H9 J) [" m- ~1 w; d或者这个，与上面效率差别不大：
4 S+ @$ X- Z& y9 {% y0 m" _

2. !using["math","sys"];! G% E# Q) P# Q( _3 ~
   
4. mvar:
   
5. 7 z/ r/ F# n+ _
6. t=clock(),
   
7. . l8 F4 ?4 U& R% |
8. oo{
   
9. 3 E% y* P# U. Y7 u
10.   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],  l- Y. v4 s. l3 O/ b
    
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. 8 N* X; M5 N* j- o% O: U) q
14. };
    
15. * w; e. t3 B' s. \% }, a1 ^  {* e
16. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">! G% E# Q) P# Q( _3 ~</font><br /> mvar:<br /> <span style=\\"display:none\\">7 z/ r/ F# n+ _</span> t=clock(),<br /> <span style=\\"display:none\\">. l8 F4 ?4 U& R% |</span> oo{<br /> <span style=\\"display:none\\">3 E% y* P# U. Y7 u</span>   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">  l- Y. v4 s. l3 O/ b</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))))))]]<br /> <span style=\\"display:none\\">8 N* X; M5 N* j- o% O: U) q</span> };<br /> <span style=\\"display:none\\">* w; e. t3 B' s. \% }, a1 ^  {* e</span> [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))))); ( k9 f; `* n5 ]% l
   
2. sum["f",0,1,0.0011,1,2,0.0009];

复制代码

结果：
* h; X$ `" u1 A$ D1008606.64947441
0.719 //时间

) c" J2 E5 t7 S0 H1 I. U, zForcal代码3：while循环
# b  }8 e, [" v) t% q

1. mvar:  L6 c/ z: m& |7 O
   
2. t=sys::clock();! k4 r. B1 G+ ]( [: W; \% j3 a
   
3. s=0,x=0,
   ) S) a0 Z2 E0 p! v0 u- X9 |
4. while{x<=1,  //while循环算法；
   / w6 c0 _' t! V( o\" W4 T7 p. c
5.    y=1, . n# D+ Q6 F: @2 U5 `* u# n
   
6.    while{y<=2, ) i5 z, {5 o2 q
   
7.        s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))), 3 A; K4 y# E, q0 P\" D7 M6 S
   
8.        y=y+0.0009 & T/ `: A6 a5 a9 I4 \
   
9.       }, * O/ t# q& d\" s\" t3 ~/ A2 @
   
10.    x=x+0.0011
    9 S7 [9 I1 _+ R+ F( @3 M
11. }, 8 C6 L$ v7 T\" C. r  y/ N
    
12. s;; }' M  F9 o# w
    
13. [sys::clock()-t]/1000;

复制代码

结果：
! n8 S& \8 w& y1008606.64947441
! U1 J& o$ ~2 m" T2 T  L4 T4 X6 e0.734 //时间

5 w+ }/ {5 w8 r, y" D, z# ~+ G大家可下载OpenFC进行测试：http://www.forcal.net/xiazai/forcal9/openfc32w.rar
8 T9 T' {% C( N. n+ Z
7 D& L+ d' W/ d2 }注意Forcal的矢量化代码第一次运行有时效率较低。

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

' x  G& X6 v, P: o  R9 N- P+ E例子2中Forcal的矢量化代码是最快的，但与matlab的矢量化代码相比仍有差距。该例子的特点是函数调用少，临时变量也少，但矩阵大。故说明Forcal的各种矩阵函数Sin、Cos及矩阵的加减运算等函数的内部设计不及matlab。
  |4 h% e" S) h# o" A
! }, g8 w& J% l! T& ]( ~% R如能在函数内部设计上下点功夫，例子2超越matlab也是可能的。在这方面，期待高手们的指点。

如果例子2速度也超越了matlab ，matlab矢量化的神秘面纱就揭开了。

顺便说一下，例子1如果用C++的运算符重载来实现，速度将比Forcal慢一些，也就是说，在涉及运算符重载时，脚本的效率有时比C++还要高些。
7 g% T  G9 Z! N0 E; |& @# t7 I

forcal

讨论有益！以下是在其他论坛的讨论帖子：
simwe：http://forum.simwe.com/thread-952532-1-3.html
1 ^' y5 [1 z6 c- i- Q9 [& ]! y8 acsdn：http://topic.csdn.net/u/20101006/21/2ed9e5c1-cc9f-4623-b1c0-ebd5d1b5a98a.html
cadn：http://topic.csdn.net/u/20101010/15/3bcf2fe0-0abd-4c29-b854-b1d007b16863.html

forcal

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

$ w# @. ^2 {# D! @2 i7 g. D9 ^我在多个帖子中有不同说明，但将这些说明再集中到一个帖子中比较麻烦，大家相互参考一下，看能否把这个问题解决了。

forcal

参考：http://bbs.emath.ac.cn/thread-2727-1-1.html
, C( ]3 W" `2 E- i( Y/ i& }' Q. P
# [: W% `# [: C, D4 R7 P' R" W关于最快速的矩阵乘实现的讨论。
3 t  z( T8 A. l& {; m% w