数学建模社区-数学中国

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

---

作者: forcal    时间: 2010-10-5 09:32
标题: 关于matlab代码矢量化的理解
代码矢量化是matlab的精髓，其基本特点是运行速度快和代码简洁，它是如何实现的？
4 c  B, q6 I+ \  v' |- H
按我的理解，代码矢量化的本质就是设计专门的函数对数组元素集中运算，这样可提高运行速度，同时兼有代码简洁的特点。
8 P4 j" C) g* R9 S# L$ d/ E& |) n* ?- A
对matlab的理解比较肤浅，但也确实看不出有更深意义的东西，望解惑。

" t2 n4 ]: h( X% R. R$ K6 L  C5 X! s大家有什么看法，愿畅所欲言。

---

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

, X3 O+ C  ?. P+ t( l: a. c
/ y; y8 u4 A. m我正在练手设计的FcMath库也打算以矩阵运算为基础，设计一些专门的函数对数组元素集中运算，运行效率确实有所提高（甚至有些涉及矩阵的算法比matlab还快），代码也简洁了，但不知这是不是矢量化？

脚本运行效率应该取决于函数调度效率、对象管理效率和函数内部算法的实现。

/ `$ F& G5 z$ M; E& m8 M我感觉，matlab的函数调度效率较低，对象管理效率这个不好说，但一些函数内部的设计比较优秀。故有些Forcal代码比matlab快，而有些慢。
  M# t4 k8 w% W  v
) A, v7 X) c/ a) b1 g4 ]7 c以下例子体现了Forcal和matlab的效率差别所在。
3 A! a- S, G7 {  n. Z
2 d3 H" I( ]! W" w9 Q8 M8 X这个matlab程序段是网友lin2009 给出的，理论结果是每个元素均为275000。

1. clear all
   
2. clc
   6 S2 g! ?8 f9 h
3. tic
   ! B8 S5 N4 S- m7 p4 V
4. k = zeros(5,5); % //生成5×5全0矩阵
   
5. % 循环计算以下程序段1000 00次：
   ! N7 o6 h! T! I  b7 Q* {
6. for m = 1:1000 00
   
7.     a = rand(5,7);
   
8.     b = rand(7,5);%//生成5×7矩阵a,7×5矩阵b,用0～1之间的随机数初始化
   / G4 U2 ]) m6 I% v/ U6 x9 ?# o
9.     k = k + a * b + a(1:5, 2:6) * b(2:6, 1:5) - a(:, 7) * b(3, :);
   
10. end
    
11. k
    + y5 `7 i& T6 X6 {9 U/ `4 r5 P
12. toc

复制代码

Forcal代码1：运行稍快的代码，比matlab约快10％吧？

1. !using["math","sys"];
   
2. mvar:
   
3. t0=clock(),
   
4. oo{k=zeros[5,5]},
   ) U# [5 o: H; D
5. i=0,((i++)<100000).while{
   
6.   oo{
   
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)]
   & K" K4 c0 E8 I  I
9.   }
   ; S" Z# j4 c4 {! X( L& k
10. },
    # z* ]  t8 \& b
11. k.outm(),
    
12. [clock()-t0]/1000;

复制代码

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

: O5 _1 s  w# u: l/ j5 l; ~Forcal代码2：比较好看些的代码，似乎也比matlab稍快吧？
/ y, u6 \$ L$ y: x: t& A

1. !using["math","sys"];
   0 X9 `: x3 u: k2 H: `
2. (:t0,k,i,a,b)=
   9 J7 [  t% C* a* |
3. {
   3 I+ F7 g3 W% |' _9 r2 v- R7 O
4.   t0=clock(),
   ( O( c: `9 A$ `
5.   oo{k=zeros[5,5]},
   5 T: {5 B# ?  S
6.   i=0,((i++)<100000).while{
   
7.     oo{
   
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)
   5 J$ f4 z3 A: Y3 |/ B; c2 w
10.     }
    
11.   },
    2 f" y& ]) s4 f3 j% Q9 ~/ @$ ~
12.   k.outm(),
    
13.   [clock()-t0]/1000
    
14. };

复制代码

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

例子2：
一段程序的Forcal实现：

1. //用C++代码描述为：
   + z5 e/ V# J5 \2 w& Z- z
2. s=0.0;  
   
3. for(x=0.0;x<=1.0;x=x+0.0011)  
   & r# Q, `! E% G3 ]: o9 C
4. {
   2 X/ f: e2 ?$ G2 h( X
5.   for(y=1.0;y<=2.0;y=y+0.0009)
   % g2 K4 U) C! `# n" M( B5 q. [
6.   {
   ; x! C; }# n  W+ m/ z" n
7.   s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   # o6 T% T- G' S4 z
8.   }
   
9. }  

复制代码

结果：
8 @$ ?8 D* K8 K* g1 e5 t1008606.64947441
4 U2 s* S5 L4 g0.609 //时间
+ y' ^8 D' G' g& M1 t' w
! ?, u( C& K) I7 |这个matlab程序段是网友yycs001给出的。
$ ~  Q; V' U8 E  m, k; p0 b

1. %file speedtest.m
   
2. function speedtest
   
3. format long
   
4. tic
   
5. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   " ?4 v; B: U1 Z( [- A. b
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，完全矢量化的代码
6 N/ p- s0 P; I' y7 |

1. !using["math","sys"];
   4 P. \* X) t* S* K* L# N( I* o  E7 N% @
2. mvar:
   2 B6 d6 u7 K1 K$ H
3. t=clock(),
   * e# U7 f, Z/ ^- c
4. oo{
   * E" v; w5 R1 {3 ~8 F' I# Z3 U
5.   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
   & ?3 d+ R5 H5 O/ K5 Q. a
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. };
   
8. [clock()-t]/1000;

复制代码

结果：
! A9 w: n' I! e+ H1 ?8 ]1008606.64947441
" q* J) d- Q$ x& q' @! _" w# U8 b0.625 //时间
* t" V- S. t; n9 F4 z7 y0 [- I
或者这个，与上面效率差别不大：
: @6 f7 O2 o$ l

1. !using["math","sys"];
   
2. mvar:
   1 ?! {2 P# Q9 y4 @
3. t=clock(),
   
4. oo{
   ; k: P$ g. I# r
5.   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
   . {" T* F5 X! r: v4 C# R
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. };
   ! m* y# p: x! i+ s, n8 @9 {
8. [clock()-t]/1000;

复制代码

Forcal代码2：求和函数sum，非矢量化代码
% D2 w; K- }8 U2 \- u9 D# i0 D

1. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   . r7 X5 Q4 |( X/ Q; m  q7 g9 _  c
2. sum["f",0,1,0.0011,1,2,0.0009];

复制代码

结果：
' O+ M* L  R5 L4 B2 y6 V9 W1008606.64947441
; ^) I  D3 h% v" {7 ~0.719 //时间

Forcal代码3：while循环
) s1 i8 q3 J$ ~( e

1. mvar:
   4 q- b( B" L2 _  `1 P9 r
2. t=sys::clock();
   
3. s=0,x=0,
   
4. while{x<=1,  //while循环算法；
     G: I5 n5 h' C1 s
5.    y=1,
   
6.    while{y<=2,
     L( _+ R- P  r1 u6 [, a8 O; A$ P
7.        s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
   
8.        y=y+0.0009
   
9.       },
   & H$ q4 n6 q7 `: S& B4 L
10.    x=x+0.0011
    ; w& y) S/ D( X4 c
11. },
    
12. s;
    ' c+ p9 Q5 }! L, j3 h2 s, n
13. [sys::clock()-t]/1000;

复制代码

结果：
" @; O* x' Y: h6 z* g; r4 ^1008606.64947441
0.734 //时间

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

注意Forcal的矢量化代码第一次运行有时效率较低。
% }$ n# @/ d) i
( ^" y0 u1 ~& ^例子1中Forcal和matlab都是矢量化代码，但matlab跑不过Forcal。该例子的特点是函数调用频繁，临时变量生成多，但矩阵很小，矩阵的各种函数运行时耗时较少。故说明Forcal函数调用＋变量管理效率优于matlab。
1 p7 }7 C; W7 ~2 _" z
  k6 Y. x+ d# S8 e  Y. ]例子2中Forcal的矢量化代码是最快的，但与matlab的矢量化代码相比仍有差距。该例子的特点是函数调用少，临时变量也少，但矩阵大。故说明Forcal的各种矩阵函数Sin、Cos及矩阵的加减运算等函数的内部设计不及matlab。
; H! Y5 V5 r1 @8 C# g* f9 A: H
如能在函数内部设计上下点功夫，例子2超越matlab也是可能的。在这方面，期待高手们的指点。

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

0 ?$ y8 N0 ~2 v7 R顺便说一下，例子1如果用C++的运算符重载来实现，速度将比Forcal慢一些，也就是说，在涉及运算符重载时，脚本的效率有时比C++还要高些。
' N" n2 y4 G8 I

---

作者: forcal    时间: 2010-10-16 16:59
讨论有益！以下是在其他论坛的讨论帖子：
* z+ ]# {5 b+ i( g; K& ~* }simwe：http://forum.simwe.com/thread-952532-1-3.html
6 g. F# x, c: @4 j# K0 @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
7 f& m% R" D, ~  S! U

---

作者: forcal    时间: 2010-10-24 16:59
参考：http://bbs.emath.ac.cn/thread-2709-1-1.html
  H$ Z& b# i( `9 [& W
% m; t# z, }+ c我在多个帖子中有不同说明，但将这些说明再集中到一个帖子中比较麻烦，大家相互参考一下，看能否把这个问题解决了。

---

作者: forcal    时间: 2010-11-2 18:06
参考：http://bbs.emath.ac.cn/thread-2727-1-1.html
, u8 h. L) N( a3 M: I$ Q
关于最快速的矩阵乘实现的讨论。

---

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

Powered by Discuz! X2.5