关于matlab代码矢量化的理解

forcal

代码矢量化是matlab的精髓，其基本特点是运行速度快和代码简洁，它是如何实现的？
0 W* J/ i. ~$ a! s- N
7 [: h' t- q& x$ t; ~1 U按我的理解，代码矢量化的本质就是设计专门的函数对数组元素集中运算，这样可提高运行速度，同时兼有代码简洁的特点。
2 W5 q6 N; Q1 l% d) @# ~, H6 A
' L4 G, Q! w; G& z0 l6 F0 Q对matlab的理解比较肤浅，但也确实看不出有更深意义的东西，望解惑。

' f; J4 X/ V5 P: Z2 ^& G大家有什么看法，愿畅所欲言。

qbist

我 正在  学  Matlab  软件！~~

zhwqqiangge

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

wznzy0822

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

master_math

,,,顶顶更健康！

forcal

* o. W  r$ _# i& |0 {
我正在练手设计的FcMath库也打算以矩阵运算为基础，设计一些专门的函数对数组元素集中运算，运行效率确实有所提高（甚至有些涉及矩阵的算法比matlab还快），代码也简洁了，但不知这是不是矢量化？
# ~/ s- ~6 c9 i( j& W) H  A
脚本运行效率应该取决于函数调度效率、对象管理效率和函数内部算法的实现。

4 I, ]5 R" V, P; B我感觉，matlab的函数调度效率较低，对象管理效率这个不好说，但一些函数内部的设计比较优秀。故有些Forcal代码比matlab快，而有些慢。

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

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

1. clear all
   $ _, X  r! f# ?0 z. v4 F
2. clc) `7 ^& i3 V; A* t& z$ e4 @& n
   
3. tic
   : M9 `. A7 @8 R# X! [( H+ O/ a
4. k = zeros(5,5); % //生成5×5全0矩阵
   $ c, X; p) |1 {6 j5 _- J) v
5. % 循环计算以下程序段1000 00次：
   \" p( y- A& O( ~\" d, v\" ~- ?- t% z5 g' _
6. for m = 1:1000 00. u# K! w8 p9 ]/ C, E& S7 H
   
7.     a = rand(5,7);
   . k- K8 X' ?! v. D. E\" u
8.     b = rand(7,5);%//生成5×7矩阵a,7×5矩阵b,用0～1之间的随机数初始化
   3 @' |  K4 d6 c! A6 O! k
9.     k = k + a * b + a(1:5, 2:6) * b(2:6, 1:5) - a(:, 7) * b(3, :);
   ! x8 `5 q$ `: v/ V1 X! L5 B$ H/ U
10. end
    * d+ G, g/ K) \: }5 h' |
11. k( h8 M) g  Y, F/ g$ X  M. Q
    
12. toc

复制代码

- ], V! s5 T: Y: u! EForcal代码1：运行稍快的代码，比matlab约快10％吧？
# q" D- G0 i) |6 y  O

2. !using["math","sys"];) O+ l, Z- S2 P$ ~
   
4. mvar:5 A# U9 Q; v1 \) F1 V) {
   
6. t0=clock(),- d& L, [1 H3 L+ p
   
8. oo{k=zeros[5,5]},
   
9. 7 [* L- Q% ~; U7 V\\" z# F/ C7 X
10. i=0,((i++)<100000).while{\\" a' C8 ?- R6 Q& ^5 J7 Q$ d
    
12.   oo{
    
13. ) R4 r+ u% S2 b% y& ^* l
14.     a=rand[5,7], b=rand[7,5],
    
15. 7 |7 m, ^. Z! D* @: P+ D
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)]4 m) E/ j# o! _8 Y5 Z- q% X3 O- j
    
18.   }8 |: w& x% ~. Y* n2 |: @0 h+ q
    
20. },
    
21. ' C/ O: ?6 E) K. S- t4 z1 I
22. k.outm(),
    
23. ( s% L3 e. r  I6 T+ V5 X
24. [clock()-t0]/1000;

!using["math","sys"];<font class=\\"jammer\\">) O+ l, Z- S2 P$ ~</font><br /> mvar:<font class=\\"jammer\\">5 A# U9 Q; v1 \) F1 V) {</font><br /> t0=clock(),<font class=\\"jammer\\">- d& L, [1 H3 L+ p</font><br /> oo{k=zeros[5,5]},<br /> <span style=\\"display:none\\">7 [* L- Q% ~; U7 V\\" z# F/ C7 X</span> i=0,((i++)<100000).while{<font class=\\"jammer\\">\\" a' C8 ?- R6 Q& ^5 J7 Q$ d</font><br />   oo{<br /> <span style=\\"display:none\\">) R4 r+ u% S2 b% y& ^* l</span>     a=rand[5,7], b=rand[7,5],<br /> <span style=\\"display:none\\">7 |7 m, ^. Z! D* @: P+ D</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)]<font class=\\"jammer\\">4 m) E/ j# o! _8 Y5 Z- q% X3 O- j</font><br />   }<font class=\\"jammer\\">8 |: w& x% ~. Y* n2 |: @0 h+ q</font><br /> },<br /> <span style=\\"display:none\\">' C/ O: ?6 E) K. S- t4 z1 I</span> k.outm(),<br /> <span style=\\"display:none\\">( s% L3 e. r  I6 T+ V5 X</span> [clock()-t0]/1000;

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

Forcal代码2：比较好看些的代码，似乎也比matlab稍快吧？
) b6 C# s1 {. v) ~4 y2 r

2. !using["math","sys"];+ L8 {  G' f/ a6 w  p
   
4. (:t0,k,i,a,b)=
   
5. 4 Q$ q& I* _, ~% Z! s* N( p; j
6. {# ?6 E3 v\\" i' b5 l& D7 P% M( w
   
8.   t0=clock(),3 I9 B& o4 ~- _& J& x  ~2 w, U
   
10.   oo{k=zeros[5,5]},
    
11. . J* u5 N\\" @$ j7 m3 n
12.   i=0,((i++)<100000).while{
    
13. 7 b. ^* R# _; V3 Y3 e  r0 h, d
14.     oo{6 }* E0 J- I# t; Z
    
16.       a=rand[5,7], b=rand[7,5],
    
17. 3 C+ p& \4 c6 _  [3 {
18.       k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)-a(neg:6)*b(3:neg)' a' |; Z6 j7 k) t( ^! D& E% u
    
20.     }6 M* }6 ]# {+ P  L/ p/ w
    
22.   },! _; F/ s* Q3 g% N$ T; w# [$ k
    
24.   k.outm(),
    
25. 6 G+ X/ d% @8 }
26.   [clock()-t0]/1000; E; e$ |! ]2 n7 p/ Y( `8 H6 {( P  Y
    
28. };

!using["math","sys"];<font class=\\"jammer\\">+ L8 {  G' f/ a6 w  p</font><br /> (:t0,k,i,a,b)=<br /> <span style=\\"display:none\\">4 Q$ q& I* _, ~% Z! s* N( p; j</span> {<font class=\\"jammer\\"># ?6 E3 v\\" i' b5 l& D7 P% M( w</font><br />   t0=clock(),<font class=\\"jammer\\">3 I9 B& o4 ~- _& J& x  ~2 w, U</font><br />   oo{k=zeros[5,5]},<br /> <span style=\\"display:none\\">. J* u5 N\\" @$ j7 m3 n</span>   i=0,((i++)<100000).while{<br /> <span style=\\"display:none\\">7 b. ^* R# _; V3 Y3 e  r0 h, d</span>     oo{<font class=\\"jammer\\">6 }* E0 J- I# t; Z</font><br />       a=rand[5,7], b=rand[7,5],<br /> <span style=\\"display:none\\">3 C+ p& \4 c6 _  [3 {</span>       k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)-a(neg:6)*b(3:neg)<font class=\\"jammer\\">' a' |; Z6 j7 k) t( ^! D& E% u</font><br />     }<font class=\\"jammer\\">6 M* }6 ]# {+ P  L/ p/ w</font><br />   },<font class=\\"jammer\\">! _; F/ s* Q3 g% N$ T; w# [$ k</font><br />   k.outm(),<br /> <span style=\\"display:none\\">6 G+ X/ d% @8 }</span>   [clock()-t0]/1000<font class=\\"jammer\\">; E; e$ |! ]2 n7 p/ Y( `8 H6 {( P  Y</font><br /> };

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

" v8 q5 z/ r+ c; u+ h# a1 T- P例子2：
一段程序的Forcal实现：
* t% i1 t( }: v6 B, ~( c& z( B

1. //用C++代码描述为：
   5 s% q& j. Z, i1 H\" ^/ O3 W- _5 s
2. s=0.0;  
   + Z9 a( @% T* m7 s; u+ P/ m
3. for(x=0.0;x<=1.0;x=x+0.0011)  
   9 Y0 Z5 q2 t1 o, S+ d
4. {) S, j' S) n, p0 U6 @* t
   
5.   for(y=1.0;y<=2.0;y=y+0.0009)9 p5 A2 a3 f. ]7 R: J
   
6.   {% b( N# O( I2 R
   
7.   s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));( C1 t/ J3 p: S- b5 C- c
   
8.   }2 q( G9 X0 @/ [% d& z
   
9. }  

复制代码

结果：
9 N, q* l6 K# l/ ^1008606.64947441
# c# m9 c: o1 O7 }" X+ U$ M! R0.609 //时间
% Q  S' X5 v' C
这个matlab程序段是网友yycs001给出的。
) P" W' q% A5 ?, S

1. %file speedtest.m8 O! r% A& ]6 |6 W3 |9 k  ]
   
2. function speedtest/ D7 O1 W% m2 L3 ~\" j' e  e3 ]
   
3. format long
   ; D7 a7 D3 u, T( A7 ~: V- W) Q1 h# V4 \
4. tic
   0 o, R6 N5 J/ F
5. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);0 S+ I1 Y6 z+ R\" c8 v
   
6. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))
   \" O0 f- {  z0 x
7. toc

复制代码

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

2. !using["math","sys"];
   
3. ! J* b$ g! Z. U9 D& B5 N. H/ ]' I
4. mvar:$ S, h. T& p- P* R\\" o' A4 ]3 f& d8 a
   
6. t=clock(),
   
7. \\" Q. }% U/ O# j: D
8. oo{/ h- t9 N\\" T& C
   
10.   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],. c% E+ U, g% ~2 Y# T. B
    
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. / Y\\" G; f2 N; _2 y/ v
14. };
    
15.   Q& p8 C; u\\" {9 _0 ^; j$ t7 ^
16. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">! J* b$ g! Z. U9 D& B5 N. H/ ]' I</span> mvar:<font class=\\"jammer\\">$ S, h. T& p- P* R\\" o' A4 ]3 f& d8 a</font><br /> t=clock(),<br /> <span style=\\"display:none\\">\\" Q. }% U/ O# j: D</span> oo{<font class=\\"jammer\\">/ h- t9 N\\" T& C</font><br />   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">. c% E+ U, g% ~2 Y# T. B</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\\">/ Y\\" G; f2 N; _2 y/ v</span> };<br /> <span style=\\"display:none\\">  Q& p8 C; u\\" {9 _0 ^; j$ t7 ^</span> [clock()-t]/1000;

结果：
1008606.64947441
0.625 //时间
, E4 t5 P9 O  g1 h: p# R' Z5 i
% h# p% b! k. h3 b或者这个，与上面效率差别不大：
4 i& A: s" F7 I" v$ K% {& U

2. !using["math","sys"];4 w1 Y9 X, n/ _/ _/ j$ F  m+ @& y! ^
   
4. mvar:/ A. y/ u5 d; C- Q( I
   
6. t=clock(),
   
7. 2 r+ a4 m0 f0 h! W* I: T
8. oo{
   
9. 4 F2 J- \: g; V. u' g7 f
10.   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],4 Z! }- j8 f$ K0 W; D% W
    
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. # e, _$ z4 j% `8 v$ X2 H/ ^
14. };( E( B# l3 Q- ~9 _
    
16. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">4 w1 Y9 X, n/ _/ _/ j$ F  m+ @& y! ^</font><br /> mvar:<font class=\\"jammer\\">/ A. y/ u5 d; C- Q( I</font><br /> t=clock(),<br /> <span style=\\"display:none\\">2 r+ a4 m0 f0 h! W* I: T</span> oo{<br /> <span style=\\"display:none\\">4 F2 J- \: g; V. u' g7 f</span>   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">4 Z! }- j8 f$ K0 W; D% W</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\\"># e, _$ z4 j% `8 v$ X2 H/ ^</span> };<font class=\\"jammer\\">( E( B# l3 Q- ~9 _</font><br /> [clock()-t]/1000;

' B0 \7 s# u. ?4 RForcal代码2：求和函数sum，非矢量化代码
! l8 [3 P- o, H0 d' A/ p" q

1. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); ( T: h( W9 F4 ~3 a& A
   
2. sum["f",0,1,0.0011,1,2,0.0009];

复制代码

结果：
' j! ?, n2 ^! P6 \( a1008606.64947441
0.719 //时间
8 h0 D1 R& H/ E) @. L' s
+ T' t5 q- ]- QForcal代码3：while循环
* X- d' N1 D/ n3 s  d

1. mvar:
   & w/ P  [: b$ ]$ m5 q5 ~7 k% q% J
2. t=sys::clock();
   : ?$ \0 x' D6 k$ d* ?# y
3. s=0,x=0,
   \" j' D, J2 j, Y- n6 U0 h- z
4. while{x<=1,  //while循环算法； / _\" ^3 A0 v7 y/ k8 P# W
   
5.    y=1,   I: Y) H3 R4 F; P$ L* M4 C
   
6.    while{y<=2,
   4 [6 I' o0 s\" S\" J* u
7.        s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
   3 L- f) \6 k) F, _# m
8.        y=y+0.0009 5 C* G0 P6 ?( T+ b
   
9.       },
   ' H  J8 V: V4 _1 t. F- |  ]
10.    x=x+0.0011
    9 I& v  i6 V+ M3 \$ C
11. }, 5 t6 u! B, d# Z- [3 `% I3 H
    
12. s;: F! y. g+ Z4 B! |; |0 a
    
13. [sys::clock()-t]/1000;

复制代码

结果：
! K/ v! S8 u, p( R) y" A1008606.64947441
3 u8 J8 @. N/ h8 @# p0.734 //时间

$ c  z& K1 C4 @/ @$ j大家可下载OpenFC进行测试：http://www.forcal.net/xiazai/forcal9/openfc32w.rar
: Q2 {& d- _/ n9 M# U7 `
注意Forcal的矢量化代码第一次运行有时效率较低。
: ]- s( u( l, D3 U$ C/ D
. @* I- V9 S% A: C, E例子1中Forcal和matlab都是矢量化代码，但matlab跑不过Forcal。该例子的特点是函数调用频繁，临时变量生成多，但矩阵很小，矩阵的各种函数运行时耗时较少。故说明Forcal函数调用＋变量管理效率优于matlab。

- {9 @( [* O2 p4 D3 L; l) Q例子2中Forcal的矢量化代码是最快的，但与matlab的矢量化代码相比仍有差距。该例子的特点是函数调用少，临时变量也少，但矩阵大。故说明Forcal的各种矩阵函数Sin、Cos及矩阵的加减运算等函数的内部设计不及matlab。
( e" }7 x: J% U* X2 b' k# N8 U
如能在函数内部设计上下点功夫，例子2超越matlab也是可能的。在这方面，期待高手们的指点。

0 X' S& G$ q- G如果例子2速度也超越了matlab ，matlab矢量化的神秘面纱就揭开了。
% G* V0 X$ }* F" c8 u0 \
, R* n/ ~. i" K& C顺便说一下，例子1如果用C++的运算符重载来实现，速度将比Forcal慢一些，也就是说，在涉及运算符重载时，脚本的效率有时比C++还要高些。

forcal

讨论有益！以下是在其他论坛的讨论帖子：
' U, [8 U% H) p5 f4 k) U% osimwe：http://forum.simwe.com/thread-952532-1-3.html
csdn：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

& D5 H5 o6 G0 W我在多个帖子中有不同说明，但将这些说明再集中到一个帖子中比较麻烦，大家相互参考一下，看能否把这个问题解决了。

forcal

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

; B7 N* Z! r$ a+ ^关于最快速的矩阵乘实现的讨论。
4 D; m+ {! I* w6 s/ U- r& {