求助：循环优化的问题

berryhaw

编的程序中有一个循环要调用几百万次，用profiler看了一下，发现这个循环比较耗时，贴出来大家帮忙看一下，能不能用向量化或矩阵的思想改进？
    程序作用：找出同时出现在一维数组a和一维数组b中的元素，并将b中的该元素置零。（b的长度与a的长度不同）
如：
b=[4 6 8]; a=[100 2 6 96 8 7 20];
for i=1:3
    dd=a(a==b(i));   
    if isempty(dd)==0, b(i)=0; end
end
运行后结果为：b=[4 0 0];

[ 本帖最后由 ChaChing 于 2009-6-13 12:38 编辑 ]

rocwoods

help intersect

berryhaw

:handshake ,谢了,函数知道的还少,还得多看看

刚才调试了一下，intersect 效率不高，时间反而比原来的运行时间长了很多！！
还有没有其他办法呀？请高人指点！

[ 本帖最后由 ChaChing 于 2009-6-13 12:35 编辑 ]

rocwoods

楼主把你用intersect的代码写出来，我看看你怎么用的.
虽然没看到你的代码，但是看了你上面的那个代码后，你说intersect效率不如那个高，几乎一定是你没有用好intersect函数，或者是没有仔细看帮助文件。
我这里b是长200万，a长为100万的两个数组，用intersect时间1秒多钟。
用你上面那个代码，跑了1000多秒还没出来呢。（估算了下，估计需要1万六七千秒！），和intersect效率差万倍左右，尤其是a,b都很长的时候更明显。

还有楼主想想intersect函数专门干什么的，专门求交集的，你写的那段代码是求交集最容易想到的，效率很低的算法。不要轻易怀疑编写MATLAB函数大牛人的智商。呵呵:lol

[ 本帖最后由 rocwoods 于 2009-6-12 17:37 编辑 ]

berryhaw

首先感谢LS的建议,我的代码为：
b=[4 6 8]; a=[100 2 6 96 8 7 20];
[c,d]=intersect(b,a);
b(d)=0;
功能：找出交集后，还要将b中出现在交集的元素置零。
帮忙看一下，:handshake

补充一下，每次求交集的两个数组都不长，b的长度小于5，a的长度小于100；但这种运算要调用几百万次。

[ 本帖最后由 ChaChing 于 2009-6-13 12:31 编辑 ]

ChaChing

原帖由 berryhaw 于 2009-6-13 09:21 发表

                               
登录/注册后可看大图

...但这种运算要调用几百万次

怎感觉楼主问了不太好, 问题好像被误解了!?
几百万次当然耗时!?

berryhaw

看来是LS误解了，程序运行当然需要时间，现在的问题是：代码能否优化一下，以提高运行速度。

[ 本帖最后由 ChaChing 于 2010-5-24 11:03 编辑 ]

ChaChing

抱歉, 个人水平专业有限, 表达可能不够清楚!
直觉楼主的问题, 不在这一小部分, 而在那"几百万次"!? 也就是LZ的for loop太大了!
或许楼主可以说说原始问题, 看看有无改善空间!
试试下面程式
tic
b=[4 6 8]; a=[100 2 6 96 8 7 20];
for i=1:3, dd=a(a==b(i)); if isempty(dd)==0, b(i)=0; end; end
toc
tic
b=[4 6 8]; a=[100 2 6 96 8 7 20];
[c,d]=intersect(b,a); b(d)=0;
toc

[ 本帖最后由 ChaChing 于 2009-6-14 12:30 编辑 ]

berryhaw

感谢各位的关注及建议。
刚才发的截图看不到，我把代码及运行结果贴上：
试验一：
>> clear all
>> tic
b=[4 6 8]; a=[100 2 6 96 8 7 20];
for i=1:3, dd=a(a==b(i)); if isempty(dd)==0, b(i)=0; end; end
toc
Elapsed time is 0.009000 seconds.
>> clear all
>> tic
b=[4 6 8]; a=[100 2 6 96 8 7 20];
[c,d]=intersect(b,a);
b(d)=0;
toc
Elapsed time is 0.133000 seconds.
试验二：
>> tic
b=[4 6 8]; a=[100 2 6 96 8 7 20];
for i=1:3, dd=a(a==b(i)); if isempty(dd)==0, b(i)=0; end; end
toc
Elapsed time is 0.013000 seconds.
>> tic
b=[4 6 8]; a=[100 2 6 96 8 7 20];
[c,d]=intersect(b,a);
b(d)=0;
toc
Elapsed time is 0.127000 seconds.
>> tic
b=[4 6 8]; a=[100 2 6 96 8 7 20];
[c,d]=intersect(b,a);
b(d)=0;
toc
Elapsed time is 0.008000 seconds.
>>
首先，可以看出试验二的第三次执行时间是不准确的。
另外，从试验一看出，intersect的效率并不高，是不是我的调用出了问题？
肯请指点！！

rocwoods

楼主说的第三次执行时间是不准确的！恰恰相反，第三次应该是更准确的。准确判断一个函数的执行时间，应该是多次调用取平均值。
事实上是这样的，MATLAB中，函数调用都要有一定的开销，这种开销一般是built-in（像sin，cos，zeros，ones，find，all，any等等但凡你 type 文件名，MATLAB会告诉你XXX  is a built-in function.这样的函数）为最高，其次是一般的m文件，工具箱里能看到源码的那些m函数以及你自己写的m文件，匿名函数，子函数、nested function等等。这些总的来说调用效率都比较高，调用效率最低的是inline函数类型。
函数调用开销相对于数值计算来说往往相当可观，即使你用built-in，因为为了使得函数具有通用性，很多函数刚开始都会有一些判断、准备工作，往往判断好几层才到了真正算法实现部分。楼主之所以得到intersect函数效率不高，原因在于你举的例子规模太小了，intersect函数执行时间都花费在了函数调用上。而你用循环写的代码由于规模太小，效率低的事实被掩盖了。楼主试试下面代码就明白了：

1. 
   
2. clear;
   
3. a = unidrnd(10^9,1000000,1);
   
4. b = unidrnd(10^9,2000000,1);
   
5. tic;[c, ia, ib] = intersect(a, b);b(ib) = 0;toc
   

复制代码

以及

1. 
   
2. clear;
   
3. a = unidrnd(10^9,1000000,1);
   
4. b = unidrnd(10^9,2000000,1);
   
5. tic
   
6. for i=1:2000000
   
7. dd=a(a==b(i));
   
8. if isempty(dd)==0, b(i)=0; end
   
9. toc
   
10. end
    

复制代码

后者估计会让你等得不耐烦，从而ctrl+c中断执行。
所以，楼主应该想办法把那调用几百万次变成调用一次或者很少次数，每次判断的数组长些。这里提倡“少拿多取”，而不是自助餐提倡的“勤拿少取”

建议楼主把原始问题描述清楚。如果不容易变动的话，楼主或许可以用b(ismember(b,a)) = 0来解决

[ 本帖最后由 rocwoods 于 2009-6-14 15:50 编辑 ]

berryhaw

LS的说明一针见血，...intersect函数执行时间都花费在了函数调用上.......提倡“少拿多取”，而不是...“勤拿少取”.....
讨论了这么多，都怪我没把问题完整的拿出来，抱歉，请大家多包涵。 完整问题如下：
数组A是一个7*1000000的矩阵，每行的形式为：
A=[100 2 6 96 8 7 20;
15 69 7 6 4 20 11;
21 101 45 48 6 5 4;
...];
b 是一维数组[4 6 8];
程序功能是找出A中的每一行与数组b的并集，并将b中的并集元素置零，并将置零后的b存到B中。
B的形式为：
B=[4 0 0; 0 0 8; 0 0 8; ......]
主要代码为：
A=unidrnd(100,1000000,7);%这里先假设A是一个随机矩阵
B=zeros(1000000,3);
for m=1:1000000, a =A(m,:); b =[4 6 8];
   for i=1:3, dd=a(a==b(i));
      if isempty(dd)==0, b(i)=0; end
   end
   B(m,:)=b;
end

[ 本帖最后由 friendchj 于 2009-6-16 22:10 编辑 ]

rocwoods

采用少拿多取的思想可以如下优化：

1. 
   
2. clear
   
3. A=unidrnd(100,1000000,7);%这里先假设A是一个随机矩阵
   
4. AA = A;
   
5. B = repmat([4,6,8],1000000,1);
   
6. BB = B;
   
7. tic;C = [any(ismember(A,4),2) any(ismember(A,6),2) any(ismember(A,8),2)];B(C) = 0;toc
   

复制代码

或者

1. 
   
2. tic;C = [any(AA == 4,2) any(AA == 6,2) any(AA == 8,2)];
   
3. BB(C) = 0;toc
   
4. isequal(B,BB)
   

复制代码

前者比后者更快，大约快10%，可见ismember效率相当高，对于大规模数组，比单独恒等逻辑判断都高。（我用的MATLAB是R2009a版本）
前者方法在我电脑上用时0.18秒左右，而楼主用循环的那个用时8秒多，速度提高了40多倍！我的电脑CPU是  扣肉2 T8100，内存2G

[ 本帖最后由 rocwoods 于 2009-6-15 02:11 编辑 ]

ChaChing

还好直觉没钝掉! 问题真在那"几百万次"上!
还有11F的字体大小有点乱, 但不知怎麽回事总修不好, 其他版主有空试试!

berryhaw

还得请高手看看这段循环怎么矢量化,我弄了几天,没搞出来,拜托了….
问题和前面的问题差不多,如下:
A是一个100000*7的矩阵,每行的形式为:
A=[1004 12 10 4 1007 1005 2;6 1009 7 1005 1004 102 1003;1008 5 1007 2 4 1005 102;5 9 1004 2 1003 1 102;…]
b是一维数组[1004 1005 1007];
程序的功能是：
找出A中每一行介于“2到上一个小于1000的数”之间的数集(若无则为空)与数组b的并集,并将b中的并集元素置零；
找出A中每一行介于“102到上一个小于1000的数”之间的数集(若无则为空)与数组b的并集,并将b中的并集元素置零；
最后将置零后的b存到B中；
（注：A中标红的元素即为在相应B中要去除的元素）
B的形式为：[1004 0 0；0 0 1007；1004 0 0；0 1005 1007；…]
我的代码为:

tic
num=1000000;
A=unidrnd(1100,num,7);%这里先假设A是一个随机矩阵
% num=4;
% A=[1004 12 10 4 1007 1005 2;
% 6 1009 7 1005 1004 102 1003;
% 1008 5 1007 2 4 1005 102;
% 5 9 1004 2 1003 1 102];
A(A==102)=2;
B=zeros(num,3);
for m=1:num
    a=A(m,:);
    b=[1004 1005 1007];
    ea=find(a==2);%2(或102)在a中的位置
    eb=length(ea);%2(或102)在a中出现的次数
    ec=find(a<1000);%a中小于1000的元素的位置
    for i=1:3 %对b数组的每一元素进行遍历;
        %找出A中每一行介于“2(或102)到上一个小于1000的数”之间的数集
        %(若无则为空)与数组b的并集,并将b中的并集元素置零；
        for j=1:eb %对同一行中的每一个元素2进行遍历
            ed=find(ec==ea(j));%ea在ec中的位置
            if ed>1
                n=ec(ed-1)+1;%n为2(或102)前面一个小于1000的元素在a中的位置
            else
                n=1;
            end
            fa=[];
            if (n<=ea(j)-1)
                fa=a(a(n:ea(j)-1)==b(i));
            end
            if (isempty(fa)==0)
                b(i)=0;
            end
        end%j循环结束
        %------------------------------------------        
    end %i循环结束
    B(m,:)=b;
end %m循环结束
toc

berryhaw

版主合并了，可是时间好象不大对，我先顶起来，呵呵。
原的谢贴被删了，在这里还是要感谢大家的帮助！！