查看: 5169|回复: 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$ ?. o+ ^7 ?( n1 S7 l\" S
clc
7 L$ w; Z0 j$ b& v7 u. |, \
tic( ^4 v N$ |3 x. i! ]
k = zeros(5,5); % //生成5×5全0矩阵
: R' @4 u! C7 C\" T6 y5 f1 S$ Z
% 循环计算以下程序段1000 00次：& D8 L u( W! w& U- @( P
for m = 1:1000 008 Z. H% Y$ N1 j- d$ h7 C9 O- x; V
a = rand(5,7);4 ?+ [% O# D7 d
b = rand(7,5);%//生成5×7矩阵a,7×5矩阵b,用0～1之间的随机数初始化: `) R2 o7 n9 m/ F4 l
k = k + a * b + a(1:5, 2:6) * b(2:6, 1:5) - a(:, 7) * b(3, :);% b' R# \9 ?( P2 i, b
end
- t: c. U/ L8 b5 A ~1 q
k
% n1 c( A\" I; R4 l7 N
toc

复制代码

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

!using["math","sys"];
6 s7 r& F% h1 P( T
mvar:7 S* O ^0 a( R, w: u* t! X. M\\" ~
t0=clock(),+ k s- M* j8 C\\" t
oo{k=zeros[5,5]},* S- a; H7 m& V+ ]
i=0,((i++)<100000).while{! V3 Y0 Q# e: i+ m9 n2 B ]/ K
oo{
: R6 U F4 c' N( q0 w\\" q
a=rand[5,7], b=rand[7,5],% j$ @) p) q N: M- a6 q6 ?% |
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)]9 V9 y& q0 Q% o: {/ u' J. ?# j
}8 G) |& E4 {9 v; [1 x
},
2 k$ a$ g: }3 P) ^1 j7 g
k.outm(),- M+ b- I' X$ ~$ m! l
[clock()-t0]/1000;

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

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

!using["math","sys"];
1 h9 |+ u5 p E0 b
(:t0,k,i,a,b)=7 ]7 L s7 U, _( t& m! p% w! l5 b
{
& A0 w8 K, b; B$ W5 a' j
t0=clock(),9 t7 [9 D2 b9 @
oo{k=zeros[5,5]},
0 M+ j: E. \. u9 v( l, } s, f1 |
i=0,((i++)<100000).while{1 I: J9 y) \( `7 g/ L2 o% R
oo{ j3 j4 _* q8 G; f9 C
a=rand[5,7], b=rand[7,5],0 m$ n2 v1 K% ?, p
k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)-a(neg:6)*b(3:neg)\\" y, s# A$ B% g% h- I$ P: h
}, ^0 P) C$ }, Q! |
},3 C+ B: ~9 w\\" |2 m
k.outm(),
; R0 n: S1 a8 P w8 F7 f; G
[clock()-t0]/1000/ n A$ ?0 m; |* @* p
};

!using["math","sys"]; 
1 h9 |+ u5 p  E0 b
(:t0,k,i,a,b)=7 ]7 L  s7 U, _( t& m! p% w! l5 b

{ 
& A0 w8 K, b; B$ W5 a' j
  t0=clock(),9 t7 [9 D2 b9 @

oo{k=zeros[5,5]}, 
0 M+ j: E. \. u9 v( l, }  s, f1 |
  i=0,((i++)<100000).while{1 I: J9 y) \( `7 g/ L2 o% R

oo{  j3 j4 _* q8 G; f9 C

a=rand[5,7], b=rand[7,5],0 m$ n2 v1 K% ?, p

k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)-a(neg:6)*b(3:neg)\\" y, s# A$ B% g% h- I$ P: h

}, ^0 P) C$ }, Q! |

},3 C+ B: ~9 w\\" |2 m

k.outm(), 
; R0 n: S1 a8 P  w8 F7 f; G
  [clock()-t0]/1000/ n  A$ ?0 m; |* @* p

};

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

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

//用C++代码描述为：
4 [ o, t) p5 q7 m d' u2 z% P
s=0.0;
' H\" n3 m\" @$ }
for(x=0.0;x<=1.0;x=x+0.0011)
\" U* h. x S& U7 M\" r8 H
{
2 r7 Y. \9 b/ E$ r2 H; Z
for(y=1.0;y<=2.0;y=y+0.0009)
( \6 `# t& E' @) d8 ^- T
{
) o/ i0 x% ?- f3 |8 I0 E
s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));/ }, C3 O\" }7 C/ J8 i
}\" E3 p7 [; q* f8 K
}

复制代码

结果：
1008606.64947441
0.609 //时间

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

%file speedtest.m
, P' T+ `3 M% T1 x* x* U% B4 M
function speedtest
& X/ d% J# o3 @6 E x1 y- N) b# f$ I
format long! x( w% @+ T\" ?+ k+ M4 V
tic- y9 a3 ]3 \* k1 j; \! H
[x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
& V% T: }5 S; P4 b- X
s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2))))))): O6 e7 I1 M; R4 y9 d w7 Y3 E
toc

复制代码

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

!using["math","sys"];
4 I) ]0 j+ i) G
mvar:
- G% A\\" Q; L1 H4 z9 J% D( {# q8 U& _
t=clock(),
3 z3 O% m1 z, R+ l* \. u1 C
oo{
% u; {, m* `\\" T8 M
ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
2 q7 E$ a6 m: ?0 A% j; J. q3 \' l
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 N+ f* r7 @/ E\\" g6 w
};; B) A8 H# A! k9 G4 C
[clock()-t]/1000;

结果：
1008606.64947441
0.625 //时间

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

!using["math","sys"];
2 P! J ~- {2 @$ O\\" v' m: r' X
mvar:6 S4 O* p$ d: E* {1 U+ d
t=clock(),
% h5 w6 a: _. z
oo{
/ ]& n: W; S8 w\\" e4 \
ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
\\" J+ R$ u3 w, H$ _% B* V. i% d0 m/ M
Sum[Sum[Cos(rn(1)-Sin(rn(1.2)*x^(y/rn(2))+Cos(rn(1)-Sin(rn(1.2)*y^(x/rn(2))))))]]) n G; P* o8 r9 f* U: l' S
};3 J# M% f) A! J' s+ B% D3 |( e
[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 G V1 w' z6 [& C8 g4 R K2 c
sum["f",0,1,0.0011,1,2,0.0009];

复制代码

结果：
1008606.64947441
0.719 //时间

Forcal代码3：while循环

mvar:
) L\" ~8 R7 n9 V/ n4 H( Y& Z% H
t=sys::clock();7 D/ P. K: K7 r0 O0 W0 e\" G j
s=0,x=0, / X! r\" K$ W\" q6 P! G7 S
while{x<=1, //while循环算法； # _. @$ P# f# l3 A5 q; A
y=1,
+ U* H- i6 H) D- O) j
while{y<=2,
* g, K1 e) w; G. w$ r
s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
/ a4 ?+ A5 V) s$ A5 j
y=y+0.0009 + X X l! o; p* r' v
},
! I ~- S9 j6 z+ A9 U0 k
x=x+0.0011
2 L! j q- n6 Y5 B/ D+ r
},
% u, Y( g# B$ X8 q* e9 C2 }- ~
s;) Z\" u: r3 ~: `4 x7 U2 i
[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

电话咨询

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