Matlab - ускорение вложенного цикла

Question

Простой вопрос, но я не так хорош с MATLAB. У меня есть векторы x, (n x 1) y, (m x 1) и w = [x;y]. Я хочу определить M (n + m x 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 . Как только вы посчитаете количество значений, попадающих в каждую ячейку (то есть между краями), совокупная сумма по этим ячейкам с помощью функции CUMSUM даст вам ваш вектор M. Причина, по которой вышеприведенный код выглядит настолько запутанным (с отрицаниями и функцией FLIPLR ) в том, что вы хотите найти значения в x , меньшие или равные каждому значению в w, но функция HISTC объединяет данные следующим образом:

n(k) считает значение x(i), если edges(k) <= x(i) < edges(k+1).

Обратите внимание, что меньше используется для верхнего предела каждой ячейки. Вы хотели бы перевернуть поведение так, чтобы вы складывали в соответствии с правилом edges(k) < x(i) <= edges(k+1), что может быть достигнуто путем отрицания значений, которые должны быть объединены, отрицания ребер, переворачивания ребер (так как ввод ребра в HISTC должно быть монотонно неубывающим), а затем переворачивать счетчик возвращаемых значений. Значение inf используется в качестве граничного значения для подсчета всего, что меньше минимального значения в w в первом бине.

Если вы хотите найти в 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, но, возможно, стоит попробовать.

edit: я должен отметить, что ни w, ни x не нужно сортировать для правильной работы этой операции.

edit: 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 )

Answer 3

Попробуйте вместо этого:

M = sum( bsxfun(@le, w', sort(w)) , 2 )

Answer 4

Некоторое время не делал MATLAB, но это должно сработать:

• Сортировка x с помощью встроенного алгоритма сортировки вверх.

• Используйте цикл с блуждающим индексом для итерации только один раз по x (j)

  j = 1;
  for i = 1:N
    while j <= n && x(j) <= w(i)
      M(i) = M(i) + 1;
      j = j+1;
    end
  end
  

• Наконец, накопить сумму

  for j =2:n
    M(j) = M(j-1) + M(j)
  end
  

Answer 5

У меня нет Matlab передо мной, поэтому я не могу подтвердить, что это на 100% функционально, но вы можете попробовать что-то вроде:

for i = 1:N
    M(i) = arrayfun(@(ary,val)length(find(ary <= val)), x, w(i))
end