关于matlab代码矢量化的理解

forcal

代码矢量化是matlab的精髓，其基本特点是运行速度快和代码简洁，它是如何实现的？
  I) e# {) @8 V3 x$ {7 x7 x: i+ c" C; p
按我的理解，代码矢量化的本质就是设计专门的函数对数组元素集中运算，这样可提高运行速度，同时兼有代码简洁的特点。

对matlab的理解比较肤浅，但也确实看不出有更深意义的东西，望解惑。
, g3 l& w" T! }
' W! |5 Q/ l) F9 }. h2 m8 h大家有什么看法，愿畅所欲言。
/ I7 O$ @  g* O1 ?, V" f! c3 m3 B. O" x

forcal

参考：http://bbs.emath.ac.cn/thread-2727-1-1.html
1 U6 k7 u3 T. b
关于最快速的矩阵乘实现的讨论。
' S/ _* }& R0 o3 P* O

forcal

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

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

forcal

讨论有益！以下是在其他论坛的讨论帖子：
simwe：http://forum.simwe.com/thread-952532-1-3.html
3 ?+ r" P0 g! ^  E/ n% ycsdn：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

. q9 M2 _$ j- F4 |+ w, Q1 w4 T我正在练手设计的FcMath库也打算以矩阵运算为基础，设计一些专门的函数对数组元素集中运算，运行效率确实有所提高（甚至有些涉及矩阵的算法比matlab还快），代码也简洁了，但不知这是不是矢量化？
8 T- |5 @6 `8 b, O1 C2 C# b0 ~2 Y. F
8 ^- Y1 q% T8 |9 W0 V' X6 T脚本运行效率应该取决于函数调度效率、对象管理效率和函数内部算法的实现。
* t. z5 e/ A: J
我感觉，matlab的函数调度效率较低，对象管理效率这个不好说，但一些函数内部的设计比较优秀。故有些Forcal代码比matlab快，而有些慢。
' V) u9 k  Z5 X6 p
以下例子体现了Forcal和matlab的效率差别所在。

; }8 R7 c$ _. u3 m5 G% {4 }2 J: O% o' v这个matlab程序段是网友lin2009 给出的，理论结果是每个元素均为275000。

1. clear all
   ! p8 K' `\" o; x- {  A5 _3 f( C
2. clc+ |: ^6 K1 o. T( h+ t9 k
   
3. tic
   / k, s# @0 G* ^
4. k = zeros(5,5); % //生成5×5全0矩阵* N1 |& H( B. G! R5 q# x
   
5. % 循环计算以下程序段1000 00次：
     }0 q1 v3 y7 Z2 s6 S
6. for m = 1:1000 00
   2 d3 z2 q7 k) z
7.     a = rand(5,7);( k\" u+ L) U$ r\" H' a
   
8.     b = rand(7,5);%//生成5×7矩阵a,7×5矩阵b,用0～1之间的随机数初始化
     A0 U2 t, z8 M% F0 a9 |
9.     k = k + a * b + a(1:5, 2:6) * b(2:6, 1:5) - a(:, 7) * b(3, :);' k  d- U6 R$ O+ C
   
10. end- M' U# b  R, @- z+ u
    
11. k  x, y# C9 s0 ~& z# c# @% \1 I) ?
    
12. toc

复制代码

Forcal代码1：运行稍快的代码，比matlab约快10％吧？

2. !using["math","sys"];
   
3. , g# ?\\" T) D7 t! |0 s$ p) m9 o
4. mvar:6 o$ ^5 \8 T* ^) z; B7 X
   
6. t0=clock(),4 K( T- O7 q, [. W( B5 |1 n' {\\" H
   
8. oo{k=zeros[5,5]},
   
9. 6 g5 J8 j( L# _. {
10. i=0,((i++)<100000).while{
    
11. 0 K3 _+ p$ z* |! n& u* a! O3 |
12.   oo{: A2 v, v2 v) p1 I. N
    
14.     a=rand[5,7], b=rand[7,5],
    
15. ; q5 }; B' R! c. Z8 a% K  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)]
    
17. + T1 w. {+ J, s4 {$ y4 Q' N0 _
18.   }0 V; J- n3 \& V' F
    
20. },
    
21. $ X5 b3 G- t) x( M  A
22. k.outm(),( I- k& m; [7 h+ U; r& ^
    
24. [clock()-t0]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">, g# ?\\" T) D7 t! |0 s$ p) m9 o</span> mvar:<font class=\\"jammer\\">6 o$ ^5 \8 T* ^) z; B7 X</font><br /> t0=clock(),<font class=\\"jammer\\">4 K( T- O7 q, [. W( B5 |1 n' {\\" H</font><br /> oo{k=zeros[5,5]},<br /> <span style=\\"display:none\\">6 g5 J8 j( L# _. {</span> i=0,((i++)<100000).while{<br /> <span style=\\"display:none\\">0 K3 _+ p$ z* |! n& u* a! O3 |</span>   oo{<font class=\\"jammer\\">: A2 v, v2 v) p1 I. N</font><br />     a=rand[5,7], b=rand[7,5],<br /> <span style=\\"display:none\\">; q5 }; B' R! c. Z8 a% K  v</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\\">+ T1 w. {+ J, s4 {$ y4 Q' N0 _</span>   }<font class=\\"jammer\\">0 V; J- n3 \& V' F</font><br /> },<br /> <span style=\\"display:none\\">$ X5 b3 G- t) x( M  A</span> k.outm(),<font class=\\"jammer\\">( I- k& m; [7 h+ U; r& ^</font><br /> [clock()-t0]/1000;

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

Forcal代码2：比较好看些的代码，似乎也比matlab稍快吧？

2. !using["math","sys"];. ?0 I$ Z# {. M3 Y9 b( e! W! X
   
4. (:t0,k,i,a,b)=* f. d7 R8 b\\" |- w' U
   
6. {9 }. r- k! X2 |
   
8.   t0=clock(),
   
9. + |5 y3 x: T. `: L9 |
10.   oo{k=zeros[5,5]},
    
11. 7 E\\" G\\" t: v' }. K$ T
12.   i=0,((i++)<100000).while{+ D  {1 K8 ^% S/ t6 Y& k  l
    
14.     oo{- T/ w2 {; Q! A8 \( l
    
16.       a=rand[5,7], b=rand[7,5],+ ]7 X4 m, M: k\\" }7 [
    
18.       k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)-a(neg:6)*b(3:neg)
    
19. 0 K6 ?. r+ N5 i( o3 N
20.     }
    
21. $ c0 C\\" ?6 p5 [6 ]
22.   },
    
23. 9 U$ G1 d1 e* _* f$ v* }
24.   k.outm(),$ p0 F: u. y\\" h: Y
    
26.   [clock()-t0]/10005 V$ ^7 Z; R! f/ f% M% n' B9 H
    
28. };

!using["math","sys"];<font class=\\"jammer\\">. ?0 I$ Z# {. M3 Y9 b( e! W! X</font><br /> (:t0,k,i,a,b)=<font class=\\"jammer\\">* f. d7 R8 b\\" |- w' U</font><br /> {<font class=\\"jammer\\">9 }. r- k! X2 |</font><br />   t0=clock(),<br /> <span style=\\"display:none\\">+ |5 y3 x: T. `: L9 |</span>   oo{k=zeros[5,5]},<br /> <span style=\\"display:none\\">7 E\\" G\\" t: v' }. K$ T</span>   i=0,((i++)<100000).while{<font class=\\"jammer\\">+ D  {1 K8 ^% S/ t6 Y& k  l</font><br />     oo{<font class=\\"jammer\\">- T/ w2 {; Q! A8 \( l</font><br />       a=rand[5,7], b=rand[7,5],<font class=\\"jammer\\">+ ]7 X4 m, M: k\\" }7 [</font><br />       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\\">0 K6 ?. r+ N5 i( o3 N</span>     }<br /> <span style=\\"display:none\\">$ c0 C\\" ?6 p5 [6 ]</span>   },<br /> <span style=\\"display:none\\">9 U$ G1 d1 e* _* f$ v* }</span>   k.outm(),<font class=\\"jammer\\">$ p0 F: u. y\\" h: Y</font><br />   [clock()-t0]/1000<font class=\\"jammer\\">5 V$ ^7 Z; R! f/ f% M% n' B9 H</font><br /> };

在我的电脑上运行时间为3.579秒。
+ q0 r% g6 |- t+ l. G% e$ _
例子2：
- e% W- k% v. d; c一段程序的Forcal实现：
7 D/ ?- b4 A  `" q) J

1. //用C++代码描述为：: C5 p  l1 p9 E
   
2. s=0.0;  
     M; Z0 o; I# U9 i. w) z
3. for(x=0.0;x<=1.0;x=x+0.0011)  
   ) Y; }7 |1 L\" m$ U$ Y# \( g
4. {8 j+ g, y. {* T9 a
   
5.   for(y=1.0;y<=2.0;y=y+0.0009)
   7 m( M- g, j\" w3 q, ~0 K! K5 H: \
6.   {
   & m9 R2 W7 `0 a( }4 E9 z. u
7.   s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));  j) G% N, |( b1 s+ C
   
8.   }; v( [( P% P* k. G
   
9. }  

复制代码

结果：
" R4 x" x$ }! @8 k; h1008606.64947441
0.609 //时间
& E( W2 Y7 r" W
, S! a8 M0 [8 A这个matlab程序段是网友yycs001给出的。

1. %file speedtest.m# o. V2 n9 z/ z6 x, K6 B
   
2. function speedtest+ f: E' Z  U3 Z* Z4 c
   
3. format long
   , i: R2 N% o% Z1 T2 T  G7 x
4. tic9 T( Q* Z& Z\" z( m0 W
   
5. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);3 N* ]* }' N; L
   
6. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))0 Y  j5 \5 g2 {/ w2 g
   
7. toc

复制代码

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

2. !using["math","sys"];! v. ^4 C. H9 N
   
4. mvar:
   
5. \\" b  B. t4 K& H' D7 I
6. t=clock(),! E6 i4 v. H! E1 x5 E, V& S! j& g3 F
   
8. oo{( a# j. L# }$ g, w
   
10.   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
11. / v8 _1 Z* A: R% F2 q% t
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]2 R% q, I' k  i6 o
    
14. };0 H* s- C5 \- \( r
    
16. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">! v. ^4 C. H9 N</font><br /> mvar:<br /> <span style=\\"display:none\\">\\" b  B. t4 K& H' D7 I</span> t=clock(),<font class=\\"jammer\\">! E6 i4 v. H! E1 x5 E, V& S! j& g3 F</font><br /> oo{<font class=\\"jammer\\">( a# j. L# }$ g, w</font><br />   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">/ v8 _1 Z* A: R% F2 q% t</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\\">2 R% q, I' k  i6 o</font><br /> };<font class=\\"jammer\\">0 H* s- C5 \- \( r</font><br /> [clock()-t]/1000;

结果：
: Y; B7 h) E0 g1008606.64947441
! F: N7 _. k: T7 A8 S$ n0.625 //时间
0 _. X7 ~. x5 _9 z) K* v' L
或者这个，与上面效率差别不大：
' ]4 }5 U" u" S3 L/ q$ I( [

2. !using["math","sys"];
   
3. ( \' n# X: j7 e5 _0 c& E/ R
4. mvar:( D7 |# N0 O; O- l, ]$ d1 Y  U
   
6. t=clock(),
   
7. ( f- w! l( K. X( @1 `$ N\\" ?3 O1 W
8. oo{
   
9. + I2 }+ y0 s& [) Q+ Q
10.   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
11. \\" m% I. g& g\\" ?8 K. [) n8 J
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. , x. P7 L\\" |: M, f, C: ~3 w# C+ m
14. };9 w, {: w* n6 u5 V) J1 \& K: F7 r& s5 B% L
    
16. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">( \' n# X: j7 e5 _0 c& E/ R</span> mvar:<font class=\\"jammer\\">( D7 |# N0 O; O- l, ]$ d1 Y  U</font><br /> t=clock(),<br /> <span style=\\"display:none\\">( f- w! l( K. X( @1 `$ N\\" ?3 O1 W</span> oo{<br /> <span style=\\"display:none\\">+ I2 }+ y0 s& [) Q+ Q</span>   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">\\" m% I. g& g\\" ?8 K. [) n8 J</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\\">, x. P7 L\\" |: M, f, C: ~3 w# C+ m</span> };<font class=\\"jammer\\">9 w, {: w* n6 u5 V) J1 \& K: F7 r& s5 B% L</font><br /> [clock()-t]/1000;

' p( G0 w9 ~- B6 _* s0 x' OForcal代码2：求和函数sum，非矢量化代码
  N& Y* ^  _+ S0 h) ^

1. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); . k# j- ?: w. \! T7 i
   
2. sum["f",0,1,0.0011,1,2,0.0009];

复制代码

结果：
% `+ K9 Z/ H. y1 n, D( z9 L1008606.64947441
" ]3 v, G8 ^* Y0.719 //时间

$ e7 y9 y: u3 ]/ `5 @Forcal代码3：while循环
6 ?% P& k; }4 l

1. mvar:
   ; J( V- P\" w\" k1 O5 n
2. t=sys::clock();9 L8 t- X7 P5 z, [3 P
   
3. s=0,x=0,
   - P$ N$ d7 T\" L( W. N4 J% C
4. while{x<=1,  //while循环算法； - j& ~8 @* ~0 s: T\" V
   
5.    y=1, 7 r0 N$ T# o) \+ I
   
6.    while{y<=2, 9 Z2 E1 _6 l5 L) \
   
7.        s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))), ' R4 ^6 }' H& h9 `3 E8 r
   
8.        y=y+0.0009
   ' \7 L, N# O$ k1 H6 ~
9.       }, ( a+ O\" Z1 b/ Q\" G- ?7 ~1 y! Y! I' O
   
10.    x=x+0.0011
    8 l9 F% C5 Y% B! w' t
11. }, ( p0 g4 g6 q, U! o' u9 M! q4 A
    
12. s;8 G3 j9 M8 y# U9 H  U
    
13. [sys::clock()-t]/1000;

复制代码

结果：
  ^" b; G. w' l8 Z. u8 A) u% F- L1008606.64947441
' g/ h" H, d$ R) w0.734 //时间

8 [5 P% ~( c5 V* U2 X, F) \大家可下载OpenFC进行测试：http://www.forcal.net/xiazai/forcal9/openfc32w.rar
$ [8 X. s' M7 r! h8 V5 v, {8 b( @3 O
! Y; t# T. \# {1 [7 j. R注意Forcal的矢量化代码第一次运行有时效率较低。
; Y2 ?. j, e7 `, e1 H. r' g
例子1中Forcal和matlab都是矢量化代码，但matlab跑不过Forcal。该例子的特点是函数调用频繁，临时变量生成多，但矩阵很小，矩阵的各种函数运行时耗时较少。故说明Forcal函数调用＋变量管理效率优于matlab。

, z9 O) g/ w) t3 Q/ d例子2中Forcal的矢量化代码是最快的，但与matlab的矢量化代码相比仍有差距。该例子的特点是函数调用少，临时变量也少，但矩阵大。故说明Forcal的各种矩阵函数Sin、Cos及矩阵的加减运算等函数的内部设计不及matlab。
& f  ?4 Q' q& M, L# N8 t
- o: f. V" I% z! q1 w8 {/ ]- P如能在函数内部设计上下点功夫，例子2超越matlab也是可能的。在这方面，期待高手们的指点。

% _3 E3 M" _# h- c如果例子2速度也超越了matlab ，matlab矢量化的神秘面纱就揭开了。
- q( E2 A  F+ [2 ^; a
顺便说一下，例子1如果用C++的运算符重载来实现，速度将比Forcal慢一些，也就是说，在涉及运算符重载时，脚本的效率有时比C++还要高些。

master_math

,,,顶顶更健康！

wznzy0822

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

zhwqqiangge

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

qbist

我 正在  学  Matlab  软件！~~