查看: 5252|回复: 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# a `\" J* |3 i\" w `
clc
: O8 B0 A; w) ^& e/ a
tic5 x: K0 } m' q. t5 g
k = zeros(5,5); % //生成5×5全0矩阵
\" L\" V9 x8 `: [3 p+ _- h( }
% 循环计算以下程序段1000 00次：
% P1 M* e1 A7 I( O0 @7 {\" L
for m = 1:1000 004 \ h4 o1 o6 P% t, z
a = rand(5,7);
) e+ I$ l c2 j! c' V( H2 a
b = rand(7,5);%//生成5×7矩阵a,7×5矩阵b,用0～1之间的随机数初始化
/ n+ {! D' w% e( R
k = k + a * b + a(1:5, 2:6) * b(2:6, 1:5) - a(:, 7) * b(3, :);
' ^/ G5 Y( A6 y! p$ {. j
end0 l. o9 i8 m- T: m4 S
k# V# X. |9 H* }& ]/ h& G c
toc

复制代码

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

!using["math","sys"]; [\\" Q9 p- u\\" X8 a' k% _8 L
mvar:# P8 Q( K7 |! F$ p1 O! g
t0=clock(),/ B- t2 ^/ X9 M
oo{k=zeros[5,5]},
- _% M9 @$ y; U5 N Y+ E) Z
i=0,((i++)<100000).while{
& Y7 [- A) ]/ W4 k
oo{
4 T: A: B\\" z! R J0 J3 ~3 Z5 W
a=rand[5,7], b=rand[7,5],
5 W# s! p' V& |& ^7 ?4 y- Y! 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)]
. F5 b- y; A; h% `+ Y5 h\\" t# o
}
! ~, \; C; |/ `. a1 F
},6 }( u$ D, ^% ]* C# [
k.outm(),
5 \. Z5 ~; _# A! H8 }% k9 _
[clock()-t0]/1000;

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

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

!using["math","sys"];
; r* S8 L4 G% v+ }' e
(:t0,k,i,a,b)=; |. u1 V* ]6 l; E; h$ U
{. Z# e x- e3 X' _# g
t0=clock(),
/ B( [# V( |* I. Z# J; q
oo{k=zeros[5,5]},, Q4 w. G: A' Y, S/ m
i=0,((i++)<100000).while{2 \6 j9 ~! ?1 `\\" C1 P }& `% T
oo{! c1 `. O) J6 o1 }& G9 H! e, @
a=rand[5,7], b=rand[7,5],! k% |0 T+ m, A5 x4 p0 ?2 k. m
k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)-a(neg:6)*b(3:neg)
- p V1 [! O& e4 Z j; {6 j( X
} w; t+ p$ R6 N: N# ~- u
},
1 X' N9 z/ [4 b* t\\" k7 _' `
k.outm(),
! u. X1 I# Z5 g5 L& K
[clock()-t0]/1000
5 U! Z+ N* }6 `- a3 v
};

!using["math","sys"]; 
; r* S8 L4 G% v+ }' e
(:t0,k,i,a,b)=; |. u1 V* ]6 l; E; h$ U

{. Z# e  x- e3 X' _# g

t0=clock(), 
/ B( [# V( |* I. Z# J; q
  oo{k=zeros[5,5]},, Q4 w. G: A' Y, S/ m

i=0,((i++)<100000).while{2 \6 j9 ~! ?1 `\\" C1 P  }& `% T

oo{! c1 `. O) J6 o1 }& G9 H! e, @

a=rand[5,7], b=rand[7,5],! k% |0 T+ m, A5 x4 p0 ?2 k. m

k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)-a(neg:6)*b(3:neg) 
- p  V1 [! O& e4 Z  j; {6 j( X
    }  w; t+ p$ R6 N: N# ~- u

}, 
1 X' N9 z/ [4 b* t\\" k7 _' `
  k.outm(), 
! u. X1 I# Z5 g5 L& K
  [clock()-t0]/1000 
5 U! Z+ N* }6 `- a3 v
};

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

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

//用C++代码描述为：
4 ^& V3 @, w4 h
s=0.0; : [3 g; X2 z7 c7 N3 h
for(x=0.0;x<=1.0;x=x+0.0011)
$ ]. t. L9 H* X/ ?) n
{
4 @$ A# w3 N; E( ~
for(y=1.0;y<=2.0;y=y+0.0009)
/ G& e; Y0 S8 p3 ]
{
7 F- t; G( n: b0 C- F; t
s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));$ |5 M g: \9 g' Q# P& l% a! R* S
}
5 S: Y6 ?+ u9 x- |* t5 \, l F\" u
}

复制代码

结果：
1008606.64947441
0.609 //时间

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

%file speedtest.m! b. j1 T1 C3 S! ^# t
function speedtest
5 E+ u0 v1 L9 j& [- i! L& C4 p
format long; a8 y\" ^. f* G2 Y$ f6 q$ E
tic
# }5 y7 c2 Q o: U
[x,y]=meshgrid(0:0.0011:1,1:0.0009:2);\" ^8 f/ v8 K+ z7 n6 E' K& \
s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))' o- y1 d7 f. k7 [9 i7 R6 }: a0 p& Y
toc

复制代码

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

!using["math","sys"];
) s$ U% L! G& h' K
mvar:% e; l* I Q& c- D
t=clock(),& C8 b2 m+ W' x2 Z6 C
oo{1 H* P. _6 p5 A- x2 h6 T/ b; _
ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
! j2 l) z4 F F2 i& p( t6 M* 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], G* q9 j7 ^2 ^8 x\\" E
};3 ~- }6 G: ]3 p( i0 L
[clock()-t]/1000;

结果：
1008606.64947441
0.625 //时间

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

!using["math","sys"];
+ e, {/ U% ?) n/ t2 U& s
mvar:$ P. C' x% K5 y& f! A
t=clock(),5 U3 F) S/ ^# p
oo{; s1 |- ?! w& X3 d$ @
ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
1 G* u I/ r1 E) B
Sum[Sum[Cos(rn(1)-Sin(rn(1.2)*x^(y/rn(2))+Cos(rn(1)-Sin(rn(1.2)*y^(x/rn(2))))))]]% O\\" `% N% b5 @. ?
};
9 M7 ]% ]! p0 r# |
[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)))));
]5 {7 T7 v$ K6 Z6 o
sum["f",0,1,0.0011,1,2,0.0009];

复制代码

结果：
1008606.64947441
0.719 //时间

Forcal代码3：while循环

mvar:
( p5 d& Q2 \, Z1 l: k, d
t=sys::clock();- p& o, `+ \+ K% @; i: W
s=0,x=0,
- C% `; F3 q( @
while{x<=1, //while循环算法； * E; \, w7 f. }4 ^3 H
y=1, 9 _ O. l6 K% w) }\" c6 o7 M
while{y<=2, ) V; o7 Q9 |! R4 S' \6 ?
s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))), # r8 i5 {! z- L
y=y+0.0009 / ^4 m& O) |; ~/ T
}, 8 u$ A' ]$ R8 l: |* t0 A
x=x+0.0011
0 ]7 e( B6 m0 M0 q4 u: g
}, 6 F; {, d$ E\" I* t
s;
' M6 F! k' N\" h. Y5 F* Y
[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++还要高些。