Vectorizar sumas de diferentes diagonales en una matriz

Question

Quiero vectorizar el siguiente código MATLAB. Creo que debe ser simple, pero me resulta confuso, sin embargo.

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 

El código calcula una tabla de puntuaciones, S, para un algoritmo de programación dinámica, desde otra tabla de puntuación, C.
La suma diagonal es generar puntajes para piezas individuales de los datos utilizados para generar C, para todas las piezas posibles (de tamaño r).

¡Gracias de antemano por cualquier respuesta! Lo siento si éste debería ser obvio ...

Nota
El built-in CONV2 resultó ser más rápido que convnfft, porque mi ojo (r) es bastante pequeño (5 < = r = < 20) . convnfft.m establece que r debe ser> 20 para que se manifieste cualquier beneficio.

Answer 1

Si lo entiendo correctamente, está tratando de calcular la suma diagonal de cada subcampo de C, donde ha eliminado la última fila y columna de C (si no debe eliminar la fila/columna, debe realizar un ciclo a m-r + 1, y necesita pasar toda la matriz C a la función en mi solución a continuación).

Usted puede hacer esta operación a través de un convolution, así:

S = conv2(C(1:end-1,1:end-1),eye(r),'valid'); 

Si C y r son grandes, es posible que desee echar un vistazo a CONVNFFT del Matlab intercambio de archivos para acelerar los cálculos.

Answer 2

Creo que es posible que tenga que reorganizar C en una matriz 3D antes de sumarlo a lo largo de una de las dimensiones. Voy a publicar con una respuesta en breve.

EDITAR

No he conseguido encontrar una manera de Vectorise limpiamente, pero encontré la función accumarray, lo que podría ser de alguna ayuda. Lo veré con más detalle cuando esté en casa.

editar # 2

encontrado una solución más simple mediante el uso de indexación lineal, pero esto podría ser intensivo de la memoria.

En C (1,1), los índices que queremos sumar son 1+ [0, m + 1, 2 * m + 2, 3 * m + 3, 4 * m + 4, ...] , o (0: r-1) + (0: m: (r-1) * m)

sum_ind = (0:r-1)+(0:m:(r-1)*m); 

crear S_offset, un (mr) por (nr) por la matriz r, de manera que S_offset (: ,:, 1) = 0, S_offset (:,:, 2) = m + 1, S_offset (:,:, 3) = 2 * m + 2, y así sucesivamente.

S_offset = permute(repmat(sum_ind, [m-r, 1, n-r]), [1, 3, 2]); 

crear S_base, se calculará una matriz de direcciones de la matriz de base de la que el offset.

S_base = reshape(1:m*n,[m n]); 
S_base = repmat(S_base(1:m-r,1:n-r), [1, 1, r]); 

Por último, utilizar S_base+S_offset para hacer frente a los valores de C.

S = sum(C(S_base+S_offset), 3); 

Puede, por supuesto, utilizar bsxfun y otros métodos para hacerlo más eficiente; aquí elegí presentarlo para mayor claridad. Todavía tengo que comparar esto para ver cómo se compara con el método de bucle doble; ¡Primero tengo que ir a casa a cenar!

Answer 3

Sobre la base de la idea de JS, y como Jonas señalado en los comentarios, esto se puede hacer en dos líneas usando IM2COL con un poco de manipulación de matrices:

B = im2col(C, [r r], 'sliding'); 
S = reshape(sum(B(1:r+1:end,:)), size(C)-r+1); 

Básicamente B contiene los elementos de todos los bloques de deslizamiento de tamaño r-by-r sobre la matriz C. Luego tomamos los elementos en la diagonal de cada uno de estos bloques B(1:r+1:end,:), calculamos su suma y redimensionamos el resultado al tamaño esperado.

---

Comparando esto con la solución basada en convolución por Jonas, esto no lleva a cabo ninguna multiplicación de matrices, solamente la indexación ...

Answer 4

Es esto lo que está buscando? Esta función agrega las diagonales y las pone en un vector similar a cómo la función 'suma' suma todas las columnas en una matriz y las pone en un vector.

function [diagSum] = diagSumCalc(squareMatrix, LLUR0_ULLR1) 
% 
% Input: squareMatrix: A square matrix. 
%  LLUR0_ULLR1: LowerLeft to UpperRight addition = 0  
%      UpperLeft to LowerRight addition = 1 
% 
% Output: diagSum: A vector of the sum of the diagnols of the matrix. 
% 
% Example: 
% 
% >> squareMatrix = [1 2 3; 
%     4 5 6; 
%     7 8 9]; 
% 
% >> diagSum = diagSumCalc(squareMatrix, 0); 
% 
% diagSum = 
% 
%  1 6 15 14 9 
% 
% >> diagSum = diagSumCalc(squareMatrix, 1); 
% 
% diagSum = 
% 
%  7 12 15 8 3 
% 
% Written by M. Phillips 
% Oct. 16th, 2013 
% MIT Open Source Copywrite 
% Contact mphillips@hmc.edu fmi. 
% 

if (nargin < 2) 
    disp('Error on input. Needs two inputs.'); 
    return; 
end 

if (LLUR0_ULLR1 ~= 0 && LLUR0_ULLR1~= 1) 
    disp('Error on input. Only accepts 0 or 1 as input for second condition.'); 
    return; 
end 

[M, N] = size(squareMatrix); 

if (M ~= N) 
    disp('Error on input. Only accepts a square matrix as input.'); 
    return; 
end 

diagSum = zeros(1, M+N-1); 

if LLUR0_ULLR1 == 1 
    squareMatrix = rot90(squareMatrix, -1); 
end 

for i = 1:length(diagSum) 
    if i <= M 
     countUp = 1; 
     countDown = i; 
     while countDown ~= 0 
      diagSum(i) = squareMatrix(countUp, countDown) + diagSum(i); 
      countUp = countUp+1; 
      countDown = countDown-1; 
     end 
    end 
    if i > M 
     countUp = i-M+1; 
     countDown = M; 
     while countUp ~= M+1 
      diagSum(i) = squareMatrix(countUp, countDown) + diagSum(i); 
      countUp = countUp+1; 
      countDown = countDown-1; 
     end 
    end 
end 

Saludos