数学建模社区-数学中国
标题:
关于matlab代码矢量化的理解
[打印本页]
作者:
forcal
时间:
2010-10-5 09:32
标题:
关于matlab代码矢量化的理解
代码矢量化是matlab的精髓,其基本特点是运行速度快和代码简洁,它是如何实现的?
* e6 U9 P1 }# ]( }2 f- x
8 k0 W3 x- F3 w4 h- ?. [
按我的理解,代码矢量化的本质就是设计专门的函数对数组元素集中运算,这样可提高运行速度,同时兼有代码简洁的特点。
* a% m7 i+ Q+ X4 c6 S' `1 h+ p2 G
7 C* o! ^: U; W4 l
对matlab的理解比较肤浅,但也确实看不出有更深意义的东西,望解惑。
$ ~; X5 i+ U" s6 }) T7 p& T1 n
* `9 o5 a. i$ S. A# d) E2 ~8 z
大家有什么看法,愿畅所欲言。
+ V6 ~7 ]+ v' P+ ]" c
作者:
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
本帖最后由 forcal 于 2010-10-12 21:46 编辑
6 y U* q' L$ t+ Y
% ?& q# b" N' N. ^
我正在练手设计的FcMath库也打算以矩阵运算为基础,设计一些专门的函数对数组元素集中运算,运行效率确实有所提高(甚至有些涉及矩阵的算法比matlab还快),代码也简洁了,但不知这是不是矢量化?
9 A: ?6 @" {; @; Y9 _! X9 w
6 m9 ~6 W5 c3 _- d4 D" Q& `) |
脚本运行效率应该取决于函数调度效率、对象管理效率和函数内部算法的实现。
. Q M$ R5 ^0 J" f
. d! p( l( }. e' G; R! b
我感觉,matlab的函数调度效率较低,对象管理效率这个不好说,但一些函数内部的设计比较优秀。故有些Forcal代码比matlab快,而有些慢。
0 I: x5 h: o, W5 a
: A% S! U1 ]+ F
以下例子体现了Forcal和matlab的效率差别所在。
3 i2 }- [! ^( b' Y0 D5 Q
: @, ?. }5 z. \2 Z
这个matlab程序段是网友lin2009 给出的,理论结果是每个元素均为275000。
clear all
0 \# l3 _: \* e' @1 _
clc
# r+ r% J, [; k8 G
tic
) x# M9 Z/ y& I
k = zeros(5,5); % //生成5×5全0矩阵
! c ]6 x" O! Y* s( G: `
% 循环计算以下程序段1000 00次:
& [% D) t# e! c4 W
for m = 1:1000 00
! M& P0 J& F8 ~' _9 L" ]
a = rand(5,7);
H% q6 B/ o1 p6 x% l
b = rand(7,5);%//生成5×7矩阵a,7×5矩阵b,用0~1之间的随机数初始化
: o1 j/ j+ f) |
k = k + a * b + a(1:5, 2:6) * b(2:6, 1:5) - a(:, 7) * b(3, :);
3 Y( _8 u9 j8 S- o1 P
end
! H n |& w( v* S! j0 p
k
+ ?6 ]) C/ f% ^: J& S( f
toc
复制代码
! E# K) C1 o2 p
Forcal代码1:运行稍快的代码,比matlab约快10%吧?
, g+ P( v9 t1 e
!using["math","sys"];
) F, ~5 E1 e7 c, Q
mvar:
* u6 u3 i+ x7 Q' m) ?. T
t0=clock(),
. [, b5 A: ^- Q4 p, R. f
oo{k=zeros[5,5]},
# i+ L2 E9 |; F2 @9 X; l
i=0,((i++)<100000).while{
' \3 \8 D' V/ O, I3 `
oo{
" |9 I' N7 v& I
a=rand[5,7], b=rand[7,5],
* P, K) I' J7 t" v* @
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)]
2 ]! B0 |) [/ s; d+ r2 |7 `9 r
}
& V- z0 o: k, d/ M3 O' b
},
! P$ E/ N1 l; r/ Q' Y
k.outm(),
" h1 F) [. L1 C) @+ v0 ^% Z
[clock()-t0]/1000;
复制代码
在我的电脑上运行时间为3.344秒。
; }$ D3 c/ z0 U1 |
5 L& D- \% l9 `( V' B
Forcal代码2:比较好看些的代码,似乎也比matlab稍快吧?
# J9 @1 q: m7 x7 W/ O5 U
!using["math","sys"];
- y# w7 u" D; S
(:t0,k,i,a,b)=
; h0 ~6 T( f7 ]$ z, N6 `
{
0 N8 A0 C, ~! L
t0=clock(),
& A* `4 A4 d7 o, R& C! Y$ Z
oo{k=zeros[5,5]},
/ g6 P0 {, S1 W
i=0,((i++)<100000).while{
' |$ X; a' n# o4 c+ a
oo{
f3 A: q, |/ a! q
a=rand[5,7], b=rand[7,5],
_4 V; u3 y0 N# e
k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)-a(neg:6)*b(3:neg)
& f' P" K% [8 w, ^) l4 c) x
}
4 v3 [7 @2 h$ a: }, ?5 e2 x. Z( Z
},
0 A2 w0 R8 W0 v
k.outm(),
6 k+ }" `+ u3 ~1 B+ @1 i; s
[clock()-t0]/1000
* Z# r" c, D/ B9 d2 [0 R
};
复制代码
在我的电脑上运行时间为3.579秒。
. D% h( v, v9 I) ]1 N% ?2 t& [
. k" B9 a: U' j( Z
例子2:
7 F6 _& ^1 ^% T+ J( V4 C/ L" j
一段程序的Forcal实现:
: `) ], w$ Z/ [5 r+ x. B; a' L
//用C++代码描述为:
. U5 A+ J: L, t7 z
s=0.0;
7 ~3 W0 A5 o! k" k0 |4 a1 _& X+ @/ G
for(x=0.0;x<=1.0;x=x+0.0011)
, g4 A( [: P) z4 R
{
4 u! n- |7 Y! f& H
for(y=1.0;y<=2.0;y=y+0.0009)
" O/ B1 N; z: o0 f1 j; \! v* A
{
6 n! {8 F" H% |% M$ ^2 V; l9 J
s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
& A0 p8 q; O" ~; ^9 A
}
" n5 ~( W: U# v5 Z. s; O
}
复制代码
结果:
8 }. S4 h4 ?4 G9 v1 Z1 W7 q! }* B0 [5 C
1008606.64947441
5 j& u+ i: y! G5 ^. N3 c; x
0.609 //时间
' j ~' @& T3 M! }5 y8 f) l) `2 }
, D7 V- X. M N4 _* Q; Z
这个matlab程序段是网友yycs001给出的。
$ V0 r: j/ s3 O7 E! \3 L1 h x
%file speedtest.m
% O. I/ {9 E2 U( E$ x
function speedtest
O6 W) l6 k5 l; O5 V l8 d9 ?
format long
! a0 K# k1 A p( R; ? q" f0 m* I
tic
# {! ~, y- B7 g8 x
[x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
0 X2 z2 A! q" p
s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))
- b$ i( p1 t2 G8 X/ W5 ~7 D) D: T5 C( F
toc
复制代码
0 A/ E w1 Y t
Forcal代码1:**数组求和函数Sum,完全矢量化的代码
% T1 O2 r" D0 k* g' p& S9 j% H
!using["math","sys"];
- p, `) I& v1 s1 w% _$ V
mvar:
' \- ~) n: a4 z- o
t=clock(),
% \- \7 M! V3 A) s) D2 G8 b* T1 E
oo{
7 B* ?5 i. M% E) @
ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
. K7 ~7 c4 I# p" _! ^$ \
Sum[Cos(rn(1)-Sin(rn(1.2)*x^(y/rn(2))+Cos(rn(1)-Sin(rn(1.2)*y^(x/rn(2)))))),0]
9 v, _3 c/ l( Z& m5 l; ]9 i
};
% l9 B& S/ E. I) B% } @; M$ s! Z
[clock()-t]/1000;
复制代码
结果:
4 ?. ^- J7 i+ F3 y
1008606.64947441
& r: Z0 W- u6 n$ E3 i- s
0.625 //时间
) X) l; w. O8 V1 G/ q" J9 R
8 O) v3 H8 e* V. i" z ^
或者这个,与上面效率差别不大:
/ R6 ^9 \6 u! |/ c" c
!using["math","sys"];
o* `4 V% `. x0 X9 M: j) U3 u9 w
mvar:
0 y" K" Y; ], e' n
t=clock(),
# W* Y9 z& t9 g* V
oo{
8 D0 Y) F- D; [8 _' K! e/ U
ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
7 q6 C3 ?* T. j" ~; E
Sum[Sum[Cos(rn(1)-Sin(rn(1.2)*x^(y/rn(2))+Cos(rn(1)-Sin(rn(1.2)*y^(x/rn(2))))))]]
r; t9 V3 U G# Z
};
) a7 z6 e5 [2 A. _. B9 p
[clock()-t]/1000;
复制代码
) N( ~/ \6 ^! B; t
Forcal代码2:求和函数sum,非矢量化代码
' D! P7 `3 A. Z$ ?
f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
" R3 q1 P% p, `
sum["f",0,1,0.0011,1,2,0.0009];
复制代码
结果:
2 I$ k" H, W$ o
1008606.64947441
4 p: d, R- l' Y
0.719 //时间
; S: d- A, ?8 c, v# U2 n; |
: J" F% y0 N0 Q, |% }4 \& O
Forcal代码3:while循环
* r% Z' k. i/ d% u0 B# \" Q$ P
mvar:
7 |/ A H8 r4 }0 @
t=sys::clock();
3 E$ t* ], x$ @2 F- X& P
s=0,x=0,
/ x, Q, `" R ~) y7 c8 o4 U7 y
while{x<=1, //while循环算法;
( t. O' N6 P9 p, b. F4 m3 f
y=1,
+ j: E! `4 s( G) _) n9 w- c+ h+ z
while{y<=2,
" H6 O! [" x- J5 [: ]8 t: x, a3 X
s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
5 _, S% e; }8 C9 ^9 M& u# p: Z) I
y=y+0.0009
' [7 p) `' [) f# f! ~
},
3 H. Q& G7 A, D0 l; {. y
x=x+0.0011
: u0 k; Y7 M3 S3 j7 ^
},
: h3 E8 g: X1 G) ^. L
s;
/ r2 d; S; T( n
[sys::clock()-t]/1000;
复制代码
结果:
+ k: P# c7 F6 `# b; \5 {7 U
1008606.64947441
8 y) u% q) w* O. V0 x( J
0.734 //时间
. y9 h+ l( [, ^! z
; z8 O- b* v* a; C6 x7 I% V
大家可下载OpenFC进行测试:
http://www.forcal.net/xiazai/forcal9/openfc32w.rar
/ U9 p$ s' [0 D
. k U+ R3 G1 x# Q1 `. o+ ?/ X& e
注意Forcal的矢量化代码第一次运行有时效率较低。
& m# I+ \, _7 f% a
* X/ r' Y. J0 S( y- C
例子1中Forcal和matlab都是矢量化代码,但matlab跑不过Forcal。该例子的特点是函数调用频繁,临时变量生成多,但矩阵很小,矩阵的各种函数运行时耗时较少。故说明Forcal函数调用+变量管理效率优于matlab。
2 Z- w8 K( S5 x5 t' T" u8 U- a. i
, `( [. d0 a- r {
例子2中Forcal的矢量化代码是最快的,但与matlab的矢量化代码相比仍有差距。该例子的特点是函数调用少,临时变量也少,但矩阵大。故说明Forcal的各种矩阵函数Sin、Cos及矩阵的加减运算等函数的内部设计不及matlab。
/ y; A/ G5 o+ i, U) w* d
$ F2 |" l/ r) t1 Y
如能在函数内部设计上下点功夫,例子2超越matlab也是可能的。在这方面,期待高手们的指点。
2 }) a) t$ ~5 {& \
J1 r$ p7 H% Q6 s
如果例子2速度也超越了matlab ,matlab矢量化的神秘面纱就揭开了。
: U% l2 r1 d7 R) B l/ h, j
9 {2 a8 `. N$ w
顺便说一下,例子1如果用C++的运算符重载来实现,速度将比Forcal慢一些,也就是说,在涉及运算符重载时,脚本的效率有时比C++还要高些。
% ?+ @) T( i1 |* _1 m
作者:
forcal
时间:
2010-10-16 16:59
讨论有益!以下是在其他论坛的讨论帖子:
' a/ d* A1 g/ t1 ^
simwe:
http://forum.simwe.com/thread-952532-1-3.html
% S0 h% ~" ]! {0 V8 s
csdn:
http://topic.csdn.net/u/20101006/21/2ed9e5c1-cc9f-4623-b1c0-ebd5d1b5a98a.html
& z, [& c" U1 n0 ~; D! k
cadn:
http://topic.csdn.net/u/20101010/15/3bcf2fe0-0abd-4c29-b854-b1d007b16863.html
, x2 \* t' O+ c/ d: D7 N) _- K" b, g
作者:
forcal
时间:
2010-10-24 16:59
参考:
http://bbs.emath.ac.cn/thread-2709-1-1.html
. Y4 t) Z. b- J% G
* s$ a [# f' [. Z
我在多个帖子中有不同说明,但将这些说明再集中到一个帖子中比较麻烦,大家相互参考一下,看能否把这个问题解决了。
作者:
forcal
时间:
2010-11-2 18:06
参考:
http://bbs.emath.ac.cn/thread-2727-1-1.html
& L1 j3 ^# l/ ~& a3 o4 l
5 P s* L- U: U* F
关于最快速的矩阵乘实现的讨论。
% Z+ F( X3 b. ?6 R: W
欢迎光临 数学建模社区-数学中国 (http://www.madio.net/)
Powered by Discuz! X2.5