performance - una - vectores en matlab

vectorizar/optimizar este código en MATLAB?

Question

Estoy construyendo mi primer programa MATLAB a gran escala, y he logrado escribir código vectorizado original para todo, hasta que llegué a intentar crear una imagen que represente la densidad vectorial en proyección estereográfica. Después de un par de intentos fallidos, fui al sitio de intercambio de archivos de Mathworks y encontré un programa de código abierto que se ajusta a mis necesidades, cortesía de Malcolm Mclean . Con una matriz de prueba, su función produce algo como esto:

Y aunque esto es casi exactamente lo que quería, su código se basa en un bucle for triply anidado. En mi estación de trabajo, una matriz de datos de prueba de tamaño 25000x2 tardó 65 segundos en esta sección de código. Esto es inaceptable ya que voy a escalar a una matriz de datos de tamaño 500000x2 en mi proyecto.

Hasta ahora he podido vectorizar el bucle más interno (que era el bucle más largo / peor), pero me gustaría continuar y deshacerme de los bucles por completo si es posible. Aquí está el código original de Malcolm que necesito para vectorizar:

dmap = zeros(height, width); % height, width: scalar with default value = 32
for ii = 0: height - 1 % 32 iterations of this loop
yi = limits(3) + ii * deltay + deltay/2; % limits(3) & deltay: scalars
for jj = 0 : width - 1 % 32 iterations of this loop
xi = limits(1) + jj * deltax + deltax/2; % limits(1) & deltax: scalars
dd = 0;
for kk = 1: length(x) % up to 500,000 iterations in this loop
dist2 = (x(kk) - xi)^2 + (y(kk) - yi)^2;
dd = dd + 1 / ( dist2 + fudge); % fudge is a scalar
end
dmap(ii+1,jj+1) = dd;
end
end

Y aquí está con los cambios que ya hice en el ciclo más interno (que fue la mayor pérdida de eficiencia). Esto reduce el tiempo de 65 segundos a 12 segundos en mi máquina para la misma matriz de prueba, que es mejor pero mucho más lenta de lo que me gustaría.

dmap = zeros(height, width); for ii = 0: height - 1 yi = limits(3) + ii * deltay + deltay/2; for jj = 0 : width - 1 xi = limits(1) + jj * deltax + deltax/2; dist2 = (x - xi) .^ 2 + (y - yi) .^ 2; dmap(ii + 1, jj + 1) = sum(1 ./ (dist2 + fudge)); end end Entonces, mi principal pregunta es si hay más cambios que pueda hacer para optimizar este código. ¿O incluso un método alternativo para abordar el problema? He considerado usar C ++ o F # en lugar de MATLAB para esta sección del programa, y ​​puedo hacerlo si no puedo alcanzar un nivel de eficiencia razonable con el código MATLAB.

También tenga en cuenta que en este punto no tengo CUALQUIER caja de herramientas adicional, si lo hiciera, entonces sé que esto sería trivial (utilizando hist3 de la caja de herramientas de estadísticas, por ejemplo).

Answer 1

Solución que consume Mem

yi = limits(3) + deltay * ( 1:height ) - .5 * deltay;
xi = limits(1) + deltax * ( 1:width ) - .5 * deltax;
dx = bsxfun( @minus, x(:), xi ) .^ 2;
dy = bsxfun( @minus, y(:), yi ) .^ 2;
dist2 = bsxfun( @plus, permute( dy, [2 3 1] ), permute( dx, [3 2 1] ) );
dmap = sum( 1./(dist2 + fudge ) , 3 );

EDITAR

manejando extremadamente grandes y al dividir la operación en bloques:

blockSize = 50000; % process up to XX elements at once
dmap = 0;
yi = limits(3) + deltay * ( 1:height ) - .5 * deltay;
xi = limits(1) + deltax * ( 1:width ) - .5 * deltax;
bi = 1;
while bi <= numel(x)
% take a block of x and y
bx = x( bi:min(end, bi + blockSize - 1) );
by = y( bi:min(end, bi + blockSize - 1) );
dx = bsxfun( @minus, bx(:), xi ) .^ 2;
dy = bsxfun( @minus, by(:), yi ) .^ 2;
dist2 = bsxfun( @plus, permute( dy, [2 3 1] ), permute( dx, [3 2 1] ) );
dmap = dmap + sum( 1./(dist2 + fudge ) , 3 );
bi = bi + blockSize;
end 

Answer 2

Alternativamente, este problema se puede resolver usando técnicas de estimación de densidad de kernel . Esto es parte de la Caja de herramientas de estadísticas, o hay una implementación de KDE por Zdravko Botev (no se necesitan cajas de herramientas).

Para el siguiente código de ejemplo, obtengo 0.3 segundos para N = 500000, o 0.7 segundos para N = 1000000.

N = 500000;
data = [randn(N,2); rand(N,1)+3.5, randn(N,1);]; % 2 overlaid distrib
tic; [bandwidth,density,X,Y] = kde2d(data); toc;
imagesc(density);

Answer 3

Este es un buen ejemplo de por qué es importante comenzar un ciclo desde 1 . La única razón por la que ii y jj se inician en 0 es para eliminar los ii * deltay y jj * deltax que, sin embargo, introducen la secuencia en la indexación de dmap , evitando la paralelización .

Ahora, reescribiendo los bucles puedes usar parfor() después de abrir un matlabpool :

dmap = zeros(height, width);
yi = limits(3) + deltay*(1:height) - .5*deltay;
matlabpool 8
parfor ii = 1: height
for jj = 1: width
xi = limits(1) + (jj-1) * deltax + deltax/2;
dist2 = (x - xi) .^ 2 + (y - yi(ii)) .^ 2;
dmap(ii, jj) = sum(1 ./ (dist2 + fudge));
end
end
matlabpool close

Tenga en cuenta que abrir y cerrar el grupo tiene una sobrecarga significativa (10 segundos en mi Intel Core Duo T9300, vista 32 Matlab 2013a).

PD. No estoy seguro de si el bucle interno en lugar del externo puede ser significativamente paralelizado. Puede intentar cambiar el parámetro por el interno y comparar las velocidades (recomendaría ir inmediatamente a la matriz grande, ya que se está ejecutando en 12 segundos y la sobrecarga es casi igual de grande).