数学建模社区-数学中国

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

---

作者: forcal    时间: 2010-10-5 09:32
标题: 关于matlab代码矢量化的理解
代码矢量化是matlab的精髓，其基本特点是运行速度快和代码简洁，它是如何实现的？
3 B1 G  t3 J0 [, S- W5 W3 F* ^* m
按我的理解，代码矢量化的本质就是设计专门的函数对数组元素集中运算，这样可提高运行速度，同时兼有代码简洁的特点。

对matlab的理解比较肤浅，但也确实看不出有更深意义的东西，望解惑。
: K* n5 |& @8 W, a$ a. i, @" s3 J
大家有什么看法，愿畅所欲言。

---

作者: 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


0 g3 ?, Q7 p; M5 `5 A1 x我正在练手设计的FcMath库也打算以矩阵运算为基础，设计一些专门的函数对数组元素集中运算，运行效率确实有所提高（甚至有些涉及矩阵的算法比matlab还快），代码也简洁了，但不知这是不是矢量化？
2 [2 @$ C* l! x0 z; a4 {, f
. w& A* E8 f4 N# S2 i! w1 L脚本运行效率应该取决于函数调度效率、对象管理效率和函数内部算法的实现。

我感觉，matlab的函数调度效率较低，对象管理效率这个不好说，但一些函数内部的设计比较优秀。故有些Forcal代码比matlab快，而有些慢。
# q/ p3 x$ V, b6 q# j( p
& k6 \" i$ D5 B( i1 X以下例子体现了Forcal和matlab的效率差别所在。
( N0 |- ^4 P; m- T: @8 i
0 c$ K4 o2 E) V* g- X. o5 \& e这个matlab程序段是网友lin2009 给出的，理论结果是每个元素均为275000。

1. clear all
   
2. clc
   " B. ]9 i- O2 u5 A3 [/ T+ y
3. tic
   . E, ~% L2 A7 i, `$ [
4. k = zeros(5,5); % //生成5×5全0矩阵
   
5. % 循环计算以下程序段1000 00次：
   
6. for m = 1:1000 00
   8 H. a/ w" ~8 y4 k1 Y( i7 C/ W
7.     a = rand(5,7);
     L# }7 X, H* u7 j' H6 U& [( w4 L
8.     b = rand(7,5);%//生成5×7矩阵a,7×5矩阵b,用0～1之间的随机数初始化
   / i4 N0 X& Y9 s# o3 F/ Q- P0 R
9.     k = k + a * b + a(1:5, 2:6) * b(2:6, 1:5) - a(:, 7) * b(3, :);
   
10. end
    
11. k
      H6 F: ~8 q- p  C/ v) _
12. toc

复制代码

. d1 v' ?, ^% c3 {* T3 hForcal代码1：运行稍快的代码，比matlab约快10％吧？
8 O: X3 q6 s' H# i+ V+ W4 Y

1. !using["math","sys"];
   
2. mvar:
   - f, q( _2 U2 Y
3. t0=clock(),
   
4. oo{k=zeros[5,5]},
   3 ?" ^$ L* X6 ~8 Z) V) V; U
5. i=0,((i++)<100000).while{
   2 m; S# g7 m$ V/ o. Y. Q2 Y
6.   oo{
   ' E- C* ~1 g2 }" p
7.     a=rand[5,7], b=rand[7,5],
   ; T  t% ~! X# @: e
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.   }
   $ j' j7 x) {* A3 N2 X$ ^- u& ~% _
10. },
    , {6 {! o  g1 p- Y+ \$ }
11. k.outm(),
    4 W% O) \! h2 j+ l' M1 k: u
12. [clock()-t0]/1000;

复制代码

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

! F9 {" j& O, {3 t* F- U! fForcal代码2：比较好看些的代码，似乎也比matlab稍快吧？

1. !using["math","sys"];
   ! }' a4 ?8 @' }: b
2. (:t0,k,i,a,b)=
   
3. {
   / i0 K0 P! D  M: A% m. K0 X
4.   t0=clock(),
   
5.   oo{k=zeros[5,5]},
   
6.   i=0,((i++)<100000).while{
   * A, m$ A7 x9 h2 n" N, G
7.     oo{
   5 |# K& ^6 A; z# c6 E# f( O' Y! e
8.       a=rand[5,7], b=rand[7,5],
   $ \/ X1 y0 f0 i: j0 S  S
9.       k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)-a(neg:6)*b(3:neg)
   
10.     }
    8 u1 x% a6 U3 `8 O/ s% Q, Y' t1 c
11.   },
    
12.   k.outm(),
    . H% p" w2 M- q! a& {+ i  K! c* v
13.   [clock()-t0]/1000
    6 r4 B. z8 V7 |9 t2 W! o0 v! _. S
14. };

复制代码

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

例子2：
1 `3 w! H$ J9 @: l一段程序的Forcal实现：
! }- z  n$ G6 i

1. //用C++代码描述为：
   
2. s=0.0;  
   
3. for(x=0.0;x<=1.0;x=x+0.0011)  
   
4. {
   
5.   for(y=1.0;y<=2.0;y=y+0.0009)
   
6.   {
   
7.   s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   
8.   }
   
9. }  

复制代码

结果：
# j# T+ {0 G% B5 H! y1008606.64947441
$ U9 N$ _" g1 ^0.609 //时间
3 E$ ?- f3 e1 t2 S$ o* w
这个matlab程序段是网友yycs001给出的。
' l& L+ F8 v" F0 A/ h! a

1. %file speedtest.m
   
2. function speedtest
   
3. format long
   8 a+ b' B/ H( ~3 ^
4. tic
   ; N$ k) f/ W" }/ K7 N
5. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   
6. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))
   
7. toc

复制代码

; O+ K9 v- }+ @$ W) tForcal代码1：**数组求和函数Sum，完全矢量化的代码
, X1 b: J3 Y8 p9 r) o4 H

1. !using["math","sys"];
   
2. mvar:
   
3. t=clock(),
   0 |: B% H) @  g
4. oo{
   8 m$ ~8 u; V8 G, _# w
5.   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
   4 b; j, L! V( O8 C& s: z% B
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]
   / ~* m6 _- u7 D6 }
7. };
   
8. [clock()-t]/1000;

复制代码

结果：
' r2 h3 T8 ^8 I( h( {% ?0 R! h1008606.64947441
0.625 //时间

或者这个，与上面效率差别不大：
% @6 M; I; s3 s

1. !using["math","sys"];
   * T& p$ b% j/ V1 v
2. mvar:
     C4 K) c; i7 \. M7 x* b0 b
3. t=clock(),
   
4. oo{
   + x: J* ]) V, L6 Y2 c
5.   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
   5 ~! R" U' \4 j0 q
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;

复制代码

* V4 `4 a( ]4 ]/ u7 z9 wForcal代码2：求和函数sum，非矢量化代码

1. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   * R2 C/ P4 {1 ?! [
2. sum["f",0,1,0.0011,1,2,0.0009];

复制代码

结果：
% M' p2 Y2 ~. {8 n/ o6 ?% x0 ^# q1008606.64947441
" _; B* I9 p0 a- d8 |8 e5 C: v0.719 //时间
% f  ~, M. A  N% U# G
+ d; I: K# a" W; {) X/ L% Y+ {4 uForcal代码3：while循环
0 g0 M8 \5 A6 ]1 [

1. mvar:
   
2. t=sys::clock();
   
3. s=0,x=0,
   ; d& Z' p& T! u( U
4. while{x<=1,  //while循环算法；
   " e% B- }% ]0 z2 I2 I
5.    y=1,
   
6.    while{y<=2,
   
7.        s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
   
8.        y=y+0.0009
   - s, ?" r6 r' v/ ^' \
9.       },
   & e" {! e0 \$ z8 e
10.    x=x+0.0011
    
11. },
    & K* O& ^2 W3 D& x  K( h$ _+ O; {
12. s;
    
13. [sys::clock()-t]/1000;

复制代码

结果：
2 L$ n: v- I! V3 S1008606.64947441
" w! E! z2 t  a8 J# i0.734 //时间

: c- ?6 ~6 A  }" g7 ?+ Y大家可下载OpenFC进行测试：http://www.forcal.net/xiazai/forcal9/openfc32w.rar

注意Forcal的矢量化代码第一次运行有时效率较低。

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

/ [3 W3 E4 X% h+ a例子2中Forcal的矢量化代码是最快的，但与matlab的矢量化代码相比仍有差距。该例子的特点是函数调用少，临时变量也少，但矩阵大。故说明Forcal的各种矩阵函数Sin、Cos及矩阵的加减运算等函数的内部设计不及matlab。
' `( S' j, P0 V
如能在函数内部设计上下点功夫，例子2超越matlab也是可能的。在这方面，期待高手们的指点。

如果例子2速度也超越了matlab ，matlab矢量化的神秘面纱就揭开了。

) K! b  z- B4 v1 }9 A: S顺便说一下，例子1如果用C++的运算符重载来实现，速度将比Forcal慢一些，也就是说，在涉及运算符重载时，脚本的效率有时比C++还要高些。

---

作者: forcal    时间: 2010-10-16 16:59
讨论有益！以下是在其他论坛的讨论帖子：
simwe：http://forum.simwe.com/thread-952532-1-3.html
csdn：http://topic.csdn.net/u/20101006/21/2ed9e5c1-cc9f-4623-b1c0-ebd5d1b5a98a.html
3 K7 D5 X: |* A0 Rcadn：http://topic.csdn.net/u/20101010/15/3bcf2fe0-0abd-4c29-b854-b1d007b16863.html
4 U1 K$ l. `* H

---

作者: 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
; i2 J6 Y5 a" o2 w
2 d* }1 K5 t$ j7 j  [. l关于最快速的矩阵乘实现的讨论。

---

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

Powered by Discuz! X2.5