数学建模社区-数学中国

标题: 关于matlab代码矢量化的理解 [打印本页]

---

作者: forcal    时间: 2010-10-5 09:32
标题: 关于matlab代码矢量化的理解
代码矢量化是matlab的精髓，其基本特点是运行速度快和代码简洁，它是如何实现的？

按我的理解，代码矢量化的本质就是设计专门的函数对数组元素集中运算，这样可提高运行速度，同时兼有代码简洁的特点。

& c! A! G" }5 c# D; W( E! w对matlab的理解比较肤浅，但也确实看不出有更深意义的东西，望解惑。
( O/ }+ C- \) Q
  P- g/ Z, s" B8 [& `! `8 T大家有什么看法，愿畅所欲言。
2 D0 `$ E  r, V( N% D4 t; q0 Q1 |

---

作者: qbist    时间: 2010-10-5 09:52
我 正在  学  Matlab  软件！~~

---

作者: zhwqqiangge    时间: 2010-10-5 10:17
??????????????????

---

作者: wznzy0822    时间: 2010-10-5 10:52
顶。。。。。。。。。。。。。。。。。。。。。。

---

作者: master_math    时间: 2010-10-5 18:32
,,,顶顶更健康！

---

作者: forcal    时间: 2010-10-12 21:36


我正在练手设计的FcMath库也打算以矩阵运算为基础，设计一些专门的函数对数组元素集中运算，运行效率确实有所提高（甚至有些涉及矩阵的算法比matlab还快），代码也简洁了，但不知这是不是矢量化？

脚本运行效率应该取决于函数调度效率、对象管理效率和函数内部算法的实现。
6 J% }& y- @7 z. V$ E/ V) T
我感觉，matlab的函数调度效率较低，对象管理效率这个不好说，但一些函数内部的设计比较优秀。故有些Forcal代码比matlab快，而有些慢。
0 h. _/ K/ J- x' O$ _/ r
& c1 m5 I) D# K  w: g5 k' d以下例子体现了Forcal和matlab的效率差别所在。
$ M* j; k9 P* u" F
这个matlab程序段是网友lin2009 给出的，理论结果是每个元素均为275000。

1. clear all
   
2. clc
   1 W+ X, R# T4 u" u1 N: I  j
3. tic
   $ c2 {9 D. q# B0 J
4. k = zeros(5,5); % //生成5×5全0矩阵
   * i, A& ~' E% o6 _+ u1 H
5. % 循环计算以下程序段1000 00次：
   ; P. E# R. x& _$ I4 \) [! W( [* A
6. for m = 1:1000 00
   
7.     a = rand(5,7);
   * h* r+ w$ `  g3 r2 u
8.     b = rand(7,5);%//生成5×7矩阵a,7×5矩阵b,用0～1之间的随机数初始化
   
9.     k = k + a * b + a(1:5, 2:6) * b(2:6, 1:5) - a(:, 7) * b(3, :);
   
10. end
    
11. k
    
12. toc

复制代码

( C+ Q* b  b( i: x0 wForcal代码1：运行稍快的代码，比matlab约快10％吧？
$ ?/ K; s9 J9 N4 S6 D5 H! ]. E3 a

1. !using["math","sys"];
   
2. mvar:
   
3. t0=clock(),
   
4. oo{k=zeros[5,5]},
   
5. i=0,((i++)<100000).while{
   - N% f: j  |& L( G  j) h( v
6.   oo{
   
7.     a=rand[5,7], b=rand[7,5],
   # I5 C8 s+ I! c& _" G$ a) O7 n8 a
8.     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)]
   
9.   }
   . F$ l, e1 E( R0 g, |6 _6 V
10. },
    
11. k.outm(),
    / Z3 K. t# u- w" {
12. [clock()-t0]/1000;

复制代码

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

8 |' |( s4 Q* h. tForcal代码2：比较好看些的代码，似乎也比matlab稍快吧？

1. !using["math","sys"];
   
2. (:t0,k,i,a,b)=
   1 S/ K) s! a+ F( Y
3. {
   
4.   t0=clock(),
   0 V, ~1 E, H3 Z" z
5.   oo{k=zeros[5,5]},
   5 p( B9 E4 R  J
6.   i=0,((i++)<100000).while{
   
7.     oo{
   
8.       a=rand[5,7], b=rand[7,5],
   2 D* K& Z: P8 u2 }- l# f( I+ A* f
9.       k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)-a(neg:6)*b(3:neg)
   
10.     }
    
11.   },
    
12.   k.outm(),
    ) a( M2 X- j9 n; ^" e3 V2 d
13.   [clock()-t0]/1000
    7 P6 V( l7 X; }& x
14. };

复制代码

在我的电脑上运行时间为3.579秒。
5 z  f9 S$ b) W4 U
8 }  e4 o% ~# P8 J; w& Y* m例子2：
. a6 `4 [' x  t3 |一段程序的Forcal实现：
( R% C) T' t" s  h7 m

1. //用C++代码描述为：
   
2. s=0.0;  
   2 [" U0 W- W8 m1 P& q% W
3. for(x=0.0;x<=1.0;x=x+0.0011)  
   ) ], f( T* @  r
4. {
   % D$ w$ h9 }0 ~: K7 ~( q' `
5.   for(y=1.0;y<=2.0;y=y+0.0009)
   
6.   {
   " q% R$ f: \  D: Y, u+ |
7.   s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   
8.   }
   
9. }  

复制代码

结果：
1008606.64947441
, M' H6 `1 E- k+ Y  l% f: d0.609 //时间

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

1. %file speedtest.m
   
2. function speedtest
   , [7 r; n- p# l2 t- O+ G" C
3. format long
   
4. tic
   6 C1 y% }3 O1 E' Y
5. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   ! L/ f! `" M7 L8 G6 J& I
6. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))
   7 K0 L" L5 `) k" b0 t0 D' ?; }
7. toc

复制代码

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

1. !using["math","sys"];
   
2. mvar:
   ; R6 ?+ e( E5 A8 C
3. t=clock(),
   8 o' O6 {% y- e/ L& ?
4. oo{
   
5.   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
   
6.   Sum[Cos(rn(1)-Sin(rn(1.2)*x^(y/rn(2))+Cos(rn(1)-Sin(rn(1.2)*y^(x/rn(2)))))),0]
   
7. };
   % e+ O: q6 c$ Q+ e) R
8. [clock()-t]/1000;

复制代码

结果：
1008606.64947441
( d9 k* o7 v( E4 k% v4 G+ ^; h0.625 //时间
9 J6 D. q5 o. M0 j5 M$ Y
: Z" U8 K" n8 ^5 D" k, W# c或者这个，与上面效率差别不大：
% t9 I3 l% r9 P, \" v/ G- F- K

1. !using["math","sys"];
   
2. mvar:
   
3. t=clock(),
   
4. oo{
   ; Y) C2 ]% B5 \, c6 N
5.   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
   # T  y" d1 T/ S0 P6 j
6.   Sum[Sum[Cos(rn(1)-Sin(rn(1.2)*x^(y/rn(2))+Cos(rn(1)-Sin(rn(1.2)*y^(x/rn(2))))))]]
   
7. };
   
8. [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)))));
   ' U- }4 ?" ^/ |) ~% G( g* h/ ^
2. sum["f",0,1,0.0011,1,2,0.0009];

复制代码

结果：
1008606.64947441
2 L+ C. Q# d' `2 t0.719 //时间

# e" y# f0 [. W0 OForcal代码3：while循环
1 ?$ q) v4 {4 N: c% i( R# \6 w- l

1. mvar:
   $ T% h% Y" E) @2 W& X( [
2. t=sys::clock();
   , y! M5 X( R" r' M
3. s=0,x=0,
   4 j* u+ s+ k; z0 C& ], `& n
4. while{x<=1,  //while循环算法；
   
5.    y=1,
   ! ?; z/ t3 Q9 _2 b5 n. J% R5 \+ B
6.    while{y<=2,
   
7.        s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
   7 B4 w9 |/ P0 G; q+ J) o6 u, K+ ]! M
8.        y=y+0.0009
   ) [0 v/ N/ ?: f, m8 E
9.       },
   , _: T% {, Q; U8 }
10.    x=x+0.0011
    & Q& D" w8 M8 O# G' r
11. },
    
12. s;
    
13. [sys::clock()-t]/1000;

复制代码

结果：
4 r, p4 t7 R! D$ F( F2 O$ g: I( L1008606.64947441
0.734 //时间
5 j, k9 r+ w4 Y9 _# f2 @+ G
4 \* Q4 p1 I5 w& ~: e大家可下载OpenFC进行测试：http://www.forcal.net/xiazai/forcal9/openfc32w.rar
3 r& v2 I% j: v+ V0 g$ n
注意Forcal的矢量化代码第一次运行有时效率较低。

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

例子2中Forcal的矢量化代码是最快的，但与matlab的矢量化代码相比仍有差距。该例子的特点是函数调用少，临时变量也少，但矩阵大。故说明Forcal的各种矩阵函数Sin、Cos及矩阵的加减运算等函数的内部设计不及matlab。
8 ^  Y: t5 `+ m; m) `
如能在函数内部设计上下点功夫，例子2超越matlab也是可能的。在这方面，期待高手们的指点。

如果例子2速度也超越了matlab ，matlab矢量化的神秘面纱就揭开了。
3 D6 n* a! C4 h2 g6 @
顺便说一下，例子1如果用C++的运算符重载来实现，速度将比Forcal慢一些，也就是说，在涉及运算符重载时，脚本的效率有时比C++还要高些。
& @6 w8 P9 W! c) \$ F6 e

---

作者: forcal    时间: 2010-10-16 16:59
讨论有益！以下是在其他论坛的讨论帖子：
& {+ V2 v* l4 b$ R+ x2 n# Wsimwe：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    时间: 2010-10-24 16:59
参考：http://bbs.emath.ac.cn/thread-2709-1-1.html

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

---

作者: forcal    时间: 2010-11-2 18:06
参考：http://bbs.emath.ac.cn/thread-2727-1-1.html
* ~1 q: j( Y( E# D% a0 [$ B! a
  I2 D2 x. W. ?( C0 F' |5 |关于最快速的矩阵乘实现的讨论。

---

欢迎光临 数学建模社区-数学中国 (http://www.madio.net/)

Powered by Discuz! X2.5