关于matlab代码矢量化的理解

forcal

代码矢量化是matlab的精髓，其基本特点是运行速度快和代码简洁，它是如何实现的？
# `3 f& h" O, T( j+ `$ D1 G
按我的理解，代码矢量化的本质就是设计专门的函数对数组元素集中运算，这样可提高运行速度，同时兼有代码简洁的特点。

2 o: h' u. c) W  F. Q4 u. I! K( |# x对matlab的理解比较肤浅，但也确实看不出有更深意义的东西，望解惑。

3 X6 S' D) d. T大家有什么看法，愿畅所欲言。

qbist

我 正在  学  Matlab  软件！~~

zhwqqiangge

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

wznzy0822

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

master_math

,,,顶顶更健康！

forcal

$ Y- Z9 e; W: m9 r+ [
1 ]. _# E/ Z0 `: H  a" G0 Q我正在练手设计的FcMath库也打算以矩阵运算为基础，设计一些专门的函数对数组元素集中运算，运行效率确实有所提高（甚至有些涉及矩阵的算法比matlab还快），代码也简洁了，但不知这是不是矢量化？
) g/ M5 G' D  q4 Z# v! R3 t
脚本运行效率应该取决于函数调度效率、对象管理效率和函数内部算法的实现。
, H/ c2 b7 `( I% s! e) C" W
我感觉，matlab的函数调度效率较低，对象管理效率这个不好说，但一些函数内部的设计比较优秀。故有些Forcal代码比matlab快，而有些慢。
' j6 y4 `. ~8 A1 u
. r) ]3 m) p- }- `. \. Y以下例子体现了Forcal和matlab的效率差别所在。

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

1. clear all
   / L% N( }3 U+ B5 v+ A5 j- `% [
2. clc. {0 s5 m+ Y* X2 [$ k( X  e! o1 C
   
3. tic
   2 \6 e\" ^( E7 N& I
4. k = zeros(5,5); % //生成5×5全0矩阵! O* K/ L: e/ U3 u
   
5. % 循环计算以下程序段1000 00次：
   * E\" i- l+ i, n+ @
6. for m = 1:1000 00& X1 Y1 I2 }6 M
   
7.     a = rand(5,7);
   $ a) f3 ~' p, ?3 \. x
8.     b = rand(7,5);%//生成5×7矩阵a,7×5矩阵b,用0～1之间的随机数初始化
   1 V5 W7 r* T0 j
9.     k = k + a * b + a(1:5, 2:6) * b(2:6, 1:5) - a(:, 7) * b(3, :);
   ! m# {8 L+ j% m& o& m
10. end$ B( @3 b- `+ x3 V+ t4 ^; F
    
11. k9 j% ?) f* v\" m' g& N
    
12. toc

复制代码

: t" s. j( P- p+ A/ {( |% j5 VForcal代码1：运行稍快的代码，比matlab约快10％吧？
( }; x0 x- X- N

2. !using["math","sys"];
   
3. 9 q% |' A' I: i
4. mvar:, b$ Z( Z6 r/ f
   
6. t0=clock(),
   
7. * W\\" [\\" Z8 l\\" l
8. oo{k=zeros[5,5]},
   
9. # g/ w/ \5 g) S+ G0 g3 d. \, o
10. i=0,((i++)<100000).while{% ^& M, f! l8 {4 F7 t
    
12.   oo{2 _0 e. _% _# _
    
14.     a=rand[5,7], b=rand[7,5],
    
15. 7 G3 M$ W, X7 M7 O3 E
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)]$ [  N& k2 h0 L$ ^
    
18.   }. o% L' S, H\\" ]! K9 q
    
20. },
    
21. ( ]$ H9 a) S8 t) `% P
22. k.outm(),
    
23. % C& p) Q! f0 ^% V2 u+ ~+ B: i! O
24. [clock()-t0]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">9 q% |' A' I: i</span> mvar:<font class=\\"jammer\\">, b$ Z( Z6 r/ f</font><br /> t0=clock(),<br /> <span style=\\"display:none\\">* W\\" [\\" Z8 l\\" l</span> oo{k=zeros[5,5]},<br /> <span style=\\"display:none\\"># g/ w/ \5 g) S+ G0 g3 d. \, o</span> i=0,((i++)<100000).while{<font class=\\"jammer\\">% ^& M, f! l8 {4 F7 t</font><br />   oo{<font class=\\"jammer\\">2 _0 e. _% _# _</font><br />     a=rand[5,7], b=rand[7,5],<br /> <span style=\\"display:none\\">7 G3 M$ W, X7 M7 O3 E</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\\">$ [  N& k2 h0 L$ ^</font><br />   }<font class=\\"jammer\\">. o% L' S, H\\" ]! K9 q</font><br /> },<br /> <span style=\\"display:none\\">( ]$ H9 a) S8 t) `% P</span> k.outm(),<br /> <span style=\\"display:none\\">% C& p) Q! f0 ^% V2 u+ ~+ B: i! O</span> [clock()-t0]/1000;

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

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

2. !using["math","sys"];8 S' t  a2 m, ~% w/ k4 c' p
   
4. (:t0,k,i,a,b)=* O1 {% l; Y6 p8 o
   
6. {) O( V' A* c* t! T- o
   
8.   t0=clock(),$ a\\" \1 \* @+ d5 P5 x8 x/ K+ Q
   
10.   oo{k=zeros[5,5]},
    
11. 5 g$ ^; a3 C* R2 H
12.   i=0,((i++)<100000).while{' y8 {. t: M, U7 S# \
    
14.     oo{9 b' V6 X' ]0 g0 y8 Z4 U8 e2 A
    
16.       a=rand[5,7], b=rand[7,5],
    
17. 5 w, Q/ J) o; }. h
18.       k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)-a(neg:6)*b(3:neg)# e2 X7 R  U1 T- [
    
20.     }
    
21. $ Z: L* R& w  V
22.   },
    
23. ) W$ R7 e7 `8 C* `, M
24.   k.outm(),1 }, `\\" L; v) b# |6 E1 l
    
26.   [clock()-t0]/1000# b+ y* L) |. W( Y! M, x& B, q& b. G
    
28. };

!using["math","sys"];<font class=\\"jammer\\">8 S' t  a2 m, ~% w/ k4 c' p</font><br /> (:t0,k,i,a,b)=<font class=\\"jammer\\">* O1 {% l; Y6 p8 o</font><br /> {<font class=\\"jammer\\">) O( V' A* c* t! T- o</font><br />   t0=clock(),<font class=\\"jammer\\">$ a\\" \1 \* @+ d5 P5 x8 x/ K+ Q</font><br />   oo{k=zeros[5,5]},<br /> <span style=\\"display:none\\">5 g$ ^; a3 C* R2 H</span>   i=0,((i++)<100000).while{<font class=\\"jammer\\">' y8 {. t: M, U7 S# \</font><br />     oo{<font class=\\"jammer\\">9 b' V6 X' ]0 g0 y8 Z4 U8 e2 A</font><br />       a=rand[5,7], b=rand[7,5],<br /> <span style=\\"display:none\\">5 w, Q/ J) o; }. h</span>       k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)-a(neg:6)*b(3:neg)<font class=\\"jammer\\"># e2 X7 R  U1 T- [</font><br />     }<br /> <span style=\\"display:none\\">$ Z: L* R& w  V</span>   },<br /> <span style=\\"display:none\\">) W$ R7 e7 `8 C* `, M</span>   k.outm(),<font class=\\"jammer\\">1 }, `\\" L; v) b# |6 E1 l</font><br />   [clock()-t0]/1000<font class=\\"jammer\\"># b+ y* L) |. W( Y! M, x& B, q& b. G</font><br /> };

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

例子2：
* p: v& {: I/ g0 ~一段程序的Forcal实现：
9 l8 g1 A" }0 v! i% W

1. //用C++代码描述为：
   . ~, h2 Q7 s+ t- R
2. s=0.0;  
   % q- h: u- ?* ]\" ~$ t5 Y
3. for(x=0.0;x<=1.0;x=x+0.0011)  1 j# v+ `5 x- I  P; c
   
4. {: ^1 s7 u' }% W; X& N- H
   
5.   for(y=1.0;y<=2.0;y=y+0.0009)
   % `2 M) Q& U4 O& [; J  \
6.   {  e  ], t% r4 U3 |
   
7.   s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));0 \+ n5 L# z  F' b* i
   
8.   }6 y& j2 Q6 W% s
   
9. }  

复制代码

结果：
( v% e1 u8 C$ r1008606.64947441
0.609 //时间

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

1. %file speedtest.m; ~; g& r: M8 N\" T/ v
   
2. function speedtest6 F: f# S' d* d& O- f' @. k
   
3. format long# F- S3 T6 V; H- v
   
4. tic
   4 N4 X. J; }5 C4 n# U% N
5. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);: I2 `9 }6 L- W2 b$ |! n
   
6. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))
   2 {5 K2 b5 J4 v6 g6 l/ J
7. toc

复制代码

4 U- d, ^7 {0 ~1 eForcal代码1：**数组求和函数Sum，完全矢量化的代码

2. !using["math","sys"];
   
3. : T5 G- ~) H9 Y6 g
4. mvar:3 z3 ~. Z& ?3 ]3 n: b  H1 W4 T5 Y
   
6. t=clock(),
   
7. 5 A$ D/ }3 y' d3 b/ d4 N& u
8. oo{3 S' L2 Y4 W$ V8 ~8 p1 s* b$ o
   
10.   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
11. 8 s; ?6 u: W# f' N
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. ) u3 ~6 `/ X; o- i
14. };0 X1 O( t/ P, I! `6 a
    
16. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">: T5 G- ~) H9 Y6 g</span> mvar:<font class=\\"jammer\\">3 z3 ~. Z& ?3 ]3 n: b  H1 W4 T5 Y</font><br /> t=clock(),<br /> <span style=\\"display:none\\">5 A$ D/ }3 y' d3 b/ d4 N& u</span> oo{<font class=\\"jammer\\">3 S' L2 Y4 W$ V8 ~8 p1 s* b$ o</font><br />   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">8 s; ?6 u: W# f' N</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]<br /> <span style=\\"display:none\\">) u3 ~6 `/ X; o- i</span> };<font class=\\"jammer\\">0 X1 O( t/ P, I! `6 a</font><br /> [clock()-t]/1000;

结果：
) M4 q) W" B8 o- M1008606.64947441
0.625 //时间
6 _" x+ h" V% |/ a# r8 P" S
% ?* g) V* @& I4 n3 e4 Y; H或者这个，与上面效率差别不大：

2. !using["math","sys"];/ u/ }3 a+ S& h3 o5 x* C
   
4. mvar:! l2 p9 s. n% ^* E6 c6 N6 q
   
6. t=clock(),/ n, h: g, ~( G& Y' m\\" [) B; R
   
8. oo{! a2 F8 q! b- n1 K/ D4 J
   
10.   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],% M4 d5 n! n0 p# t& o4 T* ?
    
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))))))]]( Y! V9 h% b) `* @. z/ i: P- r$ L
    
14. };, y4 b  R: O4 |. Y
    
16. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">/ u/ }3 a+ S& h3 o5 x* C</font><br /> mvar:<font class=\\"jammer\\">! l2 p9 s. n% ^* E6 c6 N6 q</font><br /> t=clock(),<font class=\\"jammer\\">/ n, h: g, ~( G& Y' m\\" [) B; R</font><br /> oo{<font class=\\"jammer\\">! a2 F8 q! b- n1 K/ D4 J</font><br />   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">% M4 d5 n! n0 p# t& o4 T* ?</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))))))]]<font class=\\"jammer\\">( Y! V9 h% b) `* @. z/ i: P- r$ L</font><br /> };<font class=\\"jammer\\">, y4 b  R: O4 |. Y</font><br /> [clock()-t]/1000;

Forcal代码2：求和函数sum，非矢量化代码
/ N; [$ c$ t6 ?. p5 o" @; o

1. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); 7 c9 F% U0 Q- N) _9 p) o  ?
   
2. sum["f",0,1,0.0011,1,2,0.0009];

复制代码

结果：
. k' _6 y1 c) g! z0 k1008606.64947441
  j* _" J# `. `' n, Z0.719 //时间
! m' j% u; S/ B
Forcal代码3：while循环

1. mvar:5 f. u4 }+ Y8 x4 G3 H
   
2. t=sys::clock();
   ( ^) g- B3 \( e
3. s=0,x=0,
   2 b0 b  h0 I1 s
4. while{x<=1,  //while循环算法； 1 e) T6 Q& ~+ z# N6 p9 Y& O2 q
   
5.    y=1, ' K$ k; O3 V0 l( B1 {
   
6.    while{y<=2,
   . Z5 b/ O0 @2 F1 T5 l! c
7.        s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
   ( p; k( f7 Z& i! A- C\" [% S
8.        y=y+0.0009 \" H2 m/ W( d8 G- h' N
   
9.       },
   9 g3 L; {' N) ?( T; i' l
10.    x=x+0.0011
    / p\" A- n8 s' W7 k( r) y) `
11. },
    % H) h& M1 d' q% p1 b. h: Z
12. s;
    1 \: _. w. w! i
13. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
- Q! W* ?5 J4 F, n, R  _! V4 T0.734 //时间
, I; {% g( i9 \1 y( y! J
大家可下载OpenFC进行测试：http://www.forcal.net/xiazai/forcal9/openfc32w.rar
1 y2 ]3 x4 N- \' }
注意Forcal的矢量化代码第一次运行有时效率较低。

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

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

如能在函数内部设计上下点功夫，例子2超越matlab也是可能的。在这方面，期待高手们的指点。

+ Y' O& C, a* P5 m( y如果例子2速度也超越了matlab ，matlab矢量化的神秘面纱就揭开了。
2 i1 f: v3 ^8 B4 o- R  H# V# j
- }6 y' n+ @2 Y3 i4 T: l+ j2 O顺便说一下，例子1如果用C++的运算符重载来实现，速度将比Forcal慢一些，也就是说，在涉及运算符重载时，脚本的效率有时比C++还要高些。

forcal

讨论有益！以下是在其他论坛的讨论帖子：
  y+ ^' _8 }4 [( i1 g9 g; @simwe：http://forum.simwe.com/thread-952532-1-3.html
% P1 K' G$ `- Tcsdn：http://topic.csdn.net/u/20101006/21/2ed9e5c1-cc9f-4623-b1c0-ebd5d1b5a98a.html
% ^7 f! g1 ?- |9 ^8 G) [% Ocadn：http://topic.csdn.net/u/20101010/15/3bcf2fe0-0abd-4c29-b854-b1d007b16863.html
7 {8 F7 o  V. B" K

forcal

参考：http://bbs.emath.ac.cn/thread-2709-1-1.html
3 q- s# M8 U8 }
8 G: [( T  H! c6 M) l5 O# y. C4 A+ l我在多个帖子中有不同说明，但将这些说明再集中到一个帖子中比较麻烦，大家相互参考一下，看能否把这个问题解决了。

forcal

参考：http://bbs.emath.ac.cn/thread-2727-1-1.html
% u! v* }3 Q$ t) w" \
3 G% P* |( L& D9 C# ]( Z( R关于最快速的矩阵乘实现的讨论。