数学建模社区-数学中国
标题:
关于matlab代码矢量化的理解
[打印本页]
作者:
forcal
时间:
2010-10-5 09:32
标题:
关于matlab代码矢量化的理解
代码矢量化是matlab的精髓,其基本特点是运行速度快和代码简洁,它是如何实现的?
4 c B, q6 I+ \ v' |- H
0 w. _; m: I+ q* M# y
按我的理解,代码矢量化的本质就是设计专门的函数对数组元素集中运算,这样可提高运行速度,同时兼有代码简洁的特点。
8 P4 j" C) g* R9 S# L$ d/ E& |) n* ?- A
7 p9 J# O% {5 `7 G4 N
对matlab的理解比较肤浅,但也确实看不出有更深意义的东西,望解惑。
- Q: q. s, Y `3 w m
" t2 n4 ]: h( X% R. R$ K6 L C5 X! s
大家有什么看法,愿畅所欲言。
0 r* \$ y/ M" O0 @
作者:
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 编辑
, X3 O+ C ?. P+ t( l: a. c
/ y; y8 u4 A. m
我正在练手设计的FcMath库也打算以矩阵运算为基础,设计一些专门的函数对数组元素集中运算,运行效率确实有所提高(甚至有些涉及矩阵的算法比matlab还快),代码也简洁了,但不知这是不是矢量化?
1 z! o7 w( u& {8 e
+ N! O; K: [9 N7 z- v4 Z
脚本运行效率应该取决于函数调度效率、对象管理效率和函数内部算法的实现。
* H' ^6 ^: n, ]
/ `$ 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。
clear all
' R, e I' t/ C w+ ?# U- a
clc
6 S2 g! ?8 f9 h
tic
! B8 S5 N4 S- m7 p4 V
k = zeros(5,5); % //生成5×5全0矩阵
% L7 s4 ], d. f8 T, r$ b) q
% 循环计算以下程序段1000 00次:
! N7 o6 h! T! I b7 Q* {
for m = 1:1000 00
7 T# m' U, S: v* h
a = rand(5,7);
3 ]2 I6 e* N7 J' c8 M( G
b = rand(7,5);%//生成5×7矩阵a,7×5矩阵b,用0~1之间的随机数初始化
/ G4 U2 ]) m6 I% v/ U6 x9 ?# o
k = k + a * b + a(1:5, 2:6) * b(2:6, 1:5) - a(:, 7) * b(3, :);
9 P6 R b1 I8 `$ a
end
. k# B* C9 [! T+ f# U, c
k
+ y5 `7 i& T6 X6 {9 U/ `4 r5 P
toc
复制代码
4 N8 Z. g( b B$ {+ y4 ?) {! H. e
Forcal代码1:运行稍快的代码,比matlab约快10%吧?
0 k7 m$ A/ o1 t$ O
!using["math","sys"];
4 G; R$ x2 e; I1 b0 x; c' J' s: u
mvar:
) O% |, E$ V- g
t0=clock(),
5 P2 O% Y2 ?9 {" V/ g6 ]2 c, w
oo{k=zeros[5,5]},
) U# [5 o: H; D
i=0,((i++)<100000).while{
, x% D" k0 X: g- @, R
oo{
( A' _3 p+ {1 c, Z
a=rand[5,7], b=rand[7,5],
- B4 ]6 @' u6 E% |# R4 ^+ P
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
}
; S" Z# j4 c4 {! X( L& k
},
# z* ] t8 \& b
k.outm(),
4 y1 w6 Z& P/ D9 z& t$ N; K
[clock()-t0]/1000;
复制代码
在我的电脑上运行时间为3.344秒。
# l* p9 R! V1 W* _1 e, h- y& O
: O5 _1 s w# u: l/ j5 l; ~
Forcal代码2:比较好看些的代码,似乎也比matlab稍快吧?
/ y, u6 \$ L$ y: x: t& A
!using["math","sys"];
0 X9 `: x3 u: k2 H: `
(:t0,k,i,a,b)=
9 J7 [ t% C* a* |
{
3 I+ F7 g3 W% |' _9 r2 v- R7 O
t0=clock(),
( O( c: `9 A$ `
oo{k=zeros[5,5]},
5 T: {5 B# ? S
i=0,((i++)<100000).while{
" g7 G. p0 G: b- e3 Q& `
oo{
2 `3 c, F! C- I0 K
a=rand[5,7], b=rand[7,5],
1 U/ D6 \0 u9 D+ K
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
}
T# W4 ?5 n4 L0 m
},
2 f" y& ]) s4 f3 j% Q9 ~/ @$ ~
k.outm(),
! L# W: O4 D( D6 F' _% Z0 @7 d( ]
[clock()-t0]/1000
$ a% f$ s( b7 _. z. ?
};
复制代码
在我的电脑上运行时间为3.579秒。
5 ?0 d' @# U2 H. K2 S
* A2 O( c2 ?3 r- Z# p, J+ N; u
例子2:
; H% S7 |! g; Z$ l
一段程序的Forcal实现:
# q4 Z6 h2 ]5 d- z+ _
//用C++代码描述为:
+ z5 e/ V# J5 \2 w& Z- z
s=0.0;
: h$ W5 w; h, ?9 ]8 W' p
for(x=0.0;x<=1.0;x=x+0.0011)
& r# Q, `! E% G3 ]: o9 C
{
2 X/ f: e2 ?$ G2 h( X
for(y=1.0;y<=2.0;y=y+0.0009)
% g2 K4 U) C! `# n" M( B5 q. [
{
; x! C; }# n W+ m/ z" n
s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
# o6 T% T- G' S4 z
}
% o( j2 S, K3 t& @
}
复制代码
结果:
8 @$ ?8 D* K8 K* g1 e5 t
1008606.64947441
4 U2 s* S5 L4 g
0.609 //时间
+ y' ^8 D' G' g& M1 t' w
! ?, u( C& K) I7 |
这个matlab程序段是网友yycs001给出的。
$ ~ Q; V' U8 E m, k; p0 b
%file speedtest.m
* L8 k# C4 O/ `1 X" d$ E
function speedtest
B E+ |7 B: r( @" B9 Y1 o
format long
9 S0 `, e+ h# r0 h: n- e3 i
tic
, Z) o# @. G; [. [, U9 ?
[x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
" ?4 v; B: U1 Z( [- A. b
s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))
- Z+ g+ y! r5 M4 o0 K; Z7 h
toc
复制代码
0 M- D7 q8 T+ @
Forcal代码1:**数组求和函数Sum,完全矢量化的代码
6 N/ p- s0 P; I' y7 |
!using["math","sys"];
4 P. \* X) t* S* K* L# N( I* o E7 N% @
mvar:
2 B6 d6 u7 K1 K$ H
t=clock(),
* e# U7 f, Z/ ^- c
oo{
* E" v; w5 R1 {3 ~8 F' I# Z3 U
ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
& ?3 d+ R5 H5 O/ K5 Q. a
Sum[Cos(rn(1)-Sin(rn(1.2)*x^(y/rn(2))+Cos(rn(1)-Sin(rn(1.2)*y^(x/rn(2)))))),0]
6 F: p: d% Y" G& o/ g; K+ y7 _: V
};
8 D+ a# C* n( F: r% d9 t4 L$ G9 b
[clock()-t]/1000;
复制代码
结果:
! A9 w: n' I! e+ H1 ?8 ]
1008606.64947441
" q* J) d- Q$ x& q' @! _" w# U8 b
0.625 //时间
* t" V- S. t; n9 F4 z7 y0 [- I
9 B6 h6 M& L3 v" k. j* @
或者这个,与上面效率差别不大:
: @6 f7 O2 o$ l
!using["math","sys"];
* k5 d/ M$ y/ k. j& s% L* F/ r5 M
mvar:
1 ?! {2 P# Q9 y4 @
t=clock(),
6 l% p2 W- _5 W
oo{
; k: P$ g. I# r
ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
. {" T* F5 X! r: v4 C# R
Sum[Sum[Cos(rn(1)-Sin(rn(1.2)*x^(y/rn(2))+Cos(rn(1)-Sin(rn(1.2)*y^(x/rn(2))))))]]
! h$ M+ w$ G$ W- P$ f$ F! [' ` @
};
! m* y# p: x! i+ s, n8 @9 {
[clock()-t]/1000;
复制代码
) T4 Z. P: \! {4 w# v6 H3 ^9 |
Forcal代码2:求和函数sum,非矢量化代码
% D2 w; K- }8 U2 \- u9 D# i0 D
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
sum["f",0,1,0.0011,1,2,0.0009];
复制代码
结果:
' O+ M* L R5 L4 B2 y6 V9 W
1008606.64947441
; ^) I D3 h% v" {7 ~
0.719 //时间
7 q+ S2 [. U4 C& e2 z& H' T
: H. x- `" w! E& D
Forcal代码3:while循环
) s1 i8 q3 J$ ~( e
mvar:
4 q- b( B" L2 _ `1 P9 r
t=sys::clock();
4 a" V5 v5 t7 ~9 A6 g
s=0,x=0,
1 F. E4 s' W5 W
while{x<=1, //while循环算法;
G: I5 n5 h' C1 s
y=1,
: h5 L) y! n: i3 ]9 ?( p' [
while{y<=2,
L( _+ R- P r1 u6 [, a8 O; A$ P
s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
7 e# @3 V' K% W3 T. _# k( h2 P6 b
y=y+0.0009
+ ^& A; P: U5 E8 F- D2 [$ B2 h! D
},
& H$ q4 n6 q7 `: S& B4 L
x=x+0.0011
; w& y) S/ D( X4 c
},
( d8 ^7 c$ z& Q" M: V- I/ h; q
s;
' c+ p9 Q5 }! L, j3 h2 s, n
[sys::clock()-t]/1000;
复制代码
结果:
" @; O* x' Y: h6 z* g; r4 ^
1008606.64947441
' c% ~6 `% ~1 H* k s8 t
0.734 //时间
6 L4 C% T& M3 m4 Y+ m
0 s. w/ { }1 A) r. Y# E# | d! X
大家可下载OpenFC进行测试:
http://www.forcal.net/xiazai/forcal9/openfc32w.rar
$ Y1 c0 Z' I' w. M, `8 v# y
- I: z! k5 ^% i [% {
注意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
" M/ E; Y# n/ c$ \$ h
如能在函数内部设计上下点功夫,例子2超越matlab也是可能的。在这方面,期待高手们的指点。
7 x* b2 l& j5 M- I
1 A) d" E* w" `$ b6 J5 ]) i& i. \
如果例子2速度也超越了matlab ,matlab矢量化的神秘面纱就揭开了。
; P" E* K! O% Y9 S; ] S
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
8 R8 Z* h# o/ f# u
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
: S/ R6 U, W: S a
关于最快速的矩阵乘实现的讨论。
' p7 Y- d2 p& ~' F! k) B- N$ c9 M2 X
欢迎光临 数学建模社区-数学中国 (http://www.madio.net/)
Powered by Discuz! X2.5