20wangz 发表于 2009-6-19 06:28

这段代码如何矢量化?

已知:x=f(a,b),y=g(a,b),且有a>b>1。每当给定一个a,如果b能取一个合适的值,就能得到一个x的最优值x_opt,相应地,这时的b值称为b_x_opt。自然,知道了a和b_x_opt,也就得到了一个相应的y_x_opt。现在,我要得到x_opt-a的图像、b_x_opt-a的图像和y_x_opt-a的图像。

我现在的代码如下:(具体的函数关系就不写了)
ii = 0 ;
for a = 1 : 40
    ii = ii + 1 ;
    jj = 0 ;
    for b = 1 : 0.01 : a
      jj = jj +1 ;
      x ( jj ) = f ( a , b ) ;
      y ( jj )= g ( a , b ) ;
      s1 = struct ( 'field1' , b , 'field2' , x ( jj ) , 'field3' , y ( jj) ) ;
    end
    [ x_opt ( ii ) , II ] = max ();
    b_x_opt ( ii ) = s1 ( II ).field1 ;
    y_x_opt ( ii ) = s1 ( II ).field3 ;
end
figure (1)
plot ( 1 : 40 ,x_opt )
figure ( 2 )
plot ( 1 : 40 , b_x_opt )
figure ( 3 )
plot ( 1 : 40 , y_x_opt )

这段代码似乎也能运算出结果来。但是,这段代码全是用for循环完成的,效率很低,很长时间才能得到结果。请问各位高手,如何将这段代码矢量化呢?或者,干脆另起炉灶,来实现我的要求,该如何做呢?

先谢谢各位了!

friendchj 发表于 2009-6-20 09:49

回复 楼主 20wangz 的帖子

如果函数f,g的参数a,b允许是矢量,并且向结构体赋值时,可以以矢量的形式,可以把内层的for去掉。

20wangz 发表于 2009-6-20 14:17

函数f、g都是定义域和值域为实数的方程,只是方程很长,计算复杂些。

请恕在下愚鲁,您能不能说得详细些?最好能把代码贴出来让在下学习学习,至于函数f、g,你随便写上两个就行了。

friendchj 发表于 2009-6-20 16:54

回复 板凳 20wangz 的帖子

function shiyan
ii = 0 ;
for a = 1 : 40
    ii = ii + 1 ;
%   jj = 0 ;
%   for b = 1 : 0.01 : a
%         jj = jj +1 ;
%         x ( jj ) = f ( a , b ) ;
%         y ( jj )= g ( a , b ) ;
%         s1 = struct ( 'field1' , b , 'field2' , x ( jj ) , 'field3' , y ( jj) ) ;
%   end
%   [ x_opt ( ii ) , II ] = max ();
%   b_x_opt ( ii ) = s1 ( II ).field1 ;
%   y_x_opt ( ii ) = s1 ( II ).field3 ;
      b = 1 : 0.01 : a;
      x=f(a,b);
      y=g(a,b);
      =max(x);
      b_x_opt ( ii ) = b(id);
      y_x_opt ( ii ) = y(id) ;
end
figure (1)
plot ( 1 : 40 ,x_opt )
figure ( 2 )
plot ( 1 : 40 , b_x_opt )
figure ( 3 )
plot ( 1 : 40 , y_x_opt )

function y=f(a,b)
y=a.^2+b;

function y=g(a,b)
y=a+b.^2;

%%%%%%%
注:1.红色部分若为a^2,b^2,a,b就需要使用标量
       2.程序不用结构体可能会好点,那个结构体可以用一个3行多列标量替代。

[ 本帖最后由 friendchj 于 2009-6-20 17:00 编辑 ]

20wangz 发表于 2009-6-20 20:22

解决了!太感谢了!

困扰了我很久的一个麻烦事终于解决了!问题的关键还是我对矩阵的理解不够深刻,对Matlab的矩阵操作还是不太熟悉。

当然,要是能完全不用循环就更好了。我自己再想想办法。

rocwoods 发表于 2009-6-20 21:19

在simwe给你回了。
楼主第二重循环是b按0.01步长从1增长到a来求一系列f,g值,然后再用max求极大值。
这一来求出来的极大值不精确,二来效率低。建议楼主用fminbnd来求最值。具体看看fminbnd的帮助。
第一重循环是对不同参数的循环,这步是必须的。

20wangz 发表于 2009-6-21 21:52

用Friendchj兄的方法已经能解决我的问题了。和我原来的方法相比较,得到的结果是相同的,但是运行时间只有我的1/20左右,当把b的步长取得更小时,最短的运行时间只有我的1/80左右。太爽了!非常感谢!

另外,Rocwoods兄热心地在这里和simwe都回了帖,也非常感谢!对我的问题来说,b的步长取0.01在工程中已经足够精确了。如果还需要更精确,还可以把b的步长取得更小些。再者,从物理意义上,我已经知道每给定一个a,象这样求取b的最优值是可行的,不会因为b的有限步长而求得不适当的b的最优值。

[ 本帖最后由 20wangz 于 2009-6-21 22:10 编辑 ]

rocwoods 发表于 2009-6-21 23:06

楼主你原来程序慢主要就是两方面造成的,第一,没有预分配内存。你在
for b = 1 : 0.01 : a
         jj = jj +1 ;
         x ( jj ) = f ( a , b ) ;
         y ( jj )= g ( a , b ) ;
         s1 = struct ( 'field1' , b , 'field2' , x ( jj ) , 'field3' , y ( jj) ) ;
    end
这段代码前,也就是第二重循环开始前加两句:x = zeros(size( 1 : 0.01 : a));x = zeros(size( 1 : 0.01 : a));
只做这些变化,我敢说你的代码速度就提升不少。当然这只是针对你的f,g不能接受向量输入的情况下。
根据你的回复来看,你的f,g可以改成接受向量输入的形式,所以前面预分配内存就不用做了,因为f,g返回的是向量。
还有你之前代码慢,除了没有预分配内存外,还有一点就是频繁调用f,g函数,调用函数产生的开销很不少。friendchj兄给你改写成向量后,实际上你第二重循环原来分别调用length(1:0.01:a)次f,g函数完成的工作,只需调用一次就可以了。这部分节省的时间也比较可观。
所以按friendchj兄做法可以给你程序提高几十倍的速度。
MATLAB编程向量化的一个原则通俗点说就是"少拿多取",而不是"勤拿少取".这样可以避免频繁调用函数产生的开销。对这个原则一个不错的例子就是
http://forum.vibunion.com/thread-83299-1-2.html
下面这个帖子。

[ 本帖最后由 rocwoods 于 2009-6-21 23:09 编辑 ]
页: [1]
查看完整版本: 这段代码如何矢量化?