这段代码如何矢量化？

20wangz

已知：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 ([s1.field2]);
    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

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

20wangz

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

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

friendchj

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 ([s1.field2]);
%     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);
      [x_opt( ii ),id]=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

解决了！太感谢了！

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

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

rocwoods

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

20wangz

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

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

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

rocwoods

楼主你原来程序慢主要就是两方面造成的，第一，没有预分配内存。你在

1. 
   
2. for b = 1 : 0.01 : a
   
3.          jj = jj +1 ;
   
4.          x ( jj ) = f ( a , b ) ;
   
5.          y ( jj )= g ( a , b ) ;
   
6.          s1 = struct ( 'field1' , b , 'field2' , x ( jj ) , 'field3' , y ( jj) ) ;
   
7.     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 编辑 ]