Векторизация сумм различных диагоналей в матрице
Я хочу векторизовать следующий код MATLAB. Я думаю, что это должно быть просто, но я нахожу его путающий, тем не менее.
r = some constant less than m or n
[m,n] = size(C);
S = zeros(m-r,n-r);
for i=1:m-r+1
for j=1:n-r+1
S(i,j) = sum(diag(C(i:i+r-1,j:j+r-1)));
end
end
Код вычисляет таблицу очков, S, для алгоритма динамического программирования, от другой таблицы счета, C.
Диагональное подведение итогов должно генерировать очки к отдельным частям данных, используемых для генерации C для всех возможных частей (размера r).
Заранее спасибо за любые ответы! Извините, если этот должен быть очевидным...
Примечание:
Встроенный conv2 оказался быстрее, чем convnfft, потому что мой глаз (r) является довольно маленьким (5 <= r <= 20). convnfft.m указывает, что r должен быть> 20 для любого преимущества для декларации.
3 ответа
Если я правильно понимаю, вы пытаетесь вычислить диагональную сумму каждого подмассива C, где вы удалили последнюю строку и столбец C (если вы не должны удалять строку / столбец, вам нужно выполнить цикл, чтобы m-r + 1, и вам нужно передать весь массив C функции в моем решении ниже).
Вы можете выполнить эту операцию с помощью свертки , например:
S = conv2(C(1:end-1,1:end-1),eye(r),'valid');
Если C и r большие, вы можете посмотреть CONVNFFT из файла Matlab Обмен для ускорения расчетов.
Я думаю, вам может потребоваться преобразовать C в трехмерную матрицу, прежде чем суммировать ее по одному из измерений. Я отправлю ответ в ближайшее время.
РЕДАКТИРОВАТЬ
Мне не удалось найти способ аккуратно векторизовать его, но я нашел функцию аккумулятор , которая могла бы помочь. Я рассмотрю это более подробно, когда буду дома.
РЕДАКТИРОВАТЬ №2
Нашел более простое решение с использованием линейной индексации, но это может потребовать больших объемов памяти.
В C (1,1) индексы, которые мы хотим суммировать, равны 1+ [0, m + 1, 2 * m + 2, 3 * m + 3, 4 * m + 4, ...], или (0: r-1) + (0: m: (r-1) * m)
sum_ind = (0:r-1)+(0:m:(r-1)*m);
создать S_offset , матрицу (mr) на (nr) на r, так что S_offset ( :,:, 1) = 0, S_offset (:,:, 2) = m + 1, S_offset (:,:, 3) = 2 * m + 2 и т. Д.
S_offset = permute(repmat( sum_ind, [m-r, 1, n-r] ), [1, 3, 2]);
создать S_base , матрицу адресов базового массива, от которой будет вычисляться смещение.
S_base = reshape(1:m*n,[m n]);
S_base = repmat(S_base(1:m-r,1:n-r), [1, 1, r]);
Наконец, используйте S_base + S_offset для адресации значений C.
S = sum(C(S_base+S_offset), 3);
Вы, конечно, можете использовать bsxfun и другие методы, чтобы сделать его более эффективным; здесь я решил выложить это для наглядности. Я еще не тестировал это, чтобы увидеть, как он соотносится с методом двойного цикла; Сначала мне нужно отправиться домой на ужин!
Основываясь на идее JS, и как Jonas указал в комментариях, это можно сделать в две строки, используя IM2COL с некоторыми манипуляциями с массивами:
B = im2col(C, [r r], 'sliding');
S = reshape( sum(B(1:r+1:end,:)), size(C)-r+1 );
В основном B содержит элементы всех скользящих блоков размера r-by-r над матрицей C. Затем мы берем элементы на диагонали каждого из этих блоков B(1:r+1:end,:), вычисляем их сумму и перестраиваем результат под ожидаемый размер.
Сравнивая это с решением Джонаса на основе свертки, можно сказать, что здесь не выполняется никакого умножения матриц, только индексация...