关于matlab代码矢量化的理解

forcal

代码矢量化是matlab的精髓，其基本特点是运行速度快和代码简洁，它是如何实现的？
' X9 {4 P$ U) v  z4 R
$ p# \1 E. i7 h5 H& s# n5 l  R按我的理解，代码矢量化的本质就是设计专门的函数对数组元素集中运算，这样可提高运行速度，同时兼有代码简洁的特点。
  _& ^! t$ O: d6 N
对matlab的理解比较肤浅，但也确实看不出有更深意义的东西，望解惑。
, X( l8 x! D) z9 z
4 ?% j! V0 t+ Y0 _9 g; ]% [* @大家有什么看法，愿畅所欲言。

qbist

我 正在  学  Matlab  软件！~~

zhwqqiangge

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

wznzy0822

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

master_math

,,,顶顶更健康！

forcal

# v% ?. q, s( a( k/ L
( u; A! s/ U# I9 k; n7 x8 t我正在练手设计的FcMath库也打算以矩阵运算为基础，设计一些专门的函数对数组元素集中运算，运行效率确实有所提高（甚至有些涉及矩阵的算法比matlab还快），代码也简洁了，但不知这是不是矢量化？

脚本运行效率应该取决于函数调度效率、对象管理效率和函数内部算法的实现。
' h. F" E0 ~2 t& n
我感觉，matlab的函数调度效率较低，对象管理效率这个不好说，但一些函数内部的设计比较优秀。故有些Forcal代码比matlab快，而有些慢。
3 T( ?  N  u. [8 N  {  k
以下例子体现了Forcal和matlab的效率差别所在。

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

1. clear all
   ! k! L: @' O) l: r
2. clc
   ! ^; m2 T0 ^9 E, Z\" ~7 h1 |
3. tic
   # n7 J* e  k8 b
4. k = zeros(5,5); % //生成5×5全0矩阵) ]& f, F; m( L: w/ E, a
   
5. % 循环计算以下程序段1000 00次：
   5 C& v  v) t) A( t1 G. I7 }+ A
6. for m = 1:1000 00! A$ @$ K9 F# G  v6 }9 t& Z: J
   
7.     a = rand(5,7);! J% [7 Z$ i. ~1 A  [6 i2 T, c' M
   
8.     b = rand(7,5);%//生成5×7矩阵a,7×5矩阵b,用0～1之间的随机数初始化' R4 {: v5 k. X+ u5 \
   
9.     k = k + a * b + a(1:5, 2:6) * b(2:6, 1:5) - a(:, 7) * b(3, :);( [9 x/ a3 K: F; y; X, ~
   
10. end- D% c2 u* l) U+ p. y( s\" r5 G3 {
    
11. k/ k  U  d5 i9 O+ w( g/ F6 D
    
12. toc

复制代码

" V9 m0 K3 ~9 SForcal代码1：运行稍快的代码，比matlab约快10％吧？
- O# p; V& e0 o6 {  a2 N

2. !using["math","sys"];% }& Q) _  ]1 W- x
   
4. mvar:
   
5. 1 }' `8 ^& R# L( a  ^
6. t0=clock(),) f+ [6 {* R) |; Y* Q
   
8. oo{k=zeros[5,5]},
   
9. . W! c0 K' g6 N3 u
10. i=0,((i++)<100000).while{
    
11. : ~5 {$ s3 Q2 Y
12.   oo{
    
13. \\" V3 H2 @5 b7 s( N) t: k- l
14.     a=rand[5,7], b=rand[7,5],
    
15. 7 M8 t- s( y: z
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. / h5 t- @\\" l- C
18.   }
    
19. ' X& ~0 p, Y6 v) W
20. },2 P6 M\\" O\\" p, ]2 v( S
    
22. k.outm(),
    
23. 6 `$ ^  j) ^, {8 N5 k
24. [clock()-t0]/1000;

!using["math","sys"];<font class=\\"jammer\\">% }& Q) _  ]1 W- x</font><br /> mvar:<br /> <span style=\\"display:none\\">1 }' `8 ^& R# L( a  ^</span> t0=clock(),<font class=\\"jammer\\">) f+ [6 {* R) |; Y* Q</font><br /> oo{k=zeros[5,5]},<br /> <span style=\\"display:none\\">. W! c0 K' g6 N3 u</span> i=0,((i++)<100000).while{<br /> <span style=\\"display:none\\">: ~5 {$ s3 Q2 Y</span>   oo{<br /> <span style=\\"display:none\\">\\" V3 H2 @5 b7 s( N) t: k- l</span>     a=rand[5,7], b=rand[7,5],<br /> <span style=\\"display:none\\">7 M8 t- s( y: z</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\\">/ h5 t- @\\" l- C</span>   }<br /> <span style=\\"display:none\\">' X& ~0 p, Y6 v) W</span> },<font class=\\"jammer\\">2 P6 M\\" O\\" p, ]2 v( S</font><br /> k.outm(),<br /> <span style=\\"display:none\\">6 `$ ^  j) ^, {8 N5 k</span> [clock()-t0]/1000;

在我的电脑上运行时间为3.344秒。
1 s( F1 e" I. j- V6 C- Z) s
8 N& j; [! d0 T: e" r4 U1 N. MForcal代码2：比较好看些的代码，似乎也比matlab稍快吧？
1 g: u& \( G' A, A

2. !using["math","sys"];3 {2 ^; r9 T+ ]! b! I
   
4. (:t0,k,i,a,b)=
   
5. ) q\\" c7 \4 a2 w) Q
6. {
   
7. ) }4 |; }8 f$ V. W4 ~
8.   t0=clock(),\\" \  Y1 [3 b. v3 S
   
10.   oo{k=zeros[5,5]},, b9 z4 s) X9 U4 q2 @
    
12.   i=0,((i++)<100000).while{
    
13. \\" A3 j: Q. h3 j4 H( A- s
14.     oo{
    
15. ; t' {( U1 q) G% O1 ^
16.       a=rand[5,7], b=rand[7,5],/ z! |/ }9 T# E
    
18.       k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)-a(neg:6)*b(3:neg)/ m  {; b' c# _9 l) J- W4 G\\" y
    
20.     }- y5 X% t' X% @) z1 J& V1 M/ T
    
22.   },3 L' I( G8 n# R: P; g% Z6 X' s
    
24.   k.outm(),
    
25. 2 q( d% l4 M$ ]2 t! j
26.   [clock()-t0]/10000 ]\\" K' c% P( x6 i9 n
    
28. };

!using["math","sys"];<font class=\\"jammer\\">3 {2 ^; r9 T+ ]! b! I</font><br /> (:t0,k,i,a,b)=<br /> <span style=\\"display:none\\">) q\\" c7 \4 a2 w) Q</span> {<br /> <span style=\\"display:none\\">) }4 |; }8 f$ V. W4 ~</span>   t0=clock(),<font class=\\"jammer\\">\\" \  Y1 [3 b. v3 S</font><br />   oo{k=zeros[5,5]},<font class=\\"jammer\\">, b9 z4 s) X9 U4 q2 @</font><br />   i=0,((i++)<100000).while{<br /> <span style=\\"display:none\\">\\" A3 j: Q. h3 j4 H( A- s</span>     oo{<br /> <span style=\\"display:none\\">; t' {( U1 q) G% O1 ^</span>       a=rand[5,7], b=rand[7,5],<font class=\\"jammer\\">/ z! |/ }9 T# E</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\\">/ m  {; b' c# _9 l) J- W4 G\\" y</font><br />     }<font class=\\"jammer\\">- y5 X% t' X% @) z1 J& V1 M/ T</font><br />   },<font class=\\"jammer\\">3 L' I( G8 n# R: P; g% Z6 X' s</font><br />   k.outm(),<br /> <span style=\\"display:none\\">2 q( d% l4 M$ ]2 t! j</span>   [clock()-t0]/1000<font class=\\"jammer\\">0 ]\\" K' c% P( x6 i9 n</font><br /> };

在我的电脑上运行时间为3.579秒。
  x, K7 ^7 k, `) o, V5 O; j
: T* K" P9 V# p* H例子2：
: e( K+ E% r' e) `" p- c. l: s# F7 H一段程序的Forcal实现：

1. //用C++代码描述为：* B1 L' X) k$ G/ V5 s; E
   
2. s=0.0;  
   ! r+ l. P; W. `! ^; D
3. for(x=0.0;x<=1.0;x=x+0.0011)  4 Q+ i9 d( v1 P' V
   
4. {
   $ \* m( C$ a2 m$ M# ~\" ?3 e& k
5.   for(y=1.0;y<=2.0;y=y+0.0009). v6 h8 n4 {% a. Y) L5 r
   
6.   {! a2 y\" M6 ?/ m* L& W9 c
   
7.   s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   9 h' @+ J0 \: H, F1 m: i
8.   }
   , o1 \+ F\" S7 G( w; |
9. }  

复制代码

结果：
4 J$ u8 N  N% p$ U1008606.64947441
0.609 //时间
; g4 L/ j; Y. F" c0 c1 j! F0 I
0 T  D9 K) {0 C) b: k! ^这个matlab程序段是网友yycs001给出的。
2 H9 I+ x: A: v: }7 V/ ~3 J

1. %file speedtest.m
   , [) _. R6 e- z  c9 s
2. function speedtest* i9 ~  B/ x  u: i+ D1 J
   
3. format long' Q' [* t+ R+ b; s7 }9 A
   
4. tic# |+ p% E5 S7 d\" Q$ a
   
5. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);2 q. _. y& R# C, p+ C, Z( D) a$ L
   
6. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))
   2 D4 c, S/ q% C0 w9 P8 {2 l+ O; S
7. toc

复制代码

Forcal代码1：**数组求和函数Sum，完全矢量化的代码
6 V; |( K* `* l

2. !using["math","sys"];) y5 ?/ {5 m1 ?1 _: s0 P
   
4. mvar:' k# `8 n- x- V  A0 l  K\\" T
   
6. t=clock(),
   
7. 5 D2 P& h% O$ H, c
8. oo{& Q  _: J& J: w7 i4 ]
   
10.   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],5 m\\" O8 R) u1 x7 [3 J
    
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]: X% N/ ~& G, I. j) x
    
14. };0 j5 F! Y7 d8 @6 O2 j1 n
    
16. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">) y5 ?/ {5 m1 ?1 _: s0 P</font><br /> mvar:<font class=\\"jammer\\">' k# `8 n- x- V  A0 l  K\\" T</font><br /> t=clock(),<br /> <span style=\\"display:none\\">5 D2 P& h% O$ H, c</span> oo{<font class=\\"jammer\\">& Q  _: J& J: w7 i4 ]</font><br />   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">5 m\\" O8 R) u1 x7 [3 J</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]<font class=\\"jammer\\">: X% N/ ~& G, I. j) x</font><br /> };<font class=\\"jammer\\">0 j5 F! Y7 d8 @6 O2 j1 n</font><br /> [clock()-t]/1000;

结果：
1008606.64947441
; `; G! @3 k: }2 ^0.625 //时间

6 F" G7 M; e& L或者这个，与上面效率差别不大：
* P( |) l+ B- t2 V1 m2 b% p$ m

2. !using["math","sys"];
   
3. % w* w, O$ J3 ^\\" s- C9 B
4. mvar:
   
5. % N' E2 l: d+ b  r* N9 B
6. t=clock(),
   
7. 9 F' y. U  J0 m/ D! F
8. oo{
   
9. / `+ o& y6 `! o- H! x/ S- b1 l9 ]
10.   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
11. : a1 q1 t1 j+ i  ?/ O* y; ?* N+ 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.   c: U0 c) P4 l
14. };# U6 A0 v; y- L$ |! ?5 r$ f+ ?
    
16. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">% w* w, O$ J3 ^\\" s- C9 B</span> mvar:<br /> <span style=\\"display:none\\">% N' E2 l: d+ b  r* N9 B</span> t=clock(),<br /> <span style=\\"display:none\\">9 F' y. U  J0 m/ D! F</span> oo{<br /> <span style=\\"display:none\\">/ `+ o& y6 `! o- H! x/ S- b1 l9 ]</span>   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">: a1 q1 t1 j+ i  ?/ O* y; ?* N+ W</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\\">  c: U0 c) P4 l</span> };<font class=\\"jammer\\"># U6 A0 v; y- L$ |! ?5 r$ f+ ?</font><br /> [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))))); $ K$ H9 c5 N# V
   
2. sum["f",0,1,0.0011,1,2,0.0009];

复制代码

结果：
1008606.64947441
0.719 //时间

Forcal代码3：while循环
1 d0 g& m" n6 k/ |

1. mvar:1 h) @& o\" Q5 r; H
   
2. t=sys::clock();
   $ `* r* b. B- k\" K5 s' E3 J7 r) I
3. s=0,x=0, $ ]5 S# P8 H2 \4 `1 A' \7 D
   
4. while{x<=1,  //while循环算法； . ]1 X* P9 Q* q9 x) H2 |
   
5.    y=1,
   8 j6 B  o6 n; c/ J
6.    while{y<=2, ' S3 G* e( [\" s) B- R
   
7.        s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))), / _9 ?* l6 X2 ^+ `
   
8.        y=y+0.0009 # }7 C. f\" |5 c# @& q
   
9.       }, \" g% C0 Z- o$ i; b  ~6 R( t& q\" z
   
10.    x=x+0.0011   q+ f. s1 }7 N/ f
    
11. },
    * |+ }1 P! K9 l0 x& Q6 y# N$ E
12. s;+ w5 T6 H  L% p
    
13. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
, o0 o3 f. P3 L3 K0.734 //时间
- ^  j, s8 o; P# p. j+ y
大家可下载OpenFC进行测试：http://www.forcal.net/xiazai/forcal9/openfc32w.rar
. P/ b8 v8 q% W9 y
* ?+ {2 `0 G" t/ w4 V0 m0 U* H注意Forcal的矢量化代码第一次运行有时效率较低。
  l; A, R+ ^0 s/ E# M1 N; J
例子1中Forcal和matlab都是矢量化代码，但matlab跑不过Forcal。该例子的特点是函数调用频繁，临时变量生成多，但矩阵很小，矩阵的各种函数运行时耗时较少。故说明Forcal函数调用＋变量管理效率优于matlab。

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

' o" h5 k' i3 G  f; {! X; R如能在函数内部设计上下点功夫，例子2超越matlab也是可能的。在这方面，期待高手们的指点。
* q3 _# N, W6 h) j0 B7 X5 d' x
  h& L) y4 w: P* C2 Z如果例子2速度也超越了matlab ，matlab矢量化的神秘面纱就揭开了。

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

forcal

讨论有益！以下是在其他论坛的讨论帖子：
5 g9 T9 }6 V. ]+ L8 U- K9 R; N& I5 `simwe：http://forum.simwe.com/thread-952532-1-3.html
9 Y# x' f" R3 ^* y6 W/ y9 }csdn：http://topic.csdn.net/u/20101006/21/2ed9e5c1-cc9f-4623-b1c0-ebd5d1b5a98a.html
: y% L! ]1 l" C0 d8 i; g# k  Dcadn：http://topic.csdn.net/u/20101010/15/3bcf2fe0-0abd-4c29-b854-b1d007b16863.html

forcal

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

4 S6 e( n* }  H% B8 d" f我在多个帖子中有不同说明，但将这些说明再集中到一个帖子中比较麻烦，大家相互参考一下，看能否把这个问题解决了。

forcal

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

关于最快速的矩阵乘实现的讨论。
8 I  i# H( N9 q: _9 l2 N, }5 t