查看: 5165|回复: 8

关于matlab代码矢量化的理解

[复制链接]

字体大小: 正常放大

forcal

45 主题	3 听众	282 积分

升级 91%

TA的每日心情

	难过 2012-8-27 18:22

签到天数: 1 天

[LV.1]初来乍到

电梯直达

1^#

发表于 2010-10-5 09:32 |只看该作者 |倒序浏览

|招呼Ta 关注Ta

代码矢量化是matlab的精髓，其基本特点是运行速度快和代码简洁，它是如何实现的？

按我的理解，代码矢量化的本质就是设计专门的函数对数组元素集中运算，这样可提高运行速度，同时兼有代码简洁的特点。

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

大家有什么看法，愿畅所欲言。

zan

转播0 淘帖0 分享0 收藏0 支持0 反对0 微信

使用道具举报

qbist

2 主题	3 听众	304 积分

升级 1.33%

该用户从未签到

自我介绍: 一个对未来充满信心的阳光型男孩！

2^#

发表于 2010-10-5 09:52 |只看该作者 |招呼Ta 关注Ta

我正在学 Matlab 软件！~~

我也说一句

发表

使用道具举报

zhwqqiangge

9 主题	3 听众	186 积分

升级 43%

TA的每日心情

	开心 2011-9-11 13:24

签到天数: 1 天

[LV.1]初来乍到

自我介绍: 大家好!我是新手!请多多关照!

3^#

发表于 2010-10-5 10:17 |只看该作者 |招呼Ta 关注Ta

??????????????????

我也说一句

发表

使用道具举报

wznzy0822

4 主题	3 听众	846 积分

升级 61.5%

TA的每日心情

	开心 2012-12-6 22:41

签到天数: 113 天

[LV.6]常住居民II

群组: 数学建模

群组: 数学建模培训课堂2

4^#

发表于 2010-10-5 10:52 |只看该作者 |招呼Ta 关注Ta

顶。。。。。。。。。。。。。。。。。。。。。。

我也说一句

发表

使用道具举报

master_math

19 主题	4 听众	235 积分

升级 67.5%

TA的每日心情

	开心 2016-12-19 06:10

签到天数: 32 天

[LV.5]常住居民I

群组: 数学建摸协会

5^#

发表于 2010-10-5 18:32 |只看该作者 |招呼Ta 关注Ta

,,,顶顶更健康！

我也说一句

发表

使用道具举报

forcal

45 主题	3 听众	282 积分

升级 91%

TA的每日心情

	难过 2012-8-27 18:22

签到天数: 1 天

[LV.1]初来乍到

6^#

发表于 2010-10-12 21:36 |只看该作者 |招呼Ta 关注Ta

本帖最后由 forcal 于 2010-10-12 21:46 编辑

我正在练手设计的FcMath库也打算以矩阵运算为基础，设计一些专门的函数对数组元素集中运算，运行效率确实有所提高（甚至有些涉及矩阵的算法比matlab还快），代码也简洁了，但不知这是不是矢量化？

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

我感觉，matlab的函数调度效率较低，对象管理效率这个不好说，但一些函数内部的设计比较优秀。故有些Forcal代码比matlab快，而有些慢。

以下例子体现了Forcal和matlab的效率差别所在。

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

clear all
$ T5 M5 D/ o0 d0 P8 I
clc; K2 W& L( L, v4 W; K
tic
# U, Q j\" O1 t5 o$ |& |; Z
k = zeros(5,5); % //生成5×5全0矩阵2 d+ ]( b% s3 |( s# ?5 _/ ]
% 循环计算以下程序段1000 00次：
, s0 v( J1 l- G m4 z$ E3 @5 _& P
for m = 1:1000 00
& ^2 ]3 Q; G' ~* G# Z4 Z
a = rand(5,7);$ y2 V) p) v! Y) F9 t
b = rand(7,5);%//生成5×7矩阵a,7×5矩阵b,用0～1之间的随机数初始化+ U& D$ W- z+ ]2 D
k = k + a * b + a(1:5, 2:6) * b(2:6, 1:5) - a(:, 7) * b(3, :);
! I3 D- f) y/ \' H
end# v# V\" ]/ J! ?% }* x: A6 s# R6 ^) ?
k
9 R2 V6 Q, y2 X7 t
toc

复制代码

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

!using["math","sys"];' ]% c% v- e: K3 M
mvar:
- J* I5 `/ V: h% B8 P
t0=clock(),! }2 ^7 b4 F% t+ P/ x9 P6 {
oo{k=zeros[5,5]},
) T: v7 q [ ?& q% M\\" p0 a
i=0,((i++)<100000).while{2 q0 A- y# f/ V: N0 u2 K
oo{
\\" y! R. ~5 I8 @1 S1 X
a=rand[5,7], b=rand[7,5],: l$ d' z1 ?- f9 \7 k3 P; z
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)]; p3 k/ a& b. b7 b5 B- ?* e( v2 j
}7 N( s- t) x+ o! Y/ o# P( O0 X' s
},
5 T, W3 o* m# {+ _$ G6 s0 m) }
k.outm(),
: K8 B: C4 q3 R* o& I
[clock()-t0]/1000;

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

Forcal代码2：比较好看些的代码，似乎也比matlab稍快吧？

!using["math","sys"];
9 P6 L% H; \( u% O8 J: W8 v
(:t0,k,i,a,b)=, l, z& I6 D2 [$ \8 B4 }# U
{
( H& ~4 d, Y( R' S
t0=clock(),2 j5 L6 p6 [. s' b- X. E
oo{k=zeros[5,5]},
- x* j1 S% `' ^9 H/ j8 M. n
i=0,((i++)<100000).while{
0 f4 F3 y) z9 W# m4 J
oo{7 `* _* {2 e+ \; a0 I
a=rand[5,7], b=rand[7,5],, y: |1 b, N! F$ y
k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)-a(neg:6)*b(3:neg)
K. K5 c: @& B, M& X. V0 q3 s
}! _ Z0 z' j; E+ D; g
},; e9 h) S/ x0 O& {) x1 y) I
k.outm(),0 h, M& o/ ~- \: r# O0 J
[clock()-t0]/1000# R( T0 W3 v3 C$ Z
};

!using["math","sys"]; 
9 P6 L% H; \( u% O8 J: W8 v
(:t0,k,i,a,b)=, l, z& I6 D2 [$ \8 B4 }# U

{ 
( H& ~4 d, Y( R' S
  t0=clock(),2 j5 L6 p6 [. s' b- X. E

oo{k=zeros[5,5]}, 
- x* j1 S% `' ^9 H/ j8 M. n
  i=0,((i++)<100000).while{ 
0 f4 F3 y) z9 W# m4 J
    oo{7 `* _* {2 e+ \; a0 I

a=rand[5,7], b=rand[7,5],, y: |1 b, N! F$ y

k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)-a(neg:6)*b(3:neg) 
  K. K5 c: @& B, M& X. V0 q3 s
    }! _  Z0 z' j; E+ D; g

},; e9 h) S/ x0 O& {) x1 y) I

k.outm(),0 h, M& o/ ~- \: r# O0 J

[clock()-t0]/1000# R( T0 W3 v3 C$ Z

};

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

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

//用C++代码描述为：9 W/ e. w2 _; s1 e5 A0 B% ^, v
s=0.0;
9 b7 z \# B3 d+ K
for(x=0.0;x<=1.0;x=x+0.0011)
6 q) H. p' E1 f% L\" t) L
{8 ?4 W. k5 {) @: J9 V6 Z$ @\" ?( b# i
for(y=1.0;y<=2.0;y=y+0.0009). E2 v* h& @$ y L- G& A2 _. r
{
5 i\" }' W8 o$ M9 w# l
s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));8 q, w( D& ^: _0 c9 {6 }
}
% Q: V& `5 o( y6 e7 w1 h2 ^8 M) Y; K
}

复制代码

结果：
1008606.64947441
0.609 //时间

这个matlab程序段是网友yycs001给出的。

%file speedtest.m) I8 z' A\" Z\" k3 i4 y# c\" y3 p
function speedtest
[9 j! @+ @( A
format long
3 `8 u* v4 S: x' ]) M5 L, r
tic
$ J+ L6 M6 n0 d. X2 C
[x,y]=meshgrid(0:0.0011:1,1:0.0009:2);1 G& s( R( j2 }3 [) [$ w
s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))& c. D6 T( k0 `
toc

复制代码

Forcal代码1：**数组求和函数Sum，完全矢量化的代码

!using["math","sys"];5 C1 K3 W; ^+ g; o9 i) W
mvar:: |1 l3 `$ V' u7 a0 w
t=clock(),
! L8 [& e V\\" ]# g* a* \. e. ~
oo{
% q) r+ j7 s3 d# y
ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
( u0 o6 Q# v) |2 d6 K
Sum[Cos(rn(1)-Sin(rn(1.2)*x^(y/rn(2))+Cos(rn(1)-Sin(rn(1.2)*y^(x/rn(2)))))),0]
/ A1 G+ ^& y) v4 ?1 e
};1 n# q+ T4 l% a, M2 E0 @
[clock()-t]/1000;

结果：
1008606.64947441
0.625 //时间

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

!using["math","sys"];
0 a: b9 X0 ?0 [3 V, T\\" @$ {- C
mvar:
K/ ]8 e: a% K1 [1 {
t=clock(),* i4 D% _3 w: y/ a
oo{
8 b; i+ G g! {/ R7 z1 E; x8 D
ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],) C4 k\\" @3 _/ x( N
Sum[Sum[Cos(rn(1)-Sin(rn(1.2)*x^(y/rn(2))+Cos(rn(1)-Sin(rn(1.2)*y^(x/rn(2))))))]]& N6 V5 q, M- t$ f1 d. o
};
5 u- [$ X' q4 c4 J* l, F7 w. K
[clock()-t]/1000;

Forcal代码2：求和函数sum，非矢量化代码

f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); , G; F& B/ A7 W
sum["f",0,1,0.0011,1,2,0.0009];

复制代码

结果：
1008606.64947441
0.719 //时间

Forcal代码3：while循环

mvar:' x# _5 V; _\" ]5 h2 |% k
t=sys::clock();+ ]% o Q: M6 q/ j7 A5 v
s=0,x=0,
$ s0 ?5 \# I\" [5 K- Q9 h* c8 z: `
while{x<=1, //while循环算法； p9 D+ M# o8 g1 Q
y=1,
' ]! p* y9 C, F) L- Z( A! t
while{y<=2, 5 S. J! [7 v3 c
s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))), : O) J5 C0 t1 U6 e
y=y+0.0009
* X\" k% p8 k8 m% g5 z$ M
}, ) ?8 O$ Q% q- E
x=x+0.0011 6 n/ _/ E7 c) L& h* U3 P$ g/ K
},
* o\" T8 g7 ~, L- Z* n' B& y4 ?
s;2 ^/ u+ I ~% e+ Y$ C
[sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
0.734 //时间

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

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

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

例子2中Forcal的矢量化代码是最快的，但与matlab的矢量化代码相比仍有差距。该例子的特点是函数调用少，临时变量也少，但矩阵大。故说明Forcal的各种矩阵函数Sin、Cos及矩阵的加减运算等函数的内部设计不及matlab。

如能在函数内部设计上下点功夫，例子2超越matlab也是可能的。在这方面，期待高手们的指点。

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

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

我也说一句

发表

使用道具举报

forcal

45 主题	3 听众	282 积分

升级 91%

TA的每日心情

	难过 2012-8-27 18:22

签到天数: 1 天

[LV.1]初来乍到

7^#

发表于 2010-10-16 16:59 |只看该作者 |招呼Ta 关注Ta

讨论有益！以下是在其他论坛的讨论帖子：
simwe：http://forum.simwe.com/thread-952532-1-3.html
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

我也说一句

发表

使用道具举报

forcal

45 主题	3 听众	282 积分

升级 91%

TA的每日心情

	难过 2012-8-27 18:22

签到天数: 1 天

[LV.1]初来乍到

8^#

发表于 2010-10-24 16:59 |只看该作者 |招呼Ta 关注Ta

参考：http://bbs.emath.ac.cn/thread-2709-1-1.html

我在多个帖子中有不同说明，但将这些说明再集中到一个帖子中比较麻烦，大家相互参考一下，看能否把这个问题解决了。

我也说一句

发表

使用道具举报

forcal

45 主题	3 听众	282 积分

升级 91%

TA的每日心情

	难过 2012-8-27 18:22

签到天数: 1 天

[LV.1]初来乍到

9^#

发表于 2010-11-2 18:06 |只看该作者 |招呼Ta 关注Ta

参考：http://bbs.emath.ac.cn/thread-2727-1-1.html

关于最快速的矩阵乘实现的讨论。

我也说一句

发表

使用道具举报

返回列表

帐号		密码		只需要一步，快速开始		注册地址	找回密码

关于matlab代码矢量化的理解

浏览过的版块

新人进步奖

QQ

电话咨询

关于我们| 联系我们| 诚征英才| 对外合作| 产品服务|