查看: 5167|回复: 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
0 d3 }; v, g+ n' u: @6 e1 e8 p! U
clc7 k: ^, c9 ^\" b8 ^# k* k7 N
tic
' Q7 R6 Z5 c4 z1 d! ]
k = zeros(5,5); % //生成5×5全0矩阵; r6 t! O; \; P* ~* o/ v9 E
% 循环计算以下程序段1000 00次：
v5 Y- r3 [& t4 J8 g! B- d
for m = 1:1000 00
# K4 ?1 ]3 g. S1 k F\" V& \
a = rand(5,7);
1 _/ ?0 ~9 N7 r\" t
b = rand(7,5);%//生成5×7矩阵a,7×5矩阵b,用0～1之间的随机数初始化
; W+ ^0 K( |) w\" x6 T* T( y) r
k = k + a * b + a(1:5, 2:6) * b(2:6, 1:5) - a(:, 7) * b(3, :);: o7 f* H0 w+ ~$ O1 n
end
- ^5 @) h! Q6 c0 m
k. K2 g( O4 ?6 @: I& W
toc

复制代码

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

!using["math","sys"];
7 W3 P% |+ P0 K p3 O+ M
mvar:0 {- o; g8 m4 z: t& E$ A, P5 T! A
t0=clock(),
% i- f! t+ k. K7 {3 D# j
oo{k=zeros[5,5]},! s: ` V' J* ~9 |% t
i=0,((i++)<100000).while{6 K& N. r. ?# P0 n
oo{9 L$ F) u' G! |2 [2 g2 {; [, N
a=rand[5,7], b=rand[7,5],
; ~5 D) z' d/ K8 r* B3 v5 {
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 C- B( G0 ]+ i+ P+ P- g2 r
}2 ]3 g6 X\\" R3 k4 y. b
},
8 B2 W7 Q6 Q\\" f6 ?; z
k.outm(),9 Q0 i! P9 j, w* x3 Q/ x2 {! `
[clock()-t0]/1000;

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

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

!using["math","sys"];
( |5 a, h3 U% M: ?3 K) N; @- L
(:t0,k,i,a,b)=! ~6 A' z; Z* M5 E {
{
3 E# _8 @2 q: @ i. C
t0=clock(),$ W6 u5 c9 c/ B- p
oo{k=zeros[5,5]},
8 W4 x5 ]0 |0 C\\" Q* P1 X
i=0,((i++)<100000).while{
; q! @6 X8 D0 l3 p3 K\\" o
oo{) K8 u- c8 p( f
a=rand[5,7], b=rand[7,5],+ Y6 }( g# z$ J7 F; m
k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)-a(neg:6)*b(3:neg)( M# B0 {' i4 j
}% a\\" X6 N; F$ J) u3 y
},. y) k' R, q1 S# |; ^
k.outm(),
, ^/ g# ]2 H5 H
[clock()-t0]/10006 L' D- d/ _* S( E9 v
};

!using["math","sys"]; 
( |5 a, h3 U% M: ?3 K) N; @- L
(:t0,k,i,a,b)=! ~6 A' z; Z* M5 E  {

{ 
3 E# _8 @2 q: @  i. C
  t0=clock(),$ W6 u5 c9 c/ B- p

oo{k=zeros[5,5]}, 
8 W4 x5 ]0 |0 C\\" Q* P1 X
  i=0,((i++)<100000).while{ 
; q! @6 X8 D0 l3 p3 K\\" o
    oo{) K8 u- c8 p( f

a=rand[5,7], b=rand[7,5],+ Y6 }( g# z$ J7 F; m

k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)-a(neg:6)*b(3:neg)( M# B0 {' i4 j

}% a\\" X6 N; F$ J) u3 y

},. y) k' R, q1 S# |; ^

k.outm(), 
, ^/ g# ]2 H5 H
  [clock()-t0]/10006 L' D- d/ _* S( E9 v

};

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

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

//用C++代码描述为：
( @2 \4 v* n! r
s=0.0;
3 R! o' J: T o2 M5 M4 @
for(x=0.0;x<=1.0;x=x+0.0011)
8 B4 |$ x- t K( c9 R& B
{
2 K& P) B5 E% {# J+ n) y
for(y=1.0;y<=2.0;y=y+0.0009)
) V9 t: k' |/ N% l; V\" B
{9 L, R/ g5 [& K: o, y2 T: I0 f
s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
7 Y9 h0 x* k* R6 ]% D' A( B I7 e3 z
}; x# \+ U, B' x: X2 B) q: m
}

复制代码

结果：
1008606.64947441
0.609 //时间

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

%file speedtest.m
& ~3 k\" s% k7 P/ a' B
function speedtest
4 R, s0 _( N% E9 [% N( y
format long
; V2 _. K3 }) @
tic
+ {\" E* }7 i5 Z1 V. T
[x,y]=meshgrid(0:0.0011:1,1:0.0009:2);8 P6 R# m, P6 ]4 \. u
s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))) r# @9 z: w8 b* f* h
toc

复制代码

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

!using["math","sys"];
& s4 o+ u- ^8 N: S
mvar:$ u2 K% q7 u8 @3 j' Y. l
t=clock(),
0 j6 W7 I. ~! p& c7 S
oo{
- H- c7 f) B6 K) B1 T9 P. p$ ?
ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
/ k& K4 z/ X6 ?$ n3 M
Sum[Cos(rn(1)-Sin(rn(1.2)*x^(y/rn(2))+Cos(rn(1)-Sin(rn(1.2)*y^(x/rn(2)))))),0]
8 V/ T% A' D V* q+ ?/ J
};
% ]: V9 c2 @7 L
[clock()-t]/1000;

结果：
1008606.64947441
0.625 //时间

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

!using["math","sys"];( N) c6 F5 j8 d: ?. u' F
mvar:
; @9 f! Q8 A' [7 p\\" W5 s
t=clock(),
2 j4 W! _% }+ h7 G8 X
oo{
2 m9 C7 y% d( ]% S\\" `
ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],) b1 K: v V\\" ^
Sum[Sum[Cos(rn(1)-Sin(rn(1.2)*x^(y/rn(2))+Cos(rn(1)-Sin(rn(1.2)*y^(x/rn(2))))))]]- s: E+ T! [1 @9 g Q7 U6 ~/ Q
};
8 S2 J) |/ j2 T6 P# h
[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)))));
\" }/ a* Q3 @9 K
sum["f",0,1,0.0011,1,2,0.0009];

复制代码

结果：
1008606.64947441
0.719 //时间

Forcal代码3：while循环

mvar:
. D. O) Z4 i4 y( o4 W w4 g5 N0 p
t=sys::clock();
+ }( ^( X8 R# W. r\" X; K
s=0,x=0, & L1 A2 C& X: E1 g1 a
while{x<=1, //while循环算法； ; z/ \; u& z0 A3 d. r
y=1,
% z! W3 j5 |! t
while{y<=2,
s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
# h% S% w\" }. v3 K8 M
y=y+0.0009
2 o& Y& r1 O$ s( F. n: W
}, 5 f1 I( T\" a# J# G+ P
x=x+0.0011
$ U\" Z0 ~% |5 u$ p: \! c
},
6 i. i) |9 T0 u% o& e
s;. p+ z# I& a7 K0 |# Q# g
[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++还要高些。