数学建模社区-数学中国

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

---

作者: forcal    时间: 2010-10-5 09:32
标题: 关于matlab代码矢量化的理解
代码矢量化是matlab的精髓，其基本特点是运行速度快和代码简洁，它是如何实现的？
3 h" C6 C. s$ N4 [
按我的理解，代码矢量化的本质就是设计专门的函数对数组元素集中运算，这样可提高运行速度，同时兼有代码简洁的特点。

对matlab的理解比较肤浅，但也确实看不出有更深意义的东西，望解惑。

大家有什么看法，愿畅所欲言。
! C7 C, l5 N; M* \/ ^4 H* Z

---

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

6 P, i3 Y8 V4 ~
% X& c( C" b7 M6 G7 u我正在练手设计的FcMath库也打算以矩阵运算为基础，设计一些专门的函数对数组元素集中运算，运行效率确实有所提高（甚至有些涉及矩阵的算法比matlab还快），代码也简洁了，但不知这是不是矢量化？
5 |0 e9 v8 D  o7 p
" P0 s7 F5 z, ^8 v' o1 m# J# m; \3 t脚本运行效率应该取决于函数调度效率、对象管理效率和函数内部算法的实现。
0 |6 J4 }' o4 s, ^6 D2 L! x
我感觉，matlab的函数调度效率较低，对象管理效率这个不好说，但一些函数内部的设计比较优秀。故有些Forcal代码比matlab快，而有些慢。
4 s6 _3 [8 t% t' t
, z! E, Z  F9 x. X5 h: U5 Y  b! p, g以下例子体现了Forcal和matlab的效率差别所在。

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

1. clear all
   
2. clc
   6 F5 o6 C, }1 x4 V
3. tic
   * ^( H; f0 ~& V' ]5 h9 U: v$ o" [
4. k = zeros(5,5); % //生成5×5全0矩阵
   1 i5 ^& }# h( p( T* c) z' W8 ]
5. % 循环计算以下程序段1000 00次：
   
6. for m = 1:1000 00
   . a1 E4 l. J2 j$ L
7.     a = rand(5,7);
   9 P0 a0 ^7 \2 [! h9 p
8.     b = rand(7,5);%//生成5×7矩阵a,7×5矩阵b,用0～1之间的随机数初始化
   % A+ u5 U* X: K/ `" Z2 B
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

复制代码

* v( W; k5 h- r1 S" D% d# qForcal代码1：运行稍快的代码，比matlab约快10％吧？
- p. [+ d. @$ N' V

1. !using["math","sys"];
   : o; L3 z! [9 _: W2 ~9 S
2. mvar:
   5 F( Z& q0 U8 a9 i6 ]' U
3. t0=clock(),
   
4. oo{k=zeros[5,5]},
   ( `. \, e! L* t. o" B1 b
5. i=0,((i++)<100000).while{
   
6.   oo{
   1 L, J6 J# Y8 F2 ]; Y8 ^
7.     a=rand[5,7], b=rand[7,5],
   
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)]
   & o+ }% I" t) _4 i# I) X) b" u
9.   }
   , z5 }% t* n  |( Z+ D2 }
10. },
    
11. k.outm(),
    * E! e& }( n9 N7 b% T, \" L
12. [clock()-t0]/1000;

复制代码

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

: ?$ p1 p6 W6 A* v! U0 @" P& e: ~* ZForcal代码2：比较好看些的代码，似乎也比matlab稍快吧？
/ q) a- o' O% M& D

1. !using["math","sys"];
   0 l, k3 J. _  z+ `$ G2 i- U, h8 Q8 a
2. (:t0,k,i,a,b)=
   
3. {
     K6 X$ S# e3 V, ^8 y+ l0 r' ?
4.   t0=clock(),
   # ~3 E# ^5 R' q3 b  J  ?
5.   oo{k=zeros[5,5]},
   % l9 ^- |% p* m/ l- y5 L+ J
6.   i=0,((i++)<100000).while{
   
7.     oo{
   3 `3 D- ~2 Y  J  r7 h  C+ k9 l) q
8.       a=rand[5,7], b=rand[7,5],
   
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(),
    # w# i& D5 u- R. v
13.   [clock()-t0]/1000
    2 @2 @+ o# v: G
14. };

复制代码

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

例子2：
1 r- o6 [# p3 k一段程序的Forcal实现：

1. //用C++代码描述为：
   * Q- @7 E  _; y# _
2. s=0.0;  
   4 z8 p' G" N8 z
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.   }
   8 T2 G, V0 H- N# k* \# Z
9. }  

复制代码

结果：
, k4 F% V) d- D8 d$ B1008606.64947441
0.609 //时间
; A, [/ ^+ O3 i& G+ A0 H& |
这个matlab程序段是网友yycs001给出的。
; l0 @+ q  o! f, ^, R

1. %file speedtest.m
   
2. function speedtest
   
3. format long
   
4. tic
   
5. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   * y, A4 Q5 o$ h% ~6 H( I
6. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))
   
7. toc

复制代码

Forcal代码1：**数组求和函数Sum，完全矢量化的代码
+ X* u( z3 P: ]; {, X: q$ G

1. !using["math","sys"];
   - Z6 Q# T! w: X( F& B. G
2. mvar:
   9 f. w( ~6 {' e
3. t=clock(),
   . f- j, ~* A+ W% U
4. oo{
   3 J6 }( n! L) |5 J7 Z. h
5.   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
   # I& n  W, h" d2 C1 U7 a0 u$ \
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]
   4 h- v7 o/ }1 ^# R  E
7. };
   ! q5 C5 a2 @% W; `, ~+ D
8. [clock()-t]/1000;

复制代码

结果：
5 @( g# `2 d3 B2 d! h) Y7 ?% U5 O: C! [1008606.64947441
3 V2 E" N1 u# o9 f. C  _0.625 //时间

或者这个，与上面效率差别不大：

1. !using["math","sys"];
   
2. mvar:
   
3. t=clock(),
   
4. oo{
   0 i9 l7 E* F! o4 A5 U
5.   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
   ! c5 k$ c0 Y% ?  o; H7 @
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))))))]]
   ; ]& Q: S& z; R- W5 v
7. };
   
8. [clock()-t]/1000;

复制代码

Forcal代码2：求和函数sum，非矢量化代码
% p' g: N  y2 T" k2 X0 T" A

1. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   ) h1 b+ @% D$ d# G! Y8 i
2. sum["f",0,1,0.0011,1,2,0.0009];

复制代码

结果：
1008606.64947441
1 \9 [8 S: [* e  ~# Y; W0.719 //时间

Forcal代码3：while循环

1. mvar:
   ( q9 q+ D6 P7 D6 G8 m. M; B
2. t=sys::clock();
   
3. s=0,x=0,
   
4. while{x<=1,  //while循环算法；
   % p0 `7 g5 H8 F  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))))),
   . X' t3 [0 A8 S" x( W- U4 q  F
8.        y=y+0.0009
   1 j8 m/ ?6 j- B
9.       },
   
10.    x=x+0.0011
    
11. },
    
12. s;
    
13. [sys::clock()-t]/1000;

复制代码

结果：
0 J2 t% `8 ?. \: Y9 h& K5 V1008606.64947441
3 G% Z1 V& d3 k# U/ |3 U0.734 //时间

大家可下载OpenFC进行测试：http://www.forcal.net/xiazai/forcal9/openfc32w.rar

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

  n; P) P5 b: v+ X! ]例子1中Forcal和matlab都是矢量化代码，但matlab跑不过Forcal。该例子的特点是函数调用频繁，临时变量生成多，但矩阵很小，矩阵的各种函数运行时耗时较少。故说明Forcal函数调用＋变量管理效率优于matlab。
4 M8 ~* d- m2 L! p" A) U! {" Q
例子2中Forcal的矢量化代码是最快的，但与matlab的矢量化代码相比仍有差距。该例子的特点是函数调用少，临时变量也少，但矩阵大。故说明Forcal的各种矩阵函数Sin、Cos及矩阵的加减运算等函数的内部设计不及matlab。

2 V4 Z4 q4 Q2 o如能在函数内部设计上下点功夫，例子2超越matlab也是可能的。在这方面，期待高手们的指点。
/ h, f+ B0 i8 ?4 j0 Y
& y$ s+ P! E3 l: y8 `如果例子2速度也超越了matlab ，matlab矢量化的神秘面纱就揭开了。

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

---

作者: forcal    时间: 2010-10-16 16:59
讨论有益！以下是在其他论坛的讨论帖子：
' |0 M8 T  L7 v; t( Osimwe：http://forum.simwe.com/thread-952532-1-3.html
csdn：http://topic.csdn.net/u/20101006/21/2ed9e5c1-cc9f-4623-b1c0-ebd5d1b5a98a.html
/ m9 f7 p. i+ E9 M+ x+ y$ N2 Qcadn：http://topic.csdn.net/u/20101010/15/3bcf2fe0-0abd-4c29-b854-b1d007b16863.html
5 l1 v: s& k0 U' C$ g" V! q0 @

---

作者: forcal    时间: 2010-10-24 16:59
参考：http://bbs.emath.ac.cn/thread-2709-1-1.html
( m4 T! c- J6 [+ p
+ W0 n3 C5 @7 f$ i) |3 }我在多个帖子中有不同说明，但将这些说明再集中到一个帖子中比较麻烦，大家相互参考一下，看能否把这个问题解决了。

---

作者: forcal    时间: 2010-11-2 18:06
参考：http://bbs.emath.ac.cn/thread-2727-1-1.html
: n* [$ s0 i) |
0 p! |* H% f. o7 ~* \- h- X关于最快速的矩阵乘实现的讨论。
- K+ F2 l" F- A; ?( F

---

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

Powered by Discuz! X2.5