查看: 5005|回复: 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 微信

使用道具举报

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

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

我也说一句

发表

使用道具举报

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]初来乍到

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]初来乍到

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
2 W* `+ n4 Q% i! c; w
clc
. | ?+ `! J: A, y# e: `5 w
tic
. L2 i' E3 M4 X3 \4 }
k = zeros(5,5); % //生成5×5全0矩阵
( _4 F: B, U3 c2 I2 ^
% 循环计算以下程序段1000 00次：6 L S1 p3 g5 I- U* D' Y
for m = 1:1000 00
, I2 _& X2 r# ~! F, U+ f, B
a = rand(5,7);3 a; ?9 r' k! l9 M0 ]6 J1 F7 p( X
b = rand(7,5);%//生成5×7矩阵a,7×5矩阵b,用0～1之间的随机数初始化
; a ^\" n1 T# E+ F) S
k = k + a * b + a(1:5, 2:6) * b(2:6, 1:5) - a(:, 7) * b(3, :);/ Q( F4 T e5 q. ^/ H
end9 l- u9 F9 X' j5 A
k
9 U2 Z7 v6 ^$ q7 x+ K- a
toc

复制代码

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

!using["math","sys"];! h: V* k5 P# n y. H\\" T
mvar:( o. V9 V) f4 k! O
t0=clock(),7 {1 F* C/ |6 K
oo{k=zeros[5,5]},) V( [5 _/ _, [+ x- d8 t* H
i=0,((i++)<100000).while{5 [0 c1 U1 G2 X' p
oo{+ E; }2 M) R d0 O9 Z5 I) Q) K1 |
a=rand[5,7], b=rand[7,5],0 o2 L7 Y( G: s# e. ?4 W( j7 J
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+ W/ F- c, X f- y2 ]
}( J7 ^3 T, c H: _& O
},; ]1 f, ]0 H5 O( Y6 d5 N
k.outm(),. B! L5 }' z( F( r5 y8 t, l
[clock()-t0]/1000;

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

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

!using["math","sys"];( L! _# j+ P) u. n: G& U1 v
(:t0,k,i,a,b)=+ J, @0 i% Y8 [7 c
{
r& n4 a* o: P3 ]+ r& i1 c# y+ ]6 Q
t0=clock(),
7 I; J W( v' l7 ^* p; t2 S
oo{k=zeros[5,5]},
- n% L* z8 d) U- X! c7 ^
i=0,((i++)<100000).while{. p2 Q& Z: [4 L7 ^
oo{
* T: c! N4 c/ F* d* X, z
a=rand[5,7], b=rand[7,5],
r' |3 p5 R0 J( K1 ]# c
k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)-a(neg:6)*b(3:neg): j2 b$ k! P: }& A; b4 w0 U
}, ~, i( M/ W. Q- Y& Z: i5 t( ?2 \% e
},4 K$ h8 Y6 |& q5 x2 {. ^
k.outm(),
5 ]& I6 s& {+ q# N) `5 w\\" h; G
[clock()-t0]/1000
$ ~: K7 M8 A2 s
};

!using["math","sys"];( L! _# j+ P) u. n: G& U1 v

(:t0,k,i,a,b)=+ J, @0 i% Y8 [7 c

{ 
  r& n4 a* o: P3 ]+ r& i1 c# y+ ]6 Q
  t0=clock(), 
7 I; J  W( v' l7 ^* p; t2 S
  oo{k=zeros[5,5]}, 
- n% L* z8 d) U- X! c7 ^
  i=0,((i++)<100000).while{. p2 Q& Z: [4 L7 ^

oo{ 
* T: c! N4 c/ F* d* X, z
      a=rand[5,7], b=rand[7,5], 
  r' |3 p5 R0 J( K1 ]# c
      k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)-a(neg:6)*b(3:neg): j2 b$ k! P: }& A; b4 w0 U

}, ~, i( M/ W. Q- Y& Z: i5 t( ?2 \% e

},4 K$ h8 Y6 |& q5 x2 {. ^

k.outm(), 
5 ]& I6 s& {+ q# N) `5 w\\" h; G
  [clock()-t0]/1000 
$ ~: K7 M8 A2 s
};

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

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

//用C++代码描述为：/ k7 k) C& i, \6 S
s=0.0;
, n% z- _3 u' q% x
for(x=0.0;x<=1.0;x=x+0.0011)
; Z# `5 {% @; h2 T K* t
{
2 H0 J$ [\" _5 }* y
for(y=1.0;y<=2.0;y=y+0.0009)9 [4 ^ A8 o9 M. p) z% ~+ ]
{
. ~5 h, r6 ]6 S\" s- o4 j
s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));7 b: W! X1 R: J9 I# R1 \
}) u7 m1 S+ ~/ @9 U+ K
}

复制代码

结果：
1008606.64947441
0.609 //时间

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

%file speedtest.m X# M5 N2 e Q$ A* r. [
function speedtest0 U: O\" h- j) I3 ^5 B: \
format long
1 `& g& b% f- {+ q7 Y2 K, h* k\" {
tic
/ |% S# s$ C D! e( ?
[x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
; j) O! R$ N9 Q( w; C
s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))( F0 r3 y* D+ }& m
toc

复制代码

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

!using["math","sys"];
1 l; |* f6 h9 E5 m
mvar:5 [: q: l1 k1 J- l& I+ g
t=clock(),* j& M/ H4 J! B4 D: _1 j% \( p' M
oo{
# D2 x8 R2 ^ ?; K% g
ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],1 Z- B' D+ W/ O5 L% q3 C f
Sum[Cos(rn(1)-Sin(rn(1.2)*x^(y/rn(2))+Cos(rn(1)-Sin(rn(1.2)*y^(x/rn(2)))))),0]5 `7 d2 H% V8 m) D5 T3 K; H4 d
};
. z; P5 a( d) n. F$ u
[clock()-t]/1000;

结果：
1008606.64947441
0.625 //时间

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

!using["math","sys"];/ W8 C3 Y6 }8 `( ?6 V; K' B
mvar:
1 c\\" B; w6 B8 W* E! f
t=clock(),
/ D0 @( a2 V& L. b\\" `9 @; Z- w
oo{
7 P+ G& C; X4 u, U2 G/ Z
ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
& U1 y. K6 r' |; _; V. g9 y* G
Sum[Sum[Cos(rn(1)-Sin(rn(1.2)*x^(y/rn(2))+Cos(rn(1)-Sin(rn(1.2)*y^(x/rn(2))))))]]; t, B8 h1 b* j' u _# h; n
};; T6 ^, U0 W( n/ D5 @
[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))))); 8 Z3 O9 r! `& b( m$ _; O+ v
sum["f",0,1,0.0011,1,2,0.0009];

复制代码

结果：
1008606.64947441
0.719 //时间

Forcal代码3：while循环

mvar:2 }( F* J l. K+ B! T
t=sys::clock();
# d. l5 O C/ a
s=0,x=0, 1 H. F; q z, w& p8 c\" o! @
while{x<=1, //while循环算法； / {; d. ?( ^1 n$ m9 `3 ^
y=1, - L\" q1 F; n* m0 e
while{y<=2,
) b! S+ ^# I$ p2 H. n% H
s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
- n( }* `1 B4 T6 V2 l
y=y+0.0009
' |, a6 c( z1 `4 @1 @
},
) T+ K6 m( u* L0 h* O; |
x=x+0.0011 % x0 F: p- Y8 y l
}, * [9 @( S+ w+ W4 `! Z6 C
s;1 H\" q! R+ j. w
[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++还要高些。