Matlab - 加速嵌套For循环

Question

一个简单的问题,但我对MATLAB并不是那么好.我有矢量x,(nx 1)y,(mx 1)和w = [x;y].我想将M(n + mx 1)定义为M(i)= x的元素数小于或等于w(i)(w被排序).这只是没有削减它:

N = n + m;
M = zeros(N,1);
for i = 1:N
  for j = 1:n
    if x(j) <= w(i)
      M(i) = M(i) + 1;
    end
  end
end

这不是一种特别聪明的方法,我的一些数据向量m和n大约是100000.

谢谢!

Answer 1

这可能看起来很神秘,但它应该给你与嵌套循环相同的结果:

M = histc(-x(:)',[-fliplr(w(:)') inf]);
M = cumsum(fliplr(M(1:N)))';

以上假设w已按升序排序.

说明

您的排序向量w可以被认为是用于创建具有HISTC函数的直方图的bin边.一旦计算了每个bin中(即边缘之间)的值的数量,使用CUMSUM函数对这些bin的累积总和将为您提供向量M.上面的代码看起来如此混乱(带有否定和函数FLIPLR)的原因是因为你想要找到x 小于或等于每个值的值w,但函数HISTC以下列方式存储数据:

n(k)计算值x(i)if edges(k) <= x(i) < edges(k+1).

请注意,小于每个箱的上限使用.您可能希望翻转行为,以便根据规则进行分类edges(k) < x(i) <= edges(k+1),这可以通过否定要分箱的值来实现,否定边缘,翻转边缘(因为边缘输入到HISTC必须单调非减少),然后翻转返回的bin计数.该值inf用作边值,用于计算小于w第一个bin中最低值的所有内容.

如果你想要找到的值x只是小于每个值w,代码就会简单得多:

M = histc(x(:)',[-inf w(:)']);
M = cumsum(M(1:N))';

Answer 2

内环至少可以替换为:

M(i)=sum(x<=w(i))

这将提供显着的性能改进.您可以考虑使用arrayfun:

M = arrayfun(@(wi)( sum( x <= wi ) ), w);

arrayfun不太可能在外部for循环上提供大量增益,但可能值得一试.

编辑:我应该注意,这个操作都没有w或x需要排序才能正常工作.

编辑:fwiw,我决定做一些实际的性能测试,所以我运行了这个程序:

n = 100000; m = n;

N = n + m;

x = rand(n, 1);
w = [x; rand(m, 1)];

tic;
M = zeros(N,1);
for i = 1:N
  for j = 1:n
    if x(j) <= w(i)
      M(i) = M(i) + 1;
    end
  end
end
perf = toc;
fprintf(1, 'Original : %4.3f sec\n', perf);

w = sort(w); % presorted, so don't incur time cost;
tic;
M = histc(-x(:)',[-fliplr(w(:)') inf]);
M = cumsum(fliplr(M(1:N)))';
perf = toc;
fprintf(1, 'gnovice : %4.3f sec\n', perf);

tic;
M = zeros(N,1);
for i = 1:N
    M(i)=sum(x<=w(i));
end
perf = toc;
fprintf(1, 'mine/loop : %4.3f sec\n', perf);

tic;
M = arrayfun(@(wi)( sum( x <= wi ) ), w);
perf = toc;
fprintf(1, 'mine/arrayfun : %4.3f sec\n', perf);

并获得n = 1000的这些结果:

Original : 0.086 sec
gnovice : 0.002 sec
mine/loop : 0.042 sec
mine/arrayfun : 0.070 sec

对于n = 100000:

Original : too long to tell ( >> 1m )
gnovice : 0.050 sec
mine/loop : too long to tell ( >> 1m )
mine/arrayfun : too long to tell ( >> 1m )